-
Notifications
You must be signed in to change notification settings - Fork 3
/
podi_photcalib.py
executable file
·2228 lines (1740 loc) · 88.2 KB
/
podi_photcalib.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 python3
#
# Copyright 2012-2013 Ralf Kotulla
# kotulla@uwm.edu
#
# This file is part of the ODI QuickReduce pipeline package.
#
# If you find this program or parts thereof please make sure to
# cite it appropriately (please contact the author for the most
# up-to-date reference to use). Also if you find any problems
# or have suggestiosn on how to improve the code or its
# functionality please let me know. Comments and questions are
# always welcome.
#
# The code is made publicly available. Feel free to share the link
# with whoever might be interested. However, I do ask you to not
# publish additional copies on your own website or other sources.
#
# 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.
#
"""
This module contains all functionality for performing a photometric calibration
for a given frame. It does so by matching the ODI source catalog to a catalog
of stars from the SDSS. The comparison between instrumental and reference
magnitudes then yields the calibration zeropoint.
If requested, podi_photcalib also creates diagnostic plots that allow to easily
judge the quality of the obtained calibration.
The module also contains all required functionality to access local SDSS star
catalogs or, if no local copies are available, query the SDSS online.
Command line options
--------------------
* **-multi**
Deprecated for now, do not use
* **-new**
* **-querysdss**
Test-routine that queries the SDSS and prints the returned catalog
Run as:
``podi_photcalib -querysdss ra_min ra_max dec_min dec_max (-print)``
If the -print option is given as well print the entire catalog to stdout,
otherwise only print the number of returned stars.
* **-swarp**
Compute the photometric zeropoint for a frame created with swarp.
Run as:
''podi_photcalib.py -swarp swarped_frame.fits swarp_input_frame.fits``
swarp_input_frame here is one of the frames that went into the swarped frame.
This is necessary to obtain the FILTER and EXPTIME keywords that are required
for the photometric calibration, but are typically not contained in the
swarped frame any longer.
* **-aucap**
Compute the photometric zeropoint for a frame created with the Automatic
Calibration Pipeline (AuCaP).
Run as:
``podi_photcalib.py -aucap file1.fits file2.fits file3.fits ...``
In this mode, podi_photcalib can run on any number of input frames in sequence.
Note for the -aucap and -swarp modi
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
To save time during programming, two additional flags are supported in these
modi. By default, if the source catalog for the specified input file already
exists, it is simply re-used and not created. To force re-running SExtractor,
add the -resex flag. Additionally, in the ``-aucap`` mode, if the source catalog
exists, podi_photcalib also does not re-create the diagnostic plots. To change
this behavior and re-create the plots, add the ``-replot`` flag.
Methods
-------
"""
from __future__ import print_function
import sys
import numpy
import os
from query_usno import query_usno
from podi_definitions import *
from podi_commandline import *
import astropy.io.fits as pyfits
#import date
import datetime
from astLib import astWCS
import pdb
import scipy
import scipy.stats
from wiyn_filters import *
import podi_matchcatalogs
import podi_sitesetup as sitesetup
import podi_search_ipprefcat
import podi_logging
import logging
import time
import subprocess
from podi_definitions import reference_zeropoint
arcsec = 1./3600.
number_bright_stars = 100
max_offset = 0.1
N = 200
N_brightest_ota = 70
N_brightest_ref = 150
# Compute matches in smaller blocks to not blow up memory
# Improve: Change execution to parallel !!!
blocksize = 100
def load_catalog_from_stripe82cat(ra, dec, calib_directory, sdss_filter):
"""
Load the source catalog from the custom pODI Stripe82 catalog.
Parameters
----------
ra, dec -- float
Center position of the search area. The search radius is set to 0.6 deg.
calib_directory - string
directory that contains the Stripe82 catalog
sdss_filter - string
Name of the filter for which the magnitudes should be returned.
Returns
-------
A catalog of sources within the search area, stored in an array. Columns
are as follows. Ra, Dec, mag_median, mag_mean, mag_std, mag_var
"""
# Load the standard star catalog
sdss_cat_filename = "%s/stripe82__%s.cat.npy" % (calib_directory, sdss_filter)
stdout_write("Loading standard stars (%s)..." % (sdss_cat_filename))
sdss_cat = numpy.load(sdss_cat_filename)
stdout_write(" %d stars loaded!\n" % (sdss_cat.shape[0]))
# Now take take of pointings close to RA of 0 hours
if (ra < 5):
stdout_write("Close to RA=0h, wrapping coordinates ...\n")
sdss_cat[:,0][sdss_cat[:,0] > 180] -= 360.
d_dec = 0.6
d_ra = d_dec / math.cos(math.radians(dec))
# Break down the full catalog to stars nearby
nearby = (sdss_cat[:,0] > ra-d_ra) & (sdss_cat[:,0] < ra+d_ra) & (sdss_cat[:,1] > dec-d_dec) & (sdss_cat[:,1] < dec+d_dec)
std_stars = sdss_cat[nearby]
stdout_write("Found %d nearby (RA=%.1f, DEC=%.1f, +/- %.1fdeg) standard stars!\n" % (std_stars.shape[0], ra, dec, d_dec))
"""
Output format is:
Ra, Dec, mag_median, mag_mean, mag_std, mag_var
"""
return std_stars
def ps2sdss (catalogdata, tonry=False):
'''
PS1 to SDSS conversion following Tonry et al 2012
https://www.google.com/url?q=https%3A%2F%2Farxiv.org%2Fpdf%2F1203.0297.pdf&sa=D&sntz=1&usg=AFQjCNFZ3T88aqAPI32uYeygq7tNybjUHQ
The field names of the PS1 magnitude colums are hardwired in here.
The input catalog is overwritten, so make sure you do not need the original values any more!
Correction only for griz filters!
Parameters
----------
catalogdata -- narray Catalog containing the reference magnitudes in panstarrs AB system
tonry -- boolean If true, use Tonry 2012 conversion from PS1 to SDSS. If false (default), use Finkenbeiner 2016 covnersion.
'''
# TODO: check for existance of sitesetup.catalog_mags['panstarrs'] entries
# TODO: verify the array indices do exist, i.e., are not null
logger = logging.getLogger("PS>SDSS")
# print catalogdata
ps1colorterms = {}
if tonry:
# Tonry 2012
logger.info("converting PS1 to sdss system via tonry")
ps1colorterms['g'] = [0.019, 0.145, 0.013]
ps1colorterms['r'] = [0.007, 0.004, -0.001]
ps1colorterms['i'] = [0.010, 0.011,-0.005]
ps1colorterms['z'] = [-0.012, -0.039, 0.013]
gidx = sitesetup.catalog_mags['panstarrs'].index('g')
ridx = sitesetup.catalog_mags['panstarrs'].index('r')
_g = catalogdata[:,gidx]
_r = catalogdata[:,ridx]
psgr = _g - _r
else:
# Finkbeiner 2016
# http://iopscience.iop.org/article/10.3847/0004-637X/822/2/66/meta#apj522061s2-4 Table 2
logger.info("converting PS1 to sdss system via finkbeiner")
ps1colorterms['g'] = [-0.01808,-0.13595, 0.01941,-0.00183][::-1]
ps1colorterms['r'] = [-0.01836,-0.03577, 0.02612,-0.00558][::-1]
ps1colorterms['i'] = [ 0.01170,-0.00400, 0.00066,-0.00058][::-1]
ps1colorterms['z'] = [-0.01062, 0.07529,-0.03592, 0.00890][::-1]
gidx = sitesetup.catalog_mags['panstarrs'].index('g')
iidx = sitesetup.catalog_mags['panstarrs'].index('i')
_g = catalogdata[:,gidx]
_i = catalogdata[:,iidx]
psgr = _g - _i
for filter in ps1colorterms:
colorcorrection = numpy.polyval (ps1colorterms[filter], psgr)
magidx = sitesetup.catalog_mags['panstarrs'].index(filter)
if tonry:
catalogdata[:,magidx] += colorcorrection
else:
catalogdata[:,magidx] -= colorcorrection
return
def photcalib_old(fitsfile, output_filename, calib_directory, overwrite_cat=None):
"""
Deprecated, do not use.
"""
tmp, dummy = os.path.split(sys.argv[0])
dot_config_dir = tmp + "/config/"
print(dot_config_dir)
hdulist = pyfits.open(fitsfile)
ra, dec = hdulist[1].header['CRVAL1'], hdulist[1].header['CRVAL2']
#
# Run Sextractor on the input frame
#
sex_catalogfile = fitsfile[:-5]+".photcalib.cat"
sex_config_file = dot_config_dir+"fixwcs.conf"
sex_param_file = dot_config_dir+"fixwcs.param"
sex_logfile = fitsfile[:-5]+".sextractor.log"
sex_cmd = "sex -c %s -PARAMETERS_NAME %s -CATALOG_NAME %s %s >& %s" % (sex_config_file, sex_param_file, sex_catalogfile, fitsfile, sex_logfile)
print(sex_cmd)
stdout_write("Running SExtractor to search for stars, be patient (logfile: %s) ..." % (sex_logfile))
if ((not cmdline_arg_isset("-skip_sex")) or (not os.path.isfile(sex_catalogfile)) or cmdline_arg_isset("-forcesex")):
os.system(sex_cmd)
else:
stdout_write("Re-using previous source catalog\n")
stdout_write(" done!\n")
stdout_write("\nPreparing work ...\n\n")
# Read the Sextractor output catalog
sex_cat = numpy.loadtxt(sex_catalogfile)
stdout_write("Reading %d stars from Sextractor catalog\n" % sex_cat.shape[0])
# Eliminate all stars with flags
sex_cat = sex_cat[sex_cat[:,10] == 0]
stdout_write("%d stars left after eliminating flags\n" % sex_cat.shape[0])
# Figure out which SDSS to use for calibration
filter = hdulist[0].header['FILTER']
sdss_filter = sdss_equivalents[filter]
std_stars = numpy.array([])
if (overwrite_cat != None):
# Read this catalog instead of any of the default ones:
std_stars = numpy.loadtxt(overwrite_cat, delimiter=";")
else:
if (calib_directory != None):
std_stars = load_catalog_from_stripe82cat(ra, dec, calib_directory, sdss_filter)
print(std_stars.shape)
if (std_stars.shape[0] <= 0):
if (calib_directory != None):
stdout_write("Couldn't find any stars in the Stripe82 Standard star catalog :(\n")
stdout_write("Trying to get one directly from SDSS, please wait!\n\n")
std_stars = load_catalog_from_sdss(ra, dec, sdss_filter)
if (std_stars.shape[0] <= 0):
stdout_write("No stars not found - looks like this region isn't covered by SDSS - sorry!\n\n")
sys.exit(0)
return
dump = open("dump.cat", "w")
for i in range(std_stars.shape[0]):
print(std_stars[i,0], std_stars[i,1], file=dump)
print("\n\n\n\n\n\n\n\n", file=dump)
for i in range(sex_cat.shape[0]):
print(sex_cat[i,6], sex_cat[i,7], file=dump)
dump.close()
#
# Now go through each of the extension
# Improve: Change execution to parallel !!!
#
stdout_write("\nStarting work, results in %s ...\n\n" % output_filename)
results = open(output_filename, "w")
# Now go through the reference catalog and search for matches.
# To speed things up and avoid looking for matches between distant stars, break down the area into
# a number of smaller blocks (scanning the full area in 5x5 blocks)
ref_ra_min, ref_ra_max = numpy.min(std_stars[:,0]), numpy.max(std_stars[:,0])
ref_dec_min, ref_dec_max = numpy.min(std_stars[:,1]), numpy.max(std_stars[:,1])
ra_ranges = numpy.arange(ref_ra_min, ref_ra_max, 5)
dec_ranges = numpy.arange(ref_dec_min, ref_dec_max, 5)
dummy, ra_ranges = numpy.histogram([], bins=5, range=(ref_ra_min, ref_ra_max))
dummy, dec_ranges = numpy.histogram([], bins=5, range=(ref_dec_min, ref_dec_max))
# Reorganize the SExtractor oiutput file to have the RA/DEC values be in the first two columns
# New order is then Ra, Dec, mag, mag_err, x, y, ota
sex_reorg = numpy.zeros(shape=(sex_cat.shape[0], 7))
sex_reorg[:,0:2] = sex_cat[:,6:8]
sex_reorg[:,2:6] = sex_cat[:,2:6]
sex_reorg[:,6] = sex_cat[:,1]
# Select one of the ranges and hand the parameters off to the matching routine
matched_cat = podi_matchcatalogs.match_catalogs(std_stars, sex_reorg)
if (matched_cat != None):
numpy.savetxt(results, matched_cat, delimiter=" ")
results.close()
def load_sdss_catalog_from_fits(sdss_ref_dir, ra_range, dec_range, verbose=True):
"""
Retrieve the SDSS catalog for the specified area, using the local FITS-based
catalog of the SDSS.
Parameters
----------
sdss_ref_dir : string
Directory containing the local copy of the SDSS star catalog. This
directory has to contain the catalog index in a file named
"SkyTable.fits"
ra_range : float[2] (e.g. [165.5, 165.9])
Minimum and maximum RA range, in degrees
dec_range : float[2]
as ra_range, just for declination
verbose : Bool
add some debugging output useful for error tracking
Returns
-------
source-catalog
contains the following columns: ra, dec, mag_u, magerr_u, g, r, i, z
"""
#print "# Loading catalog from SDSS offline fits catalog..."
logger = logging.getLogger("ReadSDSSCatalogFromFits")
# Load the SkyTable so we know in what files to look for the catalog"
skytable_filename = "%s/SkyTable.fits" % (sdss_ref_dir)
skytable_hdu = pyfits.open(skytable_filename)
skytable = skytable_hdu['SKY_REGION'].data
# Select entries that match our list
logger.debug("# Searching for stars in sky with ra=%s, dec=%s" % (ra_range, dec_range))
if (verbose): print("# Searching for stars in sky with ra=%s, dec=%s" % (ra_range, dec_range))
min_dec = dec_range[0]
max_dec = dec_range[1]
min_ra = ra_range[0]
max_ra = ra_range[1]
if (verbose): print(min_ra, max_ra, min_dec, max_dec)
if (max_ra > 360.):
# This wraps around the high end, shift all ra values by -180
# Now all search RAs are ok and around the 180, next also move the catalog values
selected = skytable['R_MIN'] < 180
skytable['R_MAX'][selected] += 360
skytable['R_MIN'][selected] += 360
if (min_ra < 0):
# Wrap around at the low end
selected = skytable['R_MAX'] > 180
skytable['R_MAX'][selected] -= 360
skytable['R_MIN'][selected] -= 360
if (verbose): print("# Search radius: RA=%.1f ... %.1f DEC=%.1f ... %.1f" % (min_ra, max_ra, min_dec, max_dec))
needed_catalogs = (skytable['R_MAX'] > min_ra) & (skytable['R_MIN'] < max_ra) & \
(skytable['D_MAX'] > min_dec) & (skytable['D_MIN'] < max_dec)
if (verbose): print(skytable[needed_catalogs])
files_to_read = skytable['NAME'][needed_catalogs]
logger.debug(files_to_read)
# Now we are with the skytable catalog, so close it
skytable_hdu.close()
del skytable
#
# Load all frames, one by one, and select all stars in the valid range.
# Then add them to the catalog with RAs and DECs
#
catalog_columns = ['RA', 'DEC',
'MAG_U', 'MAGERR_U',
'MAG_G', 'MAGERR_G',
'MAG_R', 'MAGERR_R',
'MAG_I', 'MAGERR_I',
'MAG_Z', 'MAGERR_Z',
]
full_catalog = numpy.zeros(shape=(0,len(catalog_columns)))
for catalogname in files_to_read:
if (verbose): print("Reading 2mass catalog")
catalogfile = "%s/%s.fits" % (sdss_ref_dir, catalogname)
hdu_cat = pyfits.open(catalogfile)
if (hdu_cat[1].header['NAXIS2'] <= 0):
hdu_cat.close()
continue
# Read the RA and DEC values
cat_ra = hdu_cat[1].data.field('RA')
cat_dec = hdu_cat[1].data.field('DEC')
# To select the right region, shift a temporary catalog
cat_ra_shifted = cat_ra
if (max_ra > 360.):
cat_ra_shifted[cat_ra < 180] += 360
elif (min_ra < 0):
cat_ra_shifted[cat_ra > 180] -= 360
in_search_range = (cat_ra_shifted > min_ra) & (cat_ra_shifted < max_ra ) & (cat_dec > min_dec) & (cat_dec < max_dec)
array_to_add = numpy.zeros(shape=(numpy.sum(in_search_range),len(catalog_columns)))
if (verbose): print(catalogfile, numpy.sum(in_search_range))
for col in range(len(catalog_columns)):
array_to_add[:,col] = hdu_cat[1].data.field(catalog_columns[col])[in_search_range]
full_catalog = numpy.append(full_catalog, array_to_add, axis=0)
if (verbose): print("# Read a total of %d stars from %d catalogs!" % (full_catalog.shape[0], len(files_to_read)))
logger.debug("# Read a total of %d stars from %d catalogs!" % (full_catalog.shape[0], len(files_to_read)))
return full_catalog
def query_sdss_catalog(ra_range, dec_range, sdss_filter, verbose=False):
"""
High-level interface to the SDSS catalog search. Depending on the settings
in podi_sitesetup (see sdss_ref_type) this forwards the query to search
either the local Stripe82 catalog, the online SDSS catalog or the local FITS
catalog.
Parameters
----------
sdss_ref_dir : string
Directory containing the local copy of the SDSS star catalog. This
directory has to contain the catalog index in a file named
"SkyTable.fits"
ra_range : float[2] (e.g. [165.5, 165.9])
Minimum and maximum RA range, in degrees
dec_range : float[2]
as ra_range, just for declination
verbose : Bool
add some debugging output useful for error tracking
Returns
-------
source-catalog
"""
#
# Get a photometric reference catalog one way or another
#
std_stars = None
logger = logging.getLogger("QuerySDSS")
if (sitesetup.sdss_ref_type == "stripe82"):
std_stars = numpy.zeros(shape=(0,0)) #load_catalog_from_stripe82cat(ra, dec, calib_directory, sdss_filter)
print(std_stars.shape)
elif (sitesetup.sdss_ref_type == 'local'):
std_stars = load_sdss_catalog_from_fits(sitesetup.sdss_ref_dir, ra_range, dec_range,
verbose=verbose)
elif (sitesetup.sdss_ref_type == 'web'):
#stdout_write("Trying to get one directly from SDSS, please wait!\n\n")
#std_stars = load_catalog_from_sdss(ra, dec, sdss_filter)
logger.info("Loading SDSS catalog from online webserver")
std_stars = load_catalog_from_sdss(ra_range, dec_range, sdss_filter,
verbose=verbose)
#print "returning",std_stars.shape[0],"stars..."
return std_stars
def estimate_zeropoint(filtername):
logger = logging.getLogger("EstimateZP")
# First, check if we know this filter
if (not filtername in filter_bandpass):
logger.warning("this (%s) is a filter I don't know!" % (filtername))
# Now get the bandpass info
bp = filter_bandpass[filtername]
_, mean, center, fwhm, _, _, fmax, fmean, area, _, _ = bp
logger.debug("%s ==> mean=%f, center=%f, fwhm=%f" % (filtername, mean, center, fwhm))
# Now find the closest two ODI filters with known zeropoints, so we can
# interpolate a zeropoint between then
known_name = [None] * len(reference_zeropoint)
# numpy.empty((len(reference_zeropoint)), dtype=numpy.str)
known_pos = numpy.empty((len(reference_zeropoint)))
known_fwhm = numpy.empty((len(reference_zeropoint)))
known_zp_per_fwhm = numpy.empty((len(reference_zeropoint)))
for idx, known_filter in enumerate(reference_zeropoint):
known_name[idx] = known_filter
_, k_mean, k_center, k_fwhm, _, _, _, _, _, _, _ = filter_bandpass[known_filter]
known_pos[idx] = k_mean
known_fwhm[idx] = k_fwhm
zp = reference_zeropoint[known_filter][0]
zp_fwhm = zp - 2.5*math.log10(k_fwhm)
#print known_filter, zp_fwhm
known_zp_per_fwhm[idx] = zp_fwhm
known_name = numpy.array(known_name)
# print "pos:", known_pos
# print "names:", known_name
# print "\n\n\n\n\n"
# compute the difference in mean wavelengths, then pick the two closest ones
delta_pos = numpy.fabs(known_pos - mean)
si = numpy.argsort(delta_pos)
known_name = known_name[si]
known_pos = known_pos[si]
known_zp_per_fwhm = known_zp_per_fwhm[si]
#print known_pos
logger.debug("Interpolating ZP between %s and %s" % (known_name[0], known_name[1]))
d_pos = known_pos[1] - known_pos[0]
d_zpfwhm = known_zp_per_fwhm[1] - known_zp_per_fwhm[0]
slope_zpfwhm = d_zpfwhm / d_pos
#print "positions:", known_name[0], known_name[1]
#print "delta-pos:", d_pos, " delta zp/fwhm=",d_zpfwhm
zp_interpol = known_zp_per_fwhm[0] + slope_zpfwhm*(mean-known_pos[0])
#print "INTERPOL:", known_zp_per_fwhm[0], slope_zpfwhm, mean, known_pos[0], mean-known_pos[0]
zp_fwhm = zp_interpol + 2.5*math.log10(fwhm)
logger.debug("ZP results: ZP/FWHM=%f, final-ZP=%f" % (zp_interpol, zp_fwhm))
# print known_name
# print known_pos
# print known_fwhm
# print zp_interpol, fwhm, zp_fwhm
return zp_fwhm
def estimate_mean_star_color(filtername, T=5000):
logger = logging.getLogger("MeanStarColor")
if (not filtername in filter_bandpass):
logger.warning("This filter (%s) is not registered" % (filtername))
return None
# compute the color of a average star (Teff=5000K) in XXX-J, where XXX is
# the specified filter
# define constants
h = 6.6626e-34 # J s
kb =1.380e-23 # J/K
c = 3e8 # m/s
def planck(l_ang, T):
lm = l_ang * 1e-10
return 2*h*c**2/(lm**5) * (1. / (numpy.exp(h*c/(lm*kb*T))-1))
spec_wl = numpy.arange(3000,15000,25)
flux = planck(spec_wl, T)
AB_calibspec = 0.11 / (spec_wl*1e-10)**2
# compute the mean spectral flux for J band
Jband_wl = (spec_wl > 11500) & (spec_wl < 13500)
Jband_mean = numpy.mean(flux[Jband_wl])
# Now get the left and right filter edge for the specified filter
_, _, _, _, left, right, _, _, _, _, _ = filter_bandpass[filtername]
Xband_wl = (spec_wl >= left) & (spec_wl <= right)
Xband_mean = numpy.mean(flux[Xband_wl])
Jmag = -2.5*numpy.log10(Jband_mean)
Xmag = -2.5*numpy.log10(Xband_mean)
logger.debug("Found X and J-magnitudes: %.3f / %.3f" % (Xmag, Jmag))
#
# Apply some correction to account for both bands being in AB magnitudes
#
Jband_AB = numpy.mean(AB_calibspec[Jband_wl])
Xband_AB = numpy.mean(AB_calibspec[Xband_wl])
logger.debug("AB zeropoint corrections: X = %.3f / J = %.3f" % (Xband_AB, Jband_AB))
#print Jmag, Xmag, -2.5*math.log10(Jband_AB), -2.5*math.log10(Xband_AB)
#print (-2.5*math.log10(Jband_AB)) - (-2.5*math.log10(Xband_AB))
AB_color_corr = (-2.5*math.log10(Jband_AB)) - (-2.5*math.log10(Xband_AB))
starcolor = (Xmag-Jmag)+AB_color_corr
logger.debug("Estimating mean star color X-J = %.3f (Teff=% 6dK)" % (starcolor, T))
return starcolor
def find_closest_sdss_counterpart(filtername):
logger = logging.getLogger("Match2SDSS")
if (not filtername in filter_bandpass):
logger.warning("This is an unknown filter, defaulting to sdss_r")
return 'r'
# Look up the filter bandpass parameters
# most importantly, we will use the mean_pos value to find a good
# filter match
_, meanpos, center, _, _, _, _, _, _, _, _ = filter_bandpass[filtername]
sdss_filters = numpy.array(['u', 'g', 'r', 'i', 'z'])
sdss_mean_pos = numpy.array([3557, 4825, 6261, 7672, 9097])
# compute the wavelength difference to all SDSS filters
d_lambda = numpy.fabs(sdss_mean_pos-meanpos)
# pick the closest
closest = numpy.argmin(d_lambda)
return sdss_filters[closest]
def photcalib(source_cat, output_filename, filtername, exptime=1,
diagplots=True, calib_directory=None, overwrite_cat=None,
plottitle=None, otalist=None,
options=None,
verbose=False,
eliminate_flags=True,
matching_radius=3,
error_cutoff=None,
detailed_return={},
saturation_limit=sitesetup.photcalib_saturation,
plot_suffix=None):
"""
Perform the photometric calibration, create the diagnostic plots and return
the results.
Parameters
----------
source_cat : ndarray
Input source array containing the catalog output from the SExtractor
run.
output_filename : string
Filename of the file containing the resulting image file. In the case of
collectcells, this is the output filename, hence the name. However, if
photcalib is run on existing images, this is the image INPUT filename.
filtername : string
Name of the ODI filter. Based on this filter name photcalib will decide
the corresponding SDSS filter.
exptime : float
Exposure time of the frame, required to normalize the photometric zero-
point.
diagplots : Bool
If True, create the diagnostic plots for the photometric calibration.
calib_directory : string
overwrite_cat : Bool
plottitle : string
Title of the diagnostic plot, e.g. containing the filename, or some
information about the Object, etc.
otalist : HDUList
If given, photcalib uses the WCS information from the individual
extensions to add OTA outlines to the global, focal-plane wide diagnostic
plot.
options : dictionary
General options, containing information about some command line flags.
eliminate_flags : Bool
If set to True, only sources with no SourceExtractor flags (that might
indicate problems with the photometry) are used for the calibration. If
set, fewer stars are used but these stars are of the best possible
quality.
matching_radius : float
Matching radius to be used when matching the ODI source catalog to the
SDSS reference catalog.
Returns
-------
ZP-median : float
median zeropoint across the entire focal plane
ZP-Std : float
standard deviation of the photometric zeropoint
ODI-SDSS-matched : ndarray
matched catalog of ODI sources and SDSS reference stars.
ZP-normalized : float
photometric zeropoint, corrected and normalized to an exposure time
of 1 second.
"""
logger = logging.getLogger("PhotCalib")
error_return_value = (99.9, 99.9, None, 99.9)
detailed_return['median'] = -99.
detailed_return['std'] = -99.
detailed_return['zp_exptime'] = -99.
detailed_return['stderrofmean'] = -99.
detailed_return['zp_upper1sigma'] = -99.
detailed_return['zp_lower1sigma'] = -99.
detailed_return['n_clipped'] = -1
detailed_return['n_raw'] = -1
detailed_return['colorterm'] = None
detailed_return['colorcorrection'] = None
detailed_return['radialZPfit'] = None
detailed_return['radialZPfit_error'] = None
detailed_return['zp_restricted'] = None
detailed_return['zp_magnitude_slope'] = None
detailed_return['aperture_mode'] = None
detailed_return['aperture_size'] = -99.
detailed_return['aperture_mag'] = None
detailed_return['aperture_magerr'] = None
detailed_return['catalog'] = "none"
detailed_return['reference_filter'] = "none"
detailed_return['reference_catalog_files'] = None
detailed_return['use_for_calibration_mask'] = None
# Eliminate all stars with flags
if (eliminate_flags):
flags = source_cat[:,SXcolumn['flags']]
no_flags_set = (flags == 0)
source_cat = source_cat[no_flags_set]
# numpy.savetxt("photcat", source_cat)
ra_min = numpy.min(source_cat[:,SXcolumn['ra']])
ra_max = numpy.max(source_cat[:,SXcolumn['ra']])
logger.debug("ra-range: %f --- %f" % (ra_min, ra_max))
# Make sure we deal with RAs around 0 properly
if (math.fabs(ra_max - ra_min) > 180):
ra_fixed = source_cat[:,SXcolumn['ra']].copy()
ra_fixed[ra_fixed > 180] -= 360
ra_min = numpy.min(ra_fixed)
ra_max = numpy.max(ra_fixed)
logger.debug("ra-range (fixed): %f --- %f" % (ra_min, ra_max))
dec_min = numpy.min(source_cat[:,SXcolumn['dec']])
dec_max = numpy.max(source_cat[:,SXcolumn['dec']])
logger.debug("Starting photometric calibration!")
#stdout_write("\nPreparing work ...\n\n")
ra_range = [ra_min, ra_max]
dec_range = [dec_min, dec_max]
# sdss_filter = sdss_equivalents[filtername]
sdss_filter = find_closest_sdss_counterpart(filtername)
logger.debug("Translating filter: %s --> %s" % (filtername, sdss_filter))
if (sdss_filter is None):
# This filter is not covered by SDSS, can't perform photometric calibration
return error_return_value
std_stars = None
ref_filenames = None
for catname in sitesetup.photcalib_order:
logger.debug("Trying to find stars in %s catalog" % (catname))
if (catname == "sdss"):
# Figure out which SDSS to use for calibration
logger.info("Querying the SDSS catalog")
_std_stars, ref_filenames = podi_search_ipprefcat.get_reference_catalog(ra_range, dec_range,
radius=None,
basedir=sitesetup.sdss_ref_dir,
cattype="sdss",
return_filenames=True)
if (_std_stars is not None and _std_stars.shape[0] > 0):
detailed_return['catalog'] = "SDSS"
# Found some SDSS stars
std_stars = _std_stars
logger.debug("Found %d stars in the SDSS" % (std_stars.shape[0]))
break
else:
logger.debug("No stars not found - looks like this region isn't covered by SDSS - sorry!")
elif (catname == "ucac4"):
detailed_return['catalog'] = "UCAC4"
# Try to get some data from UCAC4 instead
if (sitesetup.ucac4_ref_dir is not None and
not sitesetup.ucac4_ref_dir == "none" and
os.path.isdir(sitesetup.ucac4_ref_dir)):
logger.debug("Running UCAC4 query (%s)" % (sitesetup.ucac4_ref_dir))
_std_stars, ref_filenames = podi_search_ipprefcat.get_reference_catalog(
ra=ra_range,
dec=dec_range,
radius=None,
basedir=sitesetup.ucac4_ref_dir,
cattype='ucac4',
verbose=False,
return_filenames=True
)
# Sort out all stars with invalid magnitudes
valid = (_std_stars[:,2] < 20) & (numpy.fabs(_std_stars[:,3]) <= 0.9)
logger.debug("Number of UCAC4 stars: %s" % (str(_std_stars.shape)))
# numpy.savetxt("ucac.dump", _std_stars)
std_stars = _std_stars[valid]
# numpy.savetxt("ucac.dump2", std_stars)
logger.debug("Found a valid UCAC4 source catalog (%s)" % (str(std_stars.shape)))
break
else:
logger.warning("Problem with querying the UCAC4 catalog")
elif (catname == "ippref"):
logger.debug("Running IPPRef query (%s)" % (sitesetup.ippref_ref_dir))
_std_stars, ref_filenames = podi_search_ipprefcat.get_reference_catalog(ra_range, dec_range,
radius=None,
basedir=sitesetup.ippref_ref_dir,
cattype="IPPRef",
return_filenames=True)
#print "IPP:", _std_stars
if (_std_stars is not None and _std_stars.shape[0] > 0):
detailed_return['catalog'] = "IPPRef"
std_stars = _std_stars
logger.debug("Found %d reference sources in IPPRef catalog" % (std_stars.shape[0]))
break
else:
logger.warning("Problem with querying the IPPRef catalog")
elif (catname in sitesetup.catalog_directory):
# So basically we know where to find this catalog, but we are not sure yet
# the catalog covers the filter we need
if (catname not in sitesetup.catalog_mags):
# this catalog is not properly registered for the photometric calibration
logger.error("Chosen catalog for photometric calibration (%s) is known but not well described!" % (catname))
return error_return_value
if (sdss_filter not in sitesetup.catalog_mags[catname]):
# this catalog does not have the right photometric data
logger.error("Chosen catalog for photometric calibration (%s) does not contain photometry in this filter band (%s)" % (
catname, sdss_filter))
# return error_return_value
continue
logger.info("Using %s for photometric calibration" % (catname))
catalog_basedir, cat_mag = sitesetup.catalog_directory[catname]
_std_stars, ref_filenames = podi_search_ipprefcat.get_reference_catalog(
ra_range, dec_range, radius=None,
basedir=catalog_basedir,
cattype="general",
return_filenames=True
)
if (_std_stars is not None and _std_stars.shape[0] > 0):
detailed_return['catalog'] = catname
std_stars = _std_stars
logger.info("Found %d reference sources in %s catalog" % (std_stars.shape[0], catname))
break
else:
logger.warning("No suitable reference sources found in the %s catalog" % (catname))
else:
logger.error("Illegal catalog name in the photcalib order: %s" % (catname))
return error_return_value
if (std_stars is None or
(type(std_stars) == numpy.ndarray and std_stars.shape[0] <= 0)):
# No sources found, report a photometric calibration failure.
logger.warning("No photometric reference sources found :(")
return error_return_value
if (error_cutoff is None):
if (detailed_return['catalog'] in sitesetup.photcalib_error_cutoff):
error_cutoff = sitesetup.photcalib_error_cutoff[detailed_return['catalog']]
else:
error_cutoff = 0.05
detailed_return['reference_catalog_files'] = ref_filenames
#
# Now go through each of the extension
# Improve: Change execution to parallel !!!
#
logger.debug("Starting work, results in %s ..." % output_filename)
# results = open(output_filename+".photcal", "w")
odi_sdss_matched = podi_matchcatalogs.match_catalogs(source_cat, std_stars, matching_radius=matching_radius)
# numpy.savetxt("odisource.dump", source_cat)
# numpy.savetxt("matched.dump", odi_sdss_matched)
# Stars without match in SDSS have RA=-9999, let's sort them out
found_sdss_match = odi_sdss_matched[:,2] >= 0
odi_sdss_matched = odi_sdss_matched[found_sdss_match]
n_matches = numpy.sum(found_sdss_match)
# numpy.savetxt("matched.dump2", odi_sdss_matched)
logger.debug("Found %d matches between ODI source catalog and photometric reference catalog" % (