This repository has been archived by the owner on Oct 26, 2022. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathphotometry.py
1789 lines (1501 loc) · 69.5 KB
/
photometry.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#! /usr/bin/env python2
# -*- coding: utf-8 -*-
# Copyright (c) 2012 Victor Terron. All rights reserved.
# Institute of Astrophysics of Andalusia, IAA-CSIC
#
# This file is part of LEMON.
#
# LEMON is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
description = """
This module does aperture photometry on all the FITS images that it receives as
arguments. Astronomical objects are automatically detected, using SExtractor,
on the first image (which will be referred to from now on as the 'sources
image'), and then photometry is done for their celestial coordinates on the
rest of the images. The output is a LEMON database, which contains the
instrumental magnitudes (with 25 as zero point) and signal-to-noise ratio
for the stars detected in the sources image. Additionally, the instrumental
magnitudes of the stars in the sources image is also stored in the database,
in order to allow us to approximately estimate how bright each object is.
By default, the sizes of the aperture and sky annulus are determined by the
median FWHM of the images in each photometric filter, but different options
make it possible, among others, to directly specify those sizes in pixels, or
to use apertures and sky annuli that depend on the FWHM of each individual
image. The celestial coordinates of the objects on which to do photometry,
if known, can be listed in a text file, in this manner skipping the sources
detection step and working exclusively with the specified objects.
"""
import atexit
import collections
import hashlib
import itertools
import logging
import multiprocessing
import numpy
import optparse
import os
import os.path
import pwd
import shutil
import socket
import sys
import tempfile
import time
import warnings
# LEMON modules
import util
import astromatic
import customparser
import database
import defaults
import fitsimage
import json_parse
import keywords
import qphot
import seeing
import style
# Message of the warning that is issued when the width of the sky annulus is
# smaller than the number of pixels specified with the --min-sky option.
DANNULUS_TOO_THIN_MSG = (
"Whoops! Sky annulus too thin, setting it to the minimum of %.2f pixels"
)
# The Queue is global -- this works, but note that we could have
# passed its reference to the function managed by pool.map_async.
# See http://stackoverflow.com/a/3217427/184363
queue = util.Queue()
def get_fwhm(img, options):
"""Return the FWHM of the FITS image.
Attempt to read the full width at half maximum from the header of the FITS
image (keyword options.fwhmk). If the keyword cannot be found, then compute
the FWHM by calling FITSeeingImage.fwhm(). In this manner, we can always
call this method to get the FWHM of each image, without having to worry
about whether it is in the header already. The 'img' argument must be a
fitsimage.FITSImage object, while 'options' must be the optparse.Values
object returned by optparse.OptionParser.parse_args().
"""
try:
msg = "%s: reading FWHM from keyword '%s'"
args = img.path, options.fwhmk
logging.debug(msg % args)
fwhm = img.read_keyword(options.fwhmk)
msg = "%s: FWHM = %.3f (keyword '%s')"
args = img.path, fwhm, options.fwhmk
logging.debug(msg % args)
return fwhm
except KeyError:
msg = "%s: keyword '%s' not found in header"
args = img.path, options.fwhmk
logging.debug(msg % args)
if not isinstance(img, seeing.FITSeeingImage):
msg = "%s: type of argument 'img' is not FITSeeingImage ('%s')"
args = img.path, type(img)
logging.debug(msg % args)
msg = "%s: calling FITSeeingImage.__init__() with 'img'"
logging.debug(msg % img.path)
args = (img.path, options.maximum, options.margin)
kwargs = dict(coaddk=options.coaddk)
img = seeing.FITSeeingImage(*args, **kwargs)
msg = "%s: calling FITSeeingImage.fwhm() to compute FWHM"
logging.debug(msg % img.path)
mode = "mean" if options.mean else "median"
kwargs = dict(per=options.per, mode=mode)
fwhm = img.fwhm(**kwargs)
msg = "%s: FITSeeingImage.fwhm() returned %.3f"
args = img.path, fwhm
logging.debug(msg % args)
return fwhm
@util.print_exception_traceback
def parallel_photometry(args):
"""Function argument of map_async() to do photometry in parallel.
This will be the first argument passed to multiprocessing.Pool.map_async(),
which chops the iterable into a number of chunks that are submitted to the
process pool as separate tasks. 'args' must be a three-element tuple with
(1) a fitsimage.FITSImage object, (2) a database.PhotometricParameters
object and (3) 'options', the optparse.Values object returned by
optparse.OptionParser.parse_args().
This function does photometry (qphot.run()) on the astronomical objects of
the FITS image listed in options.coordinates, using the aperture, annulus
and dannulus defined by the PhotometricParameters object. The result is
another three-element tuple, which is put into the module-level 'queue'
object, a process shared queue. This tuple contains (1) a database.Image
object, (2) a database.PhotometricParameters object and (3) a qphot.QPhot
object -- therefore mapping each FITS file and the parameters used for
photometry to the measurements returned by qphot.
"""
image, pparams, options = args
logging.debug("Doing photometry on %s" % image.path)
msg = "%s: qphot aperture: %.3f"
logging.debug(msg % (image.path, pparams.aperture))
msg = "%s: qphot annulus: %.3f"
logging.debug(msg % (image.path, pparams.annulus))
msg = "%s: qphot dannulus: %.3f"
logging.debug(msg % (image.path, pparams.dannulus))
maximum = image.saturation(options.maximum, coaddk=options.coaddk)
msg = "%s: saturation level = %d ADUs'"
args = (image.path, maximum)
logging.debug(msg % args)
logging.info("Running qphot on %s" % image.path)
args = (
image,
options.coordinates,
options.epoch,
pparams.aperture,
pparams.annulus,
pparams.dannulus,
maximum,
options.datek,
options.timek,
options.exptimek,
options.uncimgk,
)
img_qphot = qphot.run(*args, cbox=options.cbox)
logging.info("Finished running qphot on %s" % image.path)
msg = "%s: qphot.run() returned %d records"
args = (image.path, len(img_qphot))
logging.debug(msg % args)
pfilter = image.pfilter(options.filterk)
logging.debug("%s: filter = %s" % (image.path, pfilter))
kwargs = dict(
date_keyword=options.datek,
time_keyword=options.timek,
exp_keyword=options.exptimek,
)
unix_time = image.date(**kwargs)
msg = "%s: observation date: %.2f (%s)"
args = (image.path, unix_time, util.utctime(unix_time))
logging.debug(msg % args)
object_ = image.read_keyword(options.objectk)
logging.debug("%s: object = %s" % (image.path, object_))
airmass = image.read_keyword(options.airmassk)
logging.debug("%s: airmass = %.4f" % (image.path, airmass))
# If not given with --gaink, read it from the FITS header
gain = options.gain or image.read_keyword(options.gaink)
msg = "%s: gain (%s) = %.4f"
gain_msg = "given by user" if options.gain else "read from header"
args = image.path, gain_msg, gain
logging.debug(msg % args)
msg = "%s: calculating coordinates of field center"
logging.debug(msg % image.path)
ra, dec = image.center_wcs()
logging.debug("%s: RA (field center) = %.8f" % (image.path, ra))
logging.debug("%s: DEC (field center) = %.8f" % (image.path, dec))
args = (image.path, pfilter, unix_time, object_, airmass, gain, ra, dec)
db_image = database.Image(*args)
queue.put((db_image, pparams, img_qphot))
msg = "%s: photometry result put into global queue"
logging.debug(msg % image.path)
parser = customparser.get_parser(description)
parser.usage = "%prog [OPTION]... SOURCES_IMG INPUT_IMGS... OUTPUT_DB"
parser.add_option(
"--overwrite",
action="store_true",
dest="overwrite",
help="overwrite output database if it already exists",
)
parser.add_option(
"--filter",
action="append",
type="passband",
dest="filters",
default=None,
help="do not do photometry on all the FITS files given "
"as input, but only on those taken in this photometric "
"filter. This option may be used multiple times in order "
"to specify more than one filter with which we want to "
"work. " + defaults.desc["filter"],
)
parser.add_option(
"--exclude",
action="append",
type="passband",
dest="excluded_filters",
default=None,
help="ignore those FITS files taken in this photometric "
"filter. This option is the opposite of --filter, and may "
"be as well used multiple times in order to specify more "
"than one photometric filter that must be discarded.",
)
parser.add_option(
"--cbox",
action="store",
type="float",
dest="cbox",
default=5,
help="the width, in pixels, of the centering box to "
"search for the accurate center of the astronomical object "
"around the input coordinates. Distortions introduced in "
"the FoV by the optics, limitations in the precision of the "
"astrometry or proper movements of objects can contribute "
"to an apparent movement of the center of the object. In "
"these cases, leaving a few pixels' margin to search for "
"the center of the star improves the photometry. If you "
"absolutely trust the celestial coordinates of each object "
"and want photometry to be done without any centering, you "
"may set this option to zero [default: %default]",
)
parser.add_option(
"--maximum",
action="store",
type="int",
dest="maximum",
default=defaults.maximum,
help=defaults.desc["maximum"],
)
parser.add_option(
"--margin",
action="store",
type="int",
dest="margin",
default=defaults.margin,
help=defaults.desc["margin"],
)
parser.add_option(
"--gain",
action="store",
type="float",
dest="gain",
default=None,
help="the gain of the CCD, in e-/ADU. Needed in order to "
"accurately calculate the SNR of each measurement. In case "
"this option is given, the value will not be read from the "
"FITS header (--gaink option)",
)
parser.add_option(
"--annuli",
action="store",
type=str,
dest="json_annuli",
default=None,
help="ignore the Aperture Photometry (FWHM and Pixels) "
"sections below an instead read the apertures and sky annuli "
"to use from a JSON file output by the 'annuli' command. "
"This file should, of course, have been generated for the "
"same set of images on which photometry is now being done. "
"For this same reason, the execution will be aborted if the "
"JSON file does not have information for the photometric "
"filters of all the input FITS files. Even when this option "
"is used, the aperture and sky annulus used for the sources "
"image are determined by the 'Aperture Photometry' sections.",
)
parser.add_option(
"--cores",
action="store",
type="int",
dest="ncores",
default=defaults.ncores,
help=defaults.desc["ncores"],
)
parser.add_option(
"-v",
"--verbose",
action="count",
dest="verbose",
default=defaults.verbosity,
help=defaults.desc["verbosity"],
)
coords_group = optparse.OptionGroup(
parser,
"List of Coordinates",
"By default, we run SExtractor on the sources image in order "
"to detect the astronomical objects on which to do photometry. "
"Alternatively, with this option it is possible to skip the "
"detection stage and directly do photometry on the objects "
"whose celestial coordinates are specified in a text file. "
"Note that this option is exclusive with the automatic "
"detection of sources, so photometry will be done *only* on "
"these coordinates.",
)
coords_group.add_option(
"--coordinates",
action="store",
type=str,
dest="coordinates",
default=None,
help="path to the file containing the celestial "
"coordinates of the objects on which photometry must "
"be done. The file must have two columns, with the "
"right ascension and declination (in decimal degrees) "
"and, optionally, two other columns with the proper "
"motion in right ascension and declination (in "
"seconds of arc per year) surrounded by brackets. The "
"coordinates of the objects with a proper motion are "
"automatically corrected for each image, in order to "
"account for the change in their position over time.",
)
coords_group.add_option(
"--epoch",
action="store",
type=int,
dest="epoch",
default=2000,
help="epoch of the coordinates [default: %default]",
)
parser.add_option_group(coords_group)
qphot_group = optparse.OptionGroup(
parser,
"Aperture Photometry (FWHM)",
"Photometry is done on the images by IRAF's qphot, the quick "
"aperture photometer, which computes accurate sky values and "
"magnitudes for a series of objects. Instead of using "
"absolute values (such as, for example, 8, 11 or 13 pixels), "
"the values of the following parameters are defined in terms "
"of the *median* FWHM of the images in each filter.\n\n"
"The full width at half-maximum of an image is defined here as "
"the median of the FWHM of all the astronomical objects detected "
"by SExtractor. Therefore, you may think of the median 'FWHM of "
"the images' as the 'median of the median of the FWHM of all the "
"astronomical objects, if you find it easier. Note also that a "
"different median FWHM is computed for each photometric filter "
"on which photometry is done.",
)
qphot_group.add_option(
"--aperture",
action="store",
type="float",
dest="aperture",
default=3.0,
help="the aperture radius, in number of times the "
"median FWHM [default: %default]",
)
qphot_group.add_option(
"--annulus",
action="store",
type="float",
dest="annulus",
default=4.5,
help="the inner radius of the sky annulus, in "
"number of times the median FWHM [default: %default]",
)
qphot_group.add_option(
"--dannulus",
action="store",
type="float",
dest="dannulus",
default=1.0,
help="the width of the sky annulus, in number "
"of times the median FWHM [default: %default]",
)
qphot_group.add_option(
"--min-sky",
action="store",
type="float",
dest="min",
default=3.0,
help="the minimum width of the sky annulus, in "
"pixels, regardless of the value specified in the "
"above parameter. This option is intended to prevent "
"small FWHMs from resulting in too thin a sky "
"annulus. [default = %default]",
)
qphot_group.add_option(
"--individual-fwhm",
action="store_true",
dest="individual_fwhm",
help="consider FITS images individually when the "
"FWHM (and, therefore, the derived aperture and sky "
"annuli) is computed. That is, instead of using the "
"median FWHM of *all* the images in the same filter, "
"when this option is present the aperture and sky "
"annuli used to do photometry on an image depend "
"exclusively on the FWHM of *that* image.",
)
parser.add_option_group(qphot_group)
qphot_fixed = optparse.OptionGroup(
parser,
"Aperture Photometry (pixels)",
"In case the exact sizes of the aperture and sky annulus are "
"known, their dimensions can be specified in pixels. If used, "
"these three options must be used together.\n\n"
"Note that setting the aperture photometry parameters to a "
"fixed size, as these options do, means that the same values "
"are used regardless of the photometric filter of the images, "
"including the sources image. You probably want to use these "
"options in conjunction with --filter, in order to do "
"photometry one photometric filter at a time.",
)
qphot_fixed.add_option(
"--aperture-pix",
action="store",
type="float",
dest="aperture_pix",
default=None,
help="the aperture radius, in pixels",
)
qphot_fixed.add_option(
"--annulus-pix",
action="store",
type="float",
dest="annulus_pix",
default=None,
help="the inner radius of the sky annulus, in pixels",
)
qphot_fixed.add_option(
"--dannulus-pix",
action="store",
type="float",
dest="dannulus_pix",
default=None,
help="the width of the sky annulus, in pixels",
)
parser.add_option_group(qphot_fixed)
fwhm_group = optparse.OptionGroup(
parser,
"FWHM",
"The full width at half maximum of each FITS image is written "
"to its header by the 'seeing' command, which is generally run "
"before photometry is done. However, it may be the case that "
"not all the images with which we have to work now are the "
"output of 'seeing' (for example, if the sources image is that "
"output by the 'mosaic' command). In these cases, we need to "
"compute the FWHM of these images, in the same way the "
"'seeing' command does.",
)
fwhm_group.add_option(
"--snr-percentile",
action="store",
type="float",
dest="per",
default=defaults.snr_percentile,
help=defaults.desc["snr_percentile"],
)
fwhm_group.add_option(
"--mean", action="store_true", dest="mean", help=defaults.desc["mean"]
)
parser.add_option_group(fwhm_group)
key_group = optparse.OptionGroup(parser, "FITS Keywords", keywords.group_description)
key_group.add_option(
"--objectk",
action="store",
type="str",
dest="objectk",
default=keywords.objectk,
help=keywords.desc["objectk"],
)
key_group.add_option(
"--filterk",
action="store",
type="str",
dest="filterk",
default=keywords.filterk,
help=keywords.desc["filterk"],
)
key_group.add_option(
"--datek",
action="store",
type="str",
dest="datek",
default=keywords.datek,
help=keywords.desc["datek"],
)
key_group.add_option(
"--timek",
action="store",
type="str",
dest="timek",
default=keywords.timek,
help=keywords.desc["timek"],
)
key_group.add_option(
"--expk",
action="store",
type="str",
dest="exptimek",
default=keywords.exptimek,
help=keywords.desc["exptimek"],
)
key_group.add_option(
"--coaddk",
action="store",
type="str",
dest="coaddk",
default=keywords.coaddk,
help=keywords.desc["coaddk"],
)
key_group.add_option(
"--gaink",
action="store",
type="str",
dest="gaink",
default=keywords.gaink,
help=keywords.desc["gaink"],
)
key_group.add_option(
"--fwhmk",
action="store",
type="str",
dest="fwhmk",
default=keywords.fwhmk,
help=keywords.desc["fwhmk"],
)
key_group.add_option(
"--airmk",
action="store",
type="str",
dest="airmassk",
default=keywords.airmassk,
help=keywords.desc["airmassk"],
)
key_group.add_option(
"--uik",
action="store",
type="str",
dest="uncimgk",
default=keywords.uncimgk,
help=keywords.desc["uncimgk"],
)
parser.add_option_group(key_group)
customparser.clear_metavars(parser)
def main(arguments=None):
"""main() function, encapsulated in a method to allow for easy invokation.
This method follows Guido van Rossum's suggestions on how to write Python
main() functions in order to make them more flexible. By encapsulating the
main code of the script in a function and making it take an optional
argument the script can be called not only from other modules, but also
from the interactive Python prompt.
Guido van van Rossum - Python main() functions:
http://www.artima.com/weblogs/viewpost.jsp?thread=4829
Keyword arguments:
arguments - the list of command line arguments passed to the script.
"""
if arguments is None:
arguments = sys.argv[1:] # ignore argv[0], the script name
# All the parameters must be cast to string before beging passed to the
# option parser, as apparently it only expects str data types. Otherwise,
# optparse may raise a TypeError exception and complains about how int or
# float objects are unsuscriptable.
arguments = [str(param) for param in arguments]
(options, args) = parser.parse_args(args=arguments)
# Adjust the logger level to WARNING, INFO or DEBUG, depending on the
# given number of -v options (none, one or two or more, respectively)
logging_level = logging.WARNING
if options.verbose == 1:
logging_level = logging.INFO
elif options.verbose >= 2:
logging_level = logging.DEBUG
logging.basicConfig(format=style.LOG_FORMAT, level=logging_level)
# Print the help and abort the execution if there are not three positional
# arguments left after parsing the options, as the user must specify the
# sources image, at least one (only one?) image on which to do photometry
# and the output LEMON database.
if len(args) < 3:
parser.print_help()
return 2 # 2 is generally used for command line syntax errors
else:
sources_img_path = args[0]
input_paths = set(args[1:-1])
output_db_path = args[-1]
assert input_paths
# If the user gives an empty string as the FITS keyword which stores the
# path to the original image, it is understood as meaning that we want
# saturation to be checked for in the same images on which photometry is
# done. If that is the case, we need to set the option to None (although an
# empty string would also work), as that is what qphot.run expects to
# receive in these cases.
if not options.uncimgk:
options.uncimgk = None
# There is not much point in using --filter and --exclude at the same time:
# the former discards all the images not taken in one or more filters, and
# the latter discards the images taken in one or more filters. We should
# always need to use one approach or the other, but never both at once.
if options.filters and options.excluded_filters:
msg = "%sError. The --filter and --exclude options are incompatible."
print msg % style.prefix
print style.error_exit_message
return 1
# The annuli JSON file (that generated by the annuli command, and specified
# with the --annuli option) must exist. The use of this file automatically
# discards whathever was specified with the Aperture Photometry (FWHM and
# pixels) options.
json_annuli = None
fixed_annuli = False
if options.json_annuli:
if not os.path.exists(options.json_annuli):
print "%sError. The file '%s' does not exist." % (
style.prefix,
options.json_annuli,
)
print style.error_exit_message
return 1
else:
json_annuli = json_parse.CandidateAnnuli.load(options.json_annuli)
msg = "%sPhotometric parameters read from the '%s' file."
args = (style.prefix, os.path.basename(options.json_annuli))
print msg % args
# Even if the annuli JSON file is used, the aperture and sky annuli
# parameters (whether FWHM-based or given in pixels) are still needed in
# order to do photometry on the sources image and extract the instrumental
# magnitude that for each astronomical object is stored in the database.
# Abort the execution if one or more of the {aperture,{,d}annulus}-pix:
# options are given, but not all three. If we are going to use fixed sizes
# for the aperture and sky annuli we need to know all of them.
fixed_pix_count = (
bool(options.aperture_pix)
+ bool(options.annulus_pix)
+ bool(options.dannulus_pix)
)
if fixed_pix_count:
if fixed_pix_count < 3:
assert 1 <= fixed_pix_count <= 2
msg = (
"%sError. The --aperture-pix, --annulus-pix and "
+ "--dannulus-pix options must be used together."
)
print msg % style.prefix
print style.error_exit_message
return 1
else:
assert fixed_pix_count == 3
fixed_annuli = True
# Abort the execution if the user gives, at the same time, the
# --aperture-pix, --annulus-pix, --dannulus-pix and --annuli options. It
# does not make sense to set the aperture and sky annuli to a fixed value
# and simultaneously specify that these very values must be read from the
# JSON file. Does not compute.
if fixed_annuli and json_annuli:
print (
"%sError. The --aperture-pix, --annulus-pix and --dannulus-pix "
"options are incompatible with --annuli." % style.prefix
)
print style.error_exit_message
return 1
if options.individual_fwhm:
# If the photometric parameters are set to a fixed value, they cannot
# be also derived from the FWHM of each image.
if fixed_annuli:
print "%sError. The --aperture-pix, --annulus-pix and " "--dannulus-pix options are incompatible with " "--individual-fwhm." % style.prefix
print style.error_exit_message
return 1
# The same applies to --annuli: if the photometric parameters are read
# from the JSON file, they cannot also depend on the FWHM of the images.
if json_annuli:
print "%sError. The --annuli option is incompatible with " "--individual-fwhm." % style.prefix
print style.error_exit_message
return 1
# The aperture, annulus and dannulus values, whether expressed in number of
# times the median FWHM or by a fixed number of pixels, must be positive
# numbers. By definition, also, the inner radius of the sky annulus must be
# greater than or equal to the aperture radius. Obviously!
fwhm_options = (options.aperture, options.annulus, options.dannulus)
pixel_options = (options.aperture_pix, options.annulus_pix, options.dannulus_pix)
if (not fixed_annuli and min(fwhm_options) <= 0) or (
fixed_annuli and min(pixel_options) <= 0
):
print "%sError. The aperture, annulus and dannulus values must be " "positive numbers." % style.prefix
print style.error_exit_message
return 1
if (not fixed_annuli and options.aperture > options.annulus) or (
fixed_annuli and options.aperture_pix > options.annulus_pix
):
print "%sError. The aperture radius (%.2f) must be smaller than or equal\n" "%sto the inner radius of the sky annulus (%.2f)" % (
style.prefix,
options.aperture_pix if fixed_annuli else options.aperture,
style.prefix,
options.annulus_pix if fixed_annuli else options.annulus,
)
print style.error_exit_message
return 1
# If the --coordinates option has been given, read the text file and store
# the four-element tuples (right ascension, declination and proper motions)
# in a list, as astromatic.Coordinates objects. Abort the execution if the
# coordinates file is empty.
if options.coordinates:
sources_coordinates = []
for args in util.load_coordinates(options.coordinates):
coords = astromatic.Coordinates(*args)
sources_coordinates.append(coords)
if not sources_coordinates:
msg = "%sError. Coordinates file '%s' is empty."
print msg % (style.prefix, options.coordinates)
print style.error_exit_message
return 1
# Each campaign must be saved to its own LEMON database, as it would not
# make much sense to merge data (since the same tables would be used) of
# astronomical objects that belong to different fields. Thus, we refuse to
# work with an existing database (which is what the LEMONdB class would do
# otherwise) unless the --overwrite option is given, in which case it is
# deleted and created again from scratch.
if os.path.exists(output_db_path):
if not options.overwrite:
print "%sError. The output database '%s' already exists." % (
style.prefix,
output_db_path,
)
print style.error_exit_message
return 1
else:
os.unlink(output_db_path)
# Loop over all the input FITS files, mapping (a) each photometric filter
# to a list of the FITS images that were observed in it, and (b) each FITS
# image to its date of observation (UTC), in Unix time.
msg = "%sExamining the headers of the %s FITS files given as input..."
print msg % (style.prefix, len(input_paths))
def get_date(img):
""" Return the date() of a FITSImage object """
return img.date(
date_keyword=options.datek,
time_keyword=options.timek,
exp_keyword=options.exptimek,
)
files = fitsimage.InputFITSFiles()
img_dates = {}
util.show_progress(0.0)
for index, img_path in enumerate(input_paths):
img = fitsimage.FITSImage(img_path)
pfilter = img.pfilter(options.filterk)
files[pfilter].append(img_path)
date = get_date(img)
img_dates[img_path] = date
percentage = (index + 1) / len(input_paths) * 100
util.show_progress(percentage)
print # progress bar doesn't include newline
msg = "%s%d different photometric filters were detected:"
print msg % (style.prefix, len(files.keys()))
for pfilter, images in sorted(files.iteritems()):
msg = "%s %s: %d files (%.2f %%)"
percentage = len(images) / len(files) * 100
print msg % (style.prefix, pfilter, len(images), percentage)
# Light curves, which are our ultimate goal, can only have one magnitude
# for each point in time. Therefore, we cannot do photometry on two or more
# images with the same observation date and photometric filter. This may
# seem (and, indeed, is) unlikely to happen, but astronomical instruments
# also have software errors - we have already come across this while
# reducing images taken with Omega 2000, a camera for the 3.5m CAHA.
#
# We might be tempted to keep one of them, such as, for example, that with
# the highest number of sources, or the best FWHM. However, the safest bet
# is to discard them all, because we cannot know which image has the right
# observation date. The only certain thing is that an error occurred. Thus,
# we better forget about these images.
msg = "%sMaking sure there are no images with the same date and filter..."
print msg % style.prefix,
sys.stdout.flush()
# Two-level dictionary: map each date of observation (UTC), in Unix time,
# to a photometric filter to a list of the corresponding FITS files. If
# there are no two or more images with the same date and filter, all the
# second-level values of the dictionary will have a length of one. We do
# not use the Counter class, from the collections module, for Python 2.6
# compatibility.
get_dict = lambda: collections.defaultdict(list)
dates_counter = collections.defaultdict(get_dict)
# 'files' maps each filter to a list of images, while 'img_dates' maps each
# image to its Unix date. There is no need, therefore, to read the date or
# photometric filter of the images from their FITS headers again.
for pfilter, images in files.iteritems():
for img_path in images:
date = img_dates[img_path]
dates_counter[date][pfilter].append(img_path)
# Find the dates and filters for which there is more than one FITS file.
# Then, remove from the InputFITSFiles object each of these images with
# a duplicate Unix time and photometric filter.
discarded = 0
for date, date_images in dates_counter.items():
for pfilter, images in date_images.items():
if len(images) > 1:
# "Making sure..." message above does not include newline.
# We need to print it, but only for the first issued warning.
if not discarded:
print
msg = "%sWarning! Multiple images have date %s and filter %s"
args = style.prefix, util.utctime(date), pfilter
warnings.warn(msg % args)
del dates_counter[date][pfilter]
for img in images:
discarded += files.remove(img)
if not discarded:
print "done."
else:
# There should be no FITS files with the same observation date and
# filter anymore, and at least two of them should have been discarded.
# Otherwise, how did we get to the 'else' clause in the first place?
if __debug__:
dates = []
for image in files:
dates.append(img_dates[image])
assert len(set(dates)) == len(dates)
assert discarded >= 2
msg = "%s%d images had duplicate dates and were discarded, %d remain."
print msg % (style.prefix, discarded, len(files))
# The --filter option allows the user to specify on which FITS files, among
# all those received as input, photometry must be done: only those files in
# any of the photometric filters contained in options.filter. The Passband
# class, which supports comparison operations, makes it possible to compare
# filters for what they really are, not how they were written: "Johnson V"
# and "johnson_v", for example, are the same filter after all, but if we
# just compared the two strings we would consider them to be different.
if options.filters:
msg = "%sIgnoring images not taken in any of the following filters:"
print msg % style.prefix
for index, pfilter in enumerate(sorted(options.filters)):
print "%s (%d) %s" % (style.prefix, index + 1, pfilter)
sys.stdout.flush()
discarded = 0
for pfilter, images in files.items():
if pfilter not in options.filters:
discarded += len(images)
del files[pfilter]
if not files:
print
msg = "%sError. No image was taken in any of the above filters."
print msg % style.prefix
print style.error_exit_message
return 1