diff --git a/csst/msc/calib_pos.py b/csst/msc/calib_pos.py index ff969c17b67b14d63454ead3a47063ef9e954bd7..ca68aabf42a6ec5c496970c186fe359227aa5763 100644 --- a/csst/msc/calib_pos.py +++ b/csst/msc/calib_pos.py @@ -3,7 +3,6 @@ import time from functools import partial from multiprocessing import Pool from subprocess import Popen - import healpy as hp import numpy as np from astropy import table @@ -14,13 +13,14 @@ from astropy.wcs import WCS from .. import PACKAGE_PATH from ..core.processor import CsstProcessor +from csst.msc.data_manager import CsstMscDataManager CONFIG_PATH = PACKAGE_PATH + "/msc/pos_calib_config/" class CsstProcMscPositionCalibration(CsstProcessor): - def join_data(self, img_list, wht_list, flg_list, path_output): + def join_data(self, img_list, wht_list, flg_list): """ Prepare data for running scamp; Combine all image data, weight files, flag files to their one frame. @@ -32,18 +32,14 @@ class CsstProcMscPositionCalibration(CsstProcessor): weith files to join together, e.g.,MSC_210304093000_0000000_06_img.fits flg_list: flag files to join together, e.g.,e.g.,MSC_210304093000_0000000_06_flg.fits - path_output: - the output dir for the joined file. - Returns ------- The joined multi-extension file(not stacked), weight, flag files. - e.g., MSC_210304093000_0000000_img.fits,MSC_210304093000_0000000_wht.fits, MSC_210304093000_0000000_flg.fits. + e.g., combined_img.fits,combined_wht.fits, combined_flg.fits. """ - img_prefix = img_list[0][0].header['FILENAME'][0:-7] - output_imgnm = path_output + img_prefix + '_img.fits' + output_imgnm = self.dm.pc_combined_file("img","fits") hdul_img = fits.HDUList() for i in range(0, len(img_list)): h0 = fits.PrimaryHDU(header=img_list[i][0].header) @@ -52,7 +48,7 @@ class CsstProcMscPositionCalibration(CsstProcessor): hdul_img.append(h1) hdul_img.writeto(output_imgnm, overwrite=True) - output_whtnm = path_output + img_prefix + '_wht.fits' + output_whtnm = self.dm.pc_combined_file("wht","fits") hdul_wht = fits.HDUList() for i in range(0, len(wht_list)): h0 = fits.PrimaryHDU(header=wht_list[i][0].header) @@ -61,7 +57,7 @@ class CsstProcMscPositionCalibration(CsstProcessor): hdul_wht.append(h1) hdul_wht.writeto(output_whtnm, overwrite=True) - output_flgnm = path_output + img_prefix + '_flg.fits' + output_flgnm = self.dm.pc_combined_file("flg","fits") hdul_flg = fits.HDUList() for i in range(0, len(flg_list)): h0 = fits.PrimaryHDU(header=flg_list[i][0].header) @@ -70,86 +66,84 @@ class CsstProcMscPositionCalibration(CsstProcessor): hdul_flg.append(h1) hdul_flg.writeto(output_flgnm, overwrite=True) - def run_sextractor(self, fn_list, path_output): + def run_sextractor(self,ccd_ids): """ Run sextractor Parameters ---------- - fn_list: - file name list, e.g.,MSC_210304093000_0000000_06_img.fits... - config_sextractor: - the path for sextractor configuration file. - path_output: - the current working dir - + ccd_ids: + ccd list... Returns ------- - The photometric catalog, with position and flux, e.g.,MSC_210304093000_0000000_06_img.acat + The photometric catalog, with position and flux, e.g.,MSC_210304093000_0000000_06_img.cat """ - fn = fn_list + cat_list = [] + fn_list = [] + for this_ccd_id in ccd_ids: + sex_cat = self.dm.l1_sci(ccd_id = this_ccd_id, suffix="img", ext="cat") + fn = self.dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="fits") + cat_list.append(sex_cat) + fn_list.append(fn) + config_sextractor = CONFIG_PATH + "new_csst_realtime.no.weight.sex" sex_comd1 = 'sex -c ' + config_sextractor + ' ' - sex_comd2 = fn + ' -CATALOG_NAME ' + fn[0:-5] + '.acat' + sex_comd2 = fn_list + ' -CATALOG_NAME ' + cat_list sex_comd3 = ' -PARAMETERS_NAME ' + CONFIG_PATH + 'csst_realtime.param' + ' -FILTER_NAME ' + CONFIG_PATH + 'csst_realtime.conv' + ' -STARNNW_NAME ' + CONFIG_PATH + 'csst_realtime.nnw' sex_comd = sex_comd1 + sex_comd2 + sex_comd3 print(sex_comd) p = Popen(sex_comd, shell=True) p.wait() - def combine_catalog(self, img_list, path_output): + def combine_catalog(self, ccd_ids): """ Combine the sextractor catalog together Parameters ----------- - img_list: - image list, in table format - path_output: - the output dir - + ccd_ids: ccd list Returns ------- - The combined catalog,e.g., MSC_210304093000_0000000.acat.fits + The combined catalog,e.g., combined_cat.fits """ - fn = path_output + img_list[0][0].header['FILENAME'][0:-7] - output_catnm = str(fn + '.acat.fits') + fn = self.dm.pc_combined_file("cat","fits") hdul = fits.HDUList() - if len(img_list) == 18: - for i in range(0, len(img_list)): - image_prefix = img_list[i][0].header['FILENAME'] - cat_nm = path_output + image_prefix + '.acat' - cat_i = fits.open(cat_nm) - hdul.append(cat_i[0]) - hdul.append(cat_i[1]) - hdul.append(cat_i[2]) - hdul.writeto(output_catnm, overwrite=True) - else: - print('the length of file list in not equal to 18, needs to check') - - def run_scamp(self, img_list, path_output): + for this_ccd_id in ccd_ids: + cat_nm = self.dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="cat") + cat_i = fits.open(cat_nm) + data=cat_i[2].data + msk=(data['MAGERR_AUTO']<0.2)&(data['CLASS_STAR']>0.5)&(data['CLASS_STAR']<=1.0) + data=data[msk] + cat_i[2].data=data + hdul.append(cat_i[0]) + hdul.append(cat_i[1]) + hdul.append(cat_i[2]) + hdul.writeto(fn, overwrite=True) + + def run_scamp(self, wcs_refine=False): """ Run scamp - - Parameters - --------- - img_list: - to join a file 'image_prefix+.acat.fits', e.g.,MSC_210304093000_0000000.acat.fits - config_scamp: - the config file path for scamp - Returns ------- - Image header updated with WCS keywords, MSC_210304093000_0000000.acat.head. + Image header updated with WCS keywords, combined_cat.head. """ - image_prefix = (img_list[0][0].header)['FILENAME'][0:-7] - config_scamp = CONFIG_PATH + "default2.scamp" - scamp_comd = 'scamp ' + image_prefix + '.acat.fits -ASTREFCAT_NAME= ' + 'ref.cat\ + cb_fn = self.dm.pc_combined_file("cat","fits") + config_scamp = CONFIG_PATH + "csst.scamp" + scamp_comd = 'scamp ' + cb_fn + ' -ASTREFCAT_NAME= ' + 'ref.cat\ -MERGEDOUTCAT_NAME ' + 'merged.cat -FULLOUTCAT_NAME ' + 'full.cat\ -c ' + config_scamp print(scamp_comd) p = Popen(scamp_comd, shell=True) p.wait() + if wcs_refine: + Popen('cp '+ self.dm.pc_combined_file("cat","head")+" "+self.dm.pc_combined_file("cat","ahead"),shell=True) + scamp_comd = 'scamp ' + cb_fn + ' -ASTREFCAT_NAME= ' + 'ref.cat\ + -MERGEDOUTCAT_NAME ' + 'merged.cat -FULLOUTCAT_NAME ' + 'full.cat\ + -c ' + config_scamp + '-AHEADER_SUFFIX '+ self.dm.pc_combined_file("cat","ahead") + print(scamp_comd) + p=Popen(scamp_comd,shell=True) + p.wait() + def convert_hdu_to_ldac(self, hdu): """ @@ -182,14 +176,14 @@ class CsstProcMscPositionCalibration(CsstProcessor): tbl2.header['EXTNAME'] = 'LDAC_OBJECTS' return tbl1, tbl2 - def get_refcat(self, img_list, path_gaia, search_radius, silent=True): + def get_refcat(self, img_list, path_gaia, search_radius, silent=True,pm_correct=True): """ Get reference catalog for scamp. The reference cat is GAIA EDR3. Parameters ---------- image_prefix: - a image to get its reference catalog, e.g.,MSC_210304093000_0000000_img.fits. + a image to get its reference catalog, e.g.,combined_img.fits. Usually the center of the image is the wcs parameters CRVAL1,CRVAL1. search_radius: circle radius for searching, units: degree. e.g., 2 degree for a 1x1 deg^2 image. @@ -199,20 +193,25 @@ class CsstProcMscPositionCalibration(CsstProcessor): Returns ------- outcat: - filename of the cross matched catalog. - This catalog is used as a reference catalog for running scamp. - e.g.,MSC_210304093000_0000000.gaialac.fits - + filename of the cross matched catalog: ref.cat """ - image_prefix = (img_list[0][0].header)['FILENAME'][0:-7] - fname = image_prefix + '_img.fits' - gaianame = image_prefix + '.gaia.fits' - gaialacnm = image_prefix + '.gaialac.fits' - outcat = gaianame + outcat = self.dm.pc_ref_cat + fnmae = cb_fn = self.dm.pc_combined_file("img","fits") hdu = fits.open(fname) header1 = hdu[0].header header2 = hdu[1].header - deltatime = 0.00 + if (pm_correct): + oday=header1['DATE-OBS'] + otime=header1['TIME-OBS'] + exptime=header1['EXPTIME'] + odaytime=header1['DATE-OBS']+'T'+header1['TIME-OBS'] + t = Time(oday, format='isot', scale='utc') + t.format = 'decimalyear' + t1=Time(2016.0, format='decimalyear', scale='utc') + deltatime= t.value -2016.0 + print('(Time - 2016.0) = deltatime =', deltatime) + else: + deltatime=0.00 ra = float(header2['CRVAL1']) dec = float(header2['CRVAL2']) c = SkyCoord(ra, dec, unit=(u.deg, u.deg)) @@ -244,7 +243,6 @@ class CsstProcMscPositionCalibration(CsstProcessor): refcat = table.vstack(refcat, join_type='inner') refcat.rename_column('ra', 'X_WORLD') refcat.rename_column('dec', 'Y_WORLD') - print('delta_time between obs_cat and ref_cat:', deltatime) mask = (refcat['pmdec'] != refcat['pmdec']) refcat['pmdec'][mask] = 0 @@ -258,21 +256,16 @@ class CsstProcMscPositionCalibration(CsstProcessor): refcat.rename_column('ra_error', 'ERRA_WORLD') refcat.rename_column('dec_error', 'ERRB_WORLD') refcat.rename_column('phot_g_mean_mag', 'MAG') - if outcat: refcat.write(outcat, format='fits', overwrite=True) - - if os.path.isfile(gaianame): - print('exist') - hdu = fits.open(gaianame) - hdu1 = self.convert_hdu_to_ldac(hdu) - hdup = fits.PrimaryHDU() - hdu = hdu1[0] - tbhdu = hdu1[1] - thdulist = fits.HDUList([hdup, hdu, tbhdu]) - if os.path.isfile(gaialacnm): os.remove(gaialacnm) - thdulist.writeto(gaialacnm) + + hdu = refcat + hdu1 = self.convert_hdu_to_ldac(hdu) + hdup = fits.PrimaryHDU() + hdu = hdu1[0] + tbhdu = hdu1[1] + thdulist = fits.HDUList([hdup, hdu, tbhdu]) + thdulist.writeto(outcat) print('##################### end #####################') - return gaialacnm def rewrite_wcs_head(self, head): """ @@ -308,7 +301,7 @@ class CsstProcMscPositionCalibration(CsstProcessor): f.close() return wcshead - def check_astrometry(self, img_list, path_output): + def check_astrometry(self, img_list): """ Check position calibration quality @@ -316,8 +309,6 @@ class CsstProcMscPositionCalibration(CsstProcessor): ------------ img_list: list of images, in table format - path_output: - work dir Returns ------- @@ -333,119 +324,142 @@ class CsstProcMscPositionCalibration(CsstProcessor): print('############## check the astrometry quality and save files ################') r1 = [] d1 = [] - image_prefix = (img_list[0][0].header)['FILENAME'][0:-7] - fn = path_output + image_prefix - wcshead = self.rewrite_wcs_head(fn + '.acat.head') - acat = fits.open(fn + '.acat.fits') - acat_change = str(fn + '.acat.change.fits') - cat_suffix = '.acat' + fn_head = self.dm.pc_combined_file("cat","head") + fn_scat = self.dm.pc_combined_file("cat","fits") + data_scat = fits.open(fn_scat) + + wcshead = self.rewrite_wcs_head(fn_head) + scat = fits.open(fn_scat) hdul = fits.HDUList() - if len(img_list) == 18: - for i in range(0, len(img_list)): - wcshdr = fits.getheader(wcshead, i, ignore_missing_simple=True) # read headers and change to RA---TPV,DEC--TPV for wcs_transfer package - wcshdr['CTYPE1'] = 'RA---TPV' - wcshdr['CTYPE2'] = 'DEC--TPV' - w = WCS(wcshdr) - # print(wcshdr) - cat_nm = path_output + (img_list[i][0].header)['FILENAME'] + cat_suffix - cat_i = fits.open(cat_nm) - sexcat = cat_i[2].data - ra_sex = sexcat['ALPHA_J2000'] - dec_sex = sexcat['DELTA_J2000'] - x = sexcat['XWIN_IMAGE'] - y = sexcat['YWIN_IMAGE'] - r, d = w.all_pix2world(x, y, 0) # convert xwin,ywin to ra,de - sexcat['ALPHA_J2000'] = r - sexcat['DELTA_J2000'] = d - cat_i[2].data = sexcat - hdul.append(cat_i[0]) - hdul.append(cat_i[1]) - hdul.append(cat_i[2]) - r1 = np.hstack((r1, r)) - d1 = np.hstack((d1, d)) - obsc = SkyCoord(ra=r1 * u.degree, dec=d1 * u.degree) - tmp_cat = np.zeros((len(obsc), 2)) - tmp_cat[:, 0] = obsc.ra - tmp_cat[:, 1] = obsc.dec - np.savetxt(path_output + 'scamp_coord.txt', tmp_cat, fmt="%.10f %.10f", delimiter="\n") - hdul.writeto(acat_change, overwrite=True) # update the cat with new ra,dec (from 1st scamp wcs.) - else: - print('the length of fitslist is not equal to 18,needs to check') - - def write_headers(self, img_list): + for i in range(0, len(img_list)): + wcshdr = fits.getheader(wcshead, i, ignore_missing_simple=True) # read headers and change to RA---TPV,DEC--TPV for wcs_transfer package + wcshdr['CTYPE1'] = 'RA---TPV' + wcshdr['CTYPE2'] = 'DEC--TPV' + w = WCS(wcshdr) + sexcat = data_scat[i*3+2].data + x = sexcat['XWIN_IMAGE'] + y = sexcat['YWIN_IMAGE'] + r, d = w.all_pix2world(x, y, 0) # convert xwin,ywin to ra,de + sexcat['ALPHA_J2000'] = r + sexcat['DELTA_J2000'] = d + r1 = np.hstack((r1, r)) + d1 = np.hstack((d1, d)) + obsc = SkyCoord(ra=r1 * u.degree, dec=d1 * u.degree) + tmp_cat = np.zeros((len(obsc), 2)) + tmp_cat[:, 0] = obsc.ra + tmp_cat[:, 1] = obsc.dec + np.savetxt(dm.pc_scamp_coord, tmp_cat, fmt="%.10f %.10f", delimiter="\n") + + gaia_cat=Table.read(self.dm.pc_ref_cat) + gaia_ra=gaia_cat['X_WORLD'] + gaia_dec=gaia_cat['Y_WORLD'] + refc=SkyCoord(ra=gaia_ra*u.degree,dec=gaia_dec*u.degree) + idx, d2d, d3d = obsc.match_to_catalog_sky(refc) + ref_uid=np.unique(idx) + obs_uid=np.full_like(ref_uid,-1) + tmpj=-1 + ccdraoff_med=ccddecoff_med=ccdra_rms=ccddec_rms=-1 + for i in ref_uid: + tmpj=tmpj+1 + iid=(idx == i) + iiid = (d2d.deg[iid] == d2d.deg[iid].min()) + obs_uid[tmpj]=iid.nonzero()[0][iiid.nonzero()[0]][0] + + uidlim=d2d[obs_uid].arcsecond <1. # set match radius=1 arcsec + if uidlim.sum()>0: + obs_uidlim=obs_uid[uidlim] + ref_uidlim=ref_uid[uidlim] + ccdraoff=(obsc[obs_uidlim].ra- refc[ref_uidlim].ra).arcsec*np.cos(obsc[obs_uidlim].dec.deg*np.pi/180.) + ccdraoff_med=np.median(ccdraoff) + ccdra_rms=np.std(ccdraoff) + ccddecoff=(obsc[obs_uidlim].dec- refc[ref_uidlim].dec).arcsec + ccddec_rms=np.std(ccddecoff) + ccddecoff_med=np.median(ccddecoff) + match_num=len(ccdraoff) + print('################# astrometry result: ##############') + if (match_num<100): + print('### bad astrometry ###') + print('median ra_off, dec_off (mas) from scamp:',ccdraoff_med*1000.,ccddecoff_med*1000.) + print('rms ra_off, dec_off (mas) from scamp:',ccdra_rms*1000.,ccddec_rms*1000.) + print('############################################') + return ccdraoff,ccddecoff,ccdraoff_med,ccddecoff_med,ccdra_rms,ccddec_rms + + def write_headers(self, ccd_ids): """ Wrtie history to header """ - head_suffix = img_list[0][0].header['FILENAME'][0:-7] + '.acat.head.fits' + head_suffix = self.dm.pc_combined_head_fits hdul2 = fits.open(head_suffix, ignore_missing_simple=True) - if len(img_list) == 18: - for i in range(0, len(img_list)): - fits_nm = img_list[i][0].header['FILENAME'] + '.head' - hdul1 = fits.open(fits_nm, mode='update', ignore_missing_simple=True) - hdr = hdul1[0].header - hdr2 = hdul2[i].header - hdr.extend(hdr2, unique=True, update=True) - WCS_S = 0 - WCS_V = '2.0.4' - WCS_P = 'default.scamp' - WCS_TOL = time.strftime('%Y-%m-%d %H:%M:%S %p') - hdr.set('WCS_S', '0', '0=done') - hdr.set('WCS_V', WCS_V, 'Version of WCS calibration') - hdr.set('WCS_P', WCS_P, 'Configure file name of WCS') - hdr.set('WCS_TOL', WCS_TOL, 'Time of last wcs calibration') - # hdul1.flush() - # hdul1.close() - else: - print('The total number of the fits files is not 18.') - - def prepare(self, path_gaia, path_output, search_radius=2.0): - self.path_gaia = path_gaia - self.path_output = path_output - self.search_radius = search_radius - - def run(self, img_list, wht_list, flg_list, fn_list, path_gaia, path_output, search_radius): + for this_ccd_id in ccd_ids: + fits_nm = self.dm.l1_sci(this_ccd_id, suffix="img", ext="head") + hdul1 = fits.open(fits_nm, mode='update', ignore_missing_simple=True) + hdr = hdul1[0].header + hdr2 = hdul2[i].header + hdr.extend(hdr2, unique=True, update=True) + WCS_S = 0 + WCS_V = '2.0.4' + WCS_P = 'default.scamp' + WCS_TOL = time.strftime('%Y-%m-%d %H:%M:%S %p') + hdr.set('WCS_S', '0', '0=done') + hdr.set('WCS_V', WCS_V, 'Version of WCS calibration') + hdr.set('WCS_P', WCS_P, 'Configure file name of WCS') + hdr.set('WCS_TOL', WCS_TOL, 'Time of last wcs calibration') + hdul1.flush() + hdul1.close() + + def make_plots(self, ccdoff): + print('##### Analyzing the scampe result, making some pltos.... ####') + plt.figure(figsize=(11,5)) + ax1=plt.subplot(121) + bin=0.05 + plt.grid(color='grey',ls='--') + plt.plot(ccdoff[0],ccdoff[1],'ko',markersize=3,alpha=0.3) + plt.xlabel(r'$\Delta$ RA (arcsec)',fontsize=12) + plt.ylabel(r'$\Delta$ Dec (arcsec)',fontsize=12) + ax2=plt.subplot(122) + plt.grid(color='grey',ls='--') + plt.hist(ccdoff[0],bins=np.arange(-1,1, bin),histtype="step",color="r",label=r'$\Delta$RA (arcsec)') + plt.hist(ccdoff[1],bins=np.arange(-1,1, bin),histtype="step",color="b",label=r'$\Delta$Dec (arcsec)') + plt.legend() + a=str(float(ccdoff[2])) + b=str(float(ccdoff[4])) + c=str(float(ccdoff[3])) + d=str(float(ccdoff[5])) + plt.text(-0.95,45,r'$\mu$='+a[0:6]+r', $\sigma$='+b[0:5]+' (arcsec)',color='red') + plt.text(-0.95,35,r'$\mu$='+c[0:6]+r', $\sigma$='+d[0:5]+' (arcsec)',color='blue') + plt.xlabel('coord_diff (arcsec)',fontsize=12) + plt.savefig(self.dm.pc_radecoff,dpi=300) + + + def prepare(self, dm): + self.dm = dm + + def run(self, img_list, wht_list, flg_list, ccd_ids, path_gaia, search_radius, pm_correct=True, wcs_refine=False, plot=False): print('preparing files for position calibration....') - self.join_data(img_list, wht_list, flg_list, path_output=path_output) + self.join_data(img_list, wht_list, flg_list) print('################## run sextractor ###################') p = Pool() - prod_x = partial(self.run_sextractor, path_output=path_output) - result = p.map(prod_x, fn_list) + prod_x = partial(self.run_sextractor) + result = p.map(prod_x, ccd_ids) p.close() p.join() print('################## sextractor done ###################') print('############### combine sextractor catalog ###############') - self.combine_catalog(img_list, path_output) + self.combine_catalog(ccd_ids) print('############### get reference catalog ###############3') - refcat = self.get_refcat(img_list, path_gaia=path_gaia, search_radius=search_radius, silent=True) - Popen('cp ' + refcat + ' ref.cat', shell=True) + self.get_refcat(img_list, path_gaia=path_gaia, search_radius=search_radius, silent=True,pm_correct=True) print('############### run scamp ##################') - self.run_scamp(img_list, path_output=path_output) + self.run_scamp(wcs_refine=wcs_refine) print('################ scamp done #################') print('Checking astrometry quality....') - self.check_astrometry(img_list, path_output) + ccdoff = self.check_astrometry(img_list) print('################ updating headers.... #############') - self.write_headers(img_list) + self.write_headers(ccd_ids) print('#### Position calibration process done ####') - def cleanup(self, img_list, path_output): - # clean up environment - image_prefix = img_list[0][0].header['FILENAME'][0:-7] - for i in range(0, len(img_list)): - fn = img_list[i][0].header['FILENAME'] + '.acat' - if os.path.isfile(path_output + fn): os.remove(path_output + fn) - if os.path.isfile(path_output + image_prefix + '.gaia.fits'): - os.remove(path_output + image_prefix + '.gaia.fits') - if os.path.isfile(path_output + image_prefix + '.gaialac.fits'): - os.remove(path_output + image_prefix + '.gaialac.fits') - if os.path.isfile(path_output + 'scamp.xml'): - os.remove(path_output + 'scamp.xml') - if os.path.isfile(path_output + 'full_1.cat'): - os.remove(path_output + 'full_1.cat') - if os.path.isfile(path_output + 'merged_1.cat'): - os.remove(path_output + 'merged_1.cat') - if os.path.isfile(path_output + image_prefix + '_img.fits.back'): - os.remove( path_output + image_prefix + '_img.fits.back') - if os.path.isfile(path_output + image_prefix + '_wht.fits'): - os.remove(path_output + image_prefix + '_wht.fits') - if os.path.isfile(path_output + image_prefix + '_flg.fits'): - os.remove(path_output + image_prefix + '_flg.fits') + if plot: + make_plot(ccdoff) + + def cleanup(self): + pass + diff --git a/csst/msc/calib_pos.py~ b/csst/msc/calib_pos.py~ new file mode 100644 index 0000000000000000000000000000000000000000..ca68aabf42a6ec5c496970c186fe359227aa5763 --- /dev/null +++ b/csst/msc/calib_pos.py~ @@ -0,0 +1,465 @@ +import os +import time +from functools import partial +from multiprocessing import Pool +from subprocess import Popen +import healpy as hp +import numpy as np +from astropy import table +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS + +from .. import PACKAGE_PATH +from ..core.processor import CsstProcessor +from csst.msc.data_manager import CsstMscDataManager + +CONFIG_PATH = PACKAGE_PATH + "/msc/pos_calib_config/" + + +class CsstProcMscPositionCalibration(CsstProcessor): + + def join_data(self, img_list, wht_list, flg_list): + """ + Prepare data for running scamp; Combine all image data, weight files, flag files to their one frame. + + Parameters + ---------- + img_list: + image files to join together, e.g.,MSC_210304093000_0000000_06_img.fits + wht_list: + weith files to join together, e.g.,MSC_210304093000_0000000_06_img.fits + flg_list: + flag files to join together, e.g.,e.g.,MSC_210304093000_0000000_06_flg.fits + Returns + ------- + The joined multi-extension file(not stacked), weight, flag files. + e.g., combined_img.fits,combined_wht.fits, combined_flg.fits. + + """ + + output_imgnm = self.dm.pc_combined_file("img","fits") + hdul_img = fits.HDUList() + for i in range(0, len(img_list)): + h0 = fits.PrimaryHDU(header=img_list[i][0].header) + h1 = fits.ImageHDU(data=img_list[i][1].data, header=img_list[i][1].header) + hdul_img.append(h0) + hdul_img.append(h1) + hdul_img.writeto(output_imgnm, overwrite=True) + + output_whtnm = self.dm.pc_combined_file("wht","fits") + hdul_wht = fits.HDUList() + for i in range(0, len(wht_list)): + h0 = fits.PrimaryHDU(header=wht_list[i][0].header) + h1 = fits.ImageHDU(data=wht_list[i][1].data, header=wht_list[i][1].header) + hdul_wht.append(h0) + hdul_wht.append(h1) + hdul_wht.writeto(output_whtnm, overwrite=True) + + output_flgnm = self.dm.pc_combined_file("flg","fits") + hdul_flg = fits.HDUList() + for i in range(0, len(flg_list)): + h0 = fits.PrimaryHDU(header=flg_list[i][0].header) + h1 = fits.ImageHDU(data=flg_list[i][1].data, header=flg_list[i][1].header) + hdul_flg.append(h0) + hdul_flg.append(h1) + hdul_flg.writeto(output_flgnm, overwrite=True) + + def run_sextractor(self,ccd_ids): + """ + Run sextractor + + Parameters + ---------- + ccd_ids: + ccd list... + Returns + ------- + The photometric catalog, with position and flux, e.g.,MSC_210304093000_0000000_06_img.cat + """ + cat_list = [] + fn_list = [] + for this_ccd_id in ccd_ids: + sex_cat = self.dm.l1_sci(ccd_id = this_ccd_id, suffix="img", ext="cat") + fn = self.dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="fits") + cat_list.append(sex_cat) + fn_list.append(fn) + + config_sextractor = CONFIG_PATH + "new_csst_realtime.no.weight.sex" + sex_comd1 = 'sex -c ' + config_sextractor + ' ' + sex_comd2 = fn_list + ' -CATALOG_NAME ' + cat_list + sex_comd3 = ' -PARAMETERS_NAME ' + CONFIG_PATH + 'csst_realtime.param' + ' -FILTER_NAME ' + CONFIG_PATH + 'csst_realtime.conv' + ' -STARNNW_NAME ' + CONFIG_PATH + 'csst_realtime.nnw' + sex_comd = sex_comd1 + sex_comd2 + sex_comd3 + print(sex_comd) + p = Popen(sex_comd, shell=True) + p.wait() + + def combine_catalog(self, ccd_ids): + """ + Combine the sextractor catalog together + + Parameters + ----------- + ccd_ids: ccd list + Returns + ------- + The combined catalog,e.g., combined_cat.fits + """ + fn = self.dm.pc_combined_file("cat","fits") + hdul = fits.HDUList() + for this_ccd_id in ccd_ids: + cat_nm = self.dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="cat") + cat_i = fits.open(cat_nm) + data=cat_i[2].data + msk=(data['MAGERR_AUTO']<0.2)&(data['CLASS_STAR']>0.5)&(data['CLASS_STAR']<=1.0) + data=data[msk] + cat_i[2].data=data + hdul.append(cat_i[0]) + hdul.append(cat_i[1]) + hdul.append(cat_i[2]) + hdul.writeto(fn, overwrite=True) + + def run_scamp(self, wcs_refine=False): + """ + Run scamp + Returns + ------- + Image header updated with WCS keywords, combined_cat.head. + """ + cb_fn = self.dm.pc_combined_file("cat","fits") + config_scamp = CONFIG_PATH + "csst.scamp" + scamp_comd = 'scamp ' + cb_fn + ' -ASTREFCAT_NAME= ' + 'ref.cat\ + -MERGEDOUTCAT_NAME ' + 'merged.cat -FULLOUTCAT_NAME ' + 'full.cat\ + -c ' + config_scamp + print(scamp_comd) + p = Popen(scamp_comd, shell=True) + p.wait() + if wcs_refine: + Popen('cp '+ self.dm.pc_combined_file("cat","head")+" "+self.dm.pc_combined_file("cat","ahead"),shell=True) + scamp_comd = 'scamp ' + cb_fn + ' -ASTREFCAT_NAME= ' + 'ref.cat\ + -MERGEDOUTCAT_NAME ' + 'merged.cat -FULLOUTCAT_NAME ' + 'full.cat\ + -c ' + config_scamp + '-AHEADER_SUFFIX '+ self.dm.pc_combined_file("cat","ahead") + print(scamp_comd) + p=Popen(scamp_comd,shell=True) + p.wait() + + + def convert_hdu_to_ldac(self, hdu): + """ + Convert an hdu table to a fits_ldac table (format used by astromatic suite) + + Parameters + ---------- + hdu: + `astropy.io.fits.BinTableHDU` or `astropy.io.fits.TableHDU` + HDUList to convert to fits_ldac HDUList + + Returns + ------- + tbl1: + `astropy.io.fits.BinTableHDU` + Header info for fits table (LDAC_IMHEAD) + tbl2: + `astropy.io.fits.BinTableHDU` + Data table (LDAC_OBJECTS) + """ + tblhdr = np.array([hdu[1].header.tostring()]) + col1 = fits.Column(name='Field Header Card', array=tblhdr, format='13200A') + cols = fits.ColDefs([col1]) + tbl1 = fits.BinTableHDU.from_columns(cols) + tbl1.header['TDIM1'] = '(80, {0})'.format(len(hdu[1].header)) + tbl1.header['EXTNAME'] = 'LDAC_IMHEAD' + + dcol = fits.ColDefs(hdu[1].data) + tbl2 = fits.BinTableHDU.from_columns(dcol) + tbl2.header['EXTNAME'] = 'LDAC_OBJECTS' + return tbl1, tbl2 + + def get_refcat(self, img_list, path_gaia, search_radius, silent=True,pm_correct=True): + """ + Get reference catalog for scamp. The reference cat is GAIA EDR3. + + Parameters + ---------- + image_prefix: + a image to get its reference catalog, e.g.,combined_img.fits. + Usually the center of the image is the wcs parameters CRVAL1,CRVAL1. + search_radius: + circle radius for searching, units: degree. e.g., 2 degree for a 1x1 deg^2 image. + For large ccd size, use larger radius. csst, r=3 deg. + path_gaia: directory of the reference catalog. + + Returns + ------- + outcat: + filename of the cross matched catalog: ref.cat + """ + outcat = self.dm.pc_ref_cat + fnmae = cb_fn = self.dm.pc_combined_file("img","fits") + hdu = fits.open(fname) + header1 = hdu[0].header + header2 = hdu[1].header + if (pm_correct): + oday=header1['DATE-OBS'] + otime=header1['TIME-OBS'] + exptime=header1['EXPTIME'] + odaytime=header1['DATE-OBS']+'T'+header1['TIME-OBS'] + t = Time(oday, format='isot', scale='utc') + t.format = 'decimalyear' + t1=Time(2016.0, format='decimalyear', scale='utc') + deltatime= t.value -2016.0 + print('(Time - 2016.0) = deltatime =', deltatime) + else: + deltatime=0.00 + ra = float(header2['CRVAL1']) + dec = float(header2['CRVAL2']) + c = SkyCoord(ra, dec, unit=(u.deg, u.deg)) + print('ra, dec ra.deg dec.deg= ', ra, dec, c.ra.deg, c.dec.deg) + ra = c.ra.deg + dec = c.dec.deg + c = SkyCoord(ra, dec, unit=(u.deg, u.deg)) + phi = c.ra.deg / (180. / np.pi) + theta = (90. - c.dec.deg) / (180 / np.pi) + vec = hp.ang2vec(theta, phi) + list1 = hp.query_disc(nside=32, vec=vec, radius=np.radians(search_radius)) + if -1 in list1: + list1.remove(-1) + ipring = np.array(list1) + pix = np.unique(ipring) + npix = pix.size + print(ipring, 'ipring', pix, 'pix', npix, 'npix') + dt = np.dtype( + [('ra', '>f8'), ('dec', '>f8'), ('ra_error', '>f8'), ('dec_error', '>f8'), ('phot_g_mean_mag', '>f8'), + ('pmra', '>f8'), ('pmra_error', '>f8'), ('pmdec', '>f8'), ('pmdec_error', '>f8')]) + refcat = table.Table(dtype=dt) + for i in pix: + print('i= %5.5d' % i, path_gaia) + fname = path_gaia + 'healpix-' + '%5.5d' % i + '.fits' + print('fname=', fname) + if not silent: print('Reading ', fname) + d = table.Table.read(fname) + refcat = [refcat, d] + refcat = table.vstack(refcat, join_type='inner') + refcat.rename_column('ra', 'X_WORLD') + refcat.rename_column('dec', 'Y_WORLD') + + mask = (refcat['pmdec'] != refcat['pmdec']) + refcat['pmdec'][mask] = 0 + mask = (refcat['pmra'] != refcat['pmra']) + refcat['pmra'][mask] = 0 + refcat['X_WORLD'] = refcat['X_WORLD'] + deltatime * refcat['pmra'] / np.cos( + refcat['Y_WORLD'] / 180. * np.pi) / 3600.0 / 1000.0 + refcat['Y_WORLD'] = refcat['Y_WORLD'] + deltatime * refcat['pmdec'] / 3600.0 / 1000.0 + refcat['ra_error'] = refcat['ra_error'] / 1000.0 / 3600.0 + refcat['dec_error'] = refcat['dec_error'] / 1000.0 / 3600.0 + refcat.rename_column('ra_error', 'ERRA_WORLD') + refcat.rename_column('dec_error', 'ERRB_WORLD') + refcat.rename_column('phot_g_mean_mag', 'MAG') + + hdu = refcat + hdu1 = self.convert_hdu_to_ldac(hdu) + hdup = fits.PrimaryHDU() + hdu = hdu1[0] + tbhdu = hdu1[1] + thdulist = fits.HDUList([hdup, hdu, tbhdu]) + thdulist.writeto(outcat) + + print('##################### end #####################') + + def rewrite_wcs_head(self, head): + """ + Rewrite the WCS head from Scamp to the standard fits header + + Parameters + ---------- + head: scamp head file names + + Returns + ------- + wcshead: a new head file in fits format + + """ + wcshead = head + '.fits' + f = open(head, 'r') + f1 = open(wcshead, 'w') + a = '' + i = 0 + for v in f.readlines(): + sp = '' + asp = '' + i += 1 + if len(v) <= 81: + sp = ' ' * (81 - len(v)) + if 'END' in v: + asp = ' ' * 80 * (36 - i % 36) + i = i + (36 - i % 36) + # print(i) + a = a + v + sp + asp + f1.write(a.replace('\n', '')) + f1.close() + f.close() + return wcshead + + def check_astrometry(self, img_list): + """ + Check position calibration quality + + Parameters + ------------ + img_list: + list of images, in table format + + Returns + ------- + ccdraoff: + ra difference between observed catalog and reference catalog + ccddecoff: + dec difference between observed catalog and reference catalog + ccdraoff_med,ccddecoff_med: + median of the ra or dec difference + ccdraoff_rms,ccddecoff_rms: + rms of the ra or dec difference + """ + print('############## check the astrometry quality and save files ################') + r1 = [] + d1 = [] + fn_head = self.dm.pc_combined_file("cat","head") + fn_scat = self.dm.pc_combined_file("cat","fits") + data_scat = fits.open(fn_scat) + + wcshead = self.rewrite_wcs_head(fn_head) + scat = fits.open(fn_scat) + hdul = fits.HDUList() + for i in range(0, len(img_list)): + wcshdr = fits.getheader(wcshead, i, ignore_missing_simple=True) # read headers and change to RA---TPV,DEC--TPV for wcs_transfer package + wcshdr['CTYPE1'] = 'RA---TPV' + wcshdr['CTYPE2'] = 'DEC--TPV' + w = WCS(wcshdr) + sexcat = data_scat[i*3+2].data + x = sexcat['XWIN_IMAGE'] + y = sexcat['YWIN_IMAGE'] + r, d = w.all_pix2world(x, y, 0) # convert xwin,ywin to ra,de + sexcat['ALPHA_J2000'] = r + sexcat['DELTA_J2000'] = d + r1 = np.hstack((r1, r)) + d1 = np.hstack((d1, d)) + obsc = SkyCoord(ra=r1 * u.degree, dec=d1 * u.degree) + tmp_cat = np.zeros((len(obsc), 2)) + tmp_cat[:, 0] = obsc.ra + tmp_cat[:, 1] = obsc.dec + np.savetxt(dm.pc_scamp_coord, tmp_cat, fmt="%.10f %.10f", delimiter="\n") + + gaia_cat=Table.read(self.dm.pc_ref_cat) + gaia_ra=gaia_cat['X_WORLD'] + gaia_dec=gaia_cat['Y_WORLD'] + refc=SkyCoord(ra=gaia_ra*u.degree,dec=gaia_dec*u.degree) + idx, d2d, d3d = obsc.match_to_catalog_sky(refc) + ref_uid=np.unique(idx) + obs_uid=np.full_like(ref_uid,-1) + tmpj=-1 + ccdraoff_med=ccddecoff_med=ccdra_rms=ccddec_rms=-1 + for i in ref_uid: + tmpj=tmpj+1 + iid=(idx == i) + iiid = (d2d.deg[iid] == d2d.deg[iid].min()) + obs_uid[tmpj]=iid.nonzero()[0][iiid.nonzero()[0]][0] + + uidlim=d2d[obs_uid].arcsecond <1. # set match radius=1 arcsec + if uidlim.sum()>0: + obs_uidlim=obs_uid[uidlim] + ref_uidlim=ref_uid[uidlim] + ccdraoff=(obsc[obs_uidlim].ra- refc[ref_uidlim].ra).arcsec*np.cos(obsc[obs_uidlim].dec.deg*np.pi/180.) + ccdraoff_med=np.median(ccdraoff) + ccdra_rms=np.std(ccdraoff) + ccddecoff=(obsc[obs_uidlim].dec- refc[ref_uidlim].dec).arcsec + ccddec_rms=np.std(ccddecoff) + ccddecoff_med=np.median(ccddecoff) + match_num=len(ccdraoff) + print('################# astrometry result: ##############') + if (match_num<100): + print('### bad astrometry ###') + print('median ra_off, dec_off (mas) from scamp:',ccdraoff_med*1000.,ccddecoff_med*1000.) + print('rms ra_off, dec_off (mas) from scamp:',ccdra_rms*1000.,ccddec_rms*1000.) + print('############################################') + return ccdraoff,ccddecoff,ccdraoff_med,ccddecoff_med,ccdra_rms,ccddec_rms + + def write_headers(self, ccd_ids): + """ + Wrtie history to header + """ + head_suffix = self.dm.pc_combined_head_fits + hdul2 = fits.open(head_suffix, ignore_missing_simple=True) + for this_ccd_id in ccd_ids: + fits_nm = self.dm.l1_sci(this_ccd_id, suffix="img", ext="head") + hdul1 = fits.open(fits_nm, mode='update', ignore_missing_simple=True) + hdr = hdul1[0].header + hdr2 = hdul2[i].header + hdr.extend(hdr2, unique=True, update=True) + WCS_S = 0 + WCS_V = '2.0.4' + WCS_P = 'default.scamp' + WCS_TOL = time.strftime('%Y-%m-%d %H:%M:%S %p') + hdr.set('WCS_S', '0', '0=done') + hdr.set('WCS_V', WCS_V, 'Version of WCS calibration') + hdr.set('WCS_P', WCS_P, 'Configure file name of WCS') + hdr.set('WCS_TOL', WCS_TOL, 'Time of last wcs calibration') + hdul1.flush() + hdul1.close() + + def make_plots(self, ccdoff): + print('##### Analyzing the scampe result, making some pltos.... ####') + plt.figure(figsize=(11,5)) + ax1=plt.subplot(121) + bin=0.05 + plt.grid(color='grey',ls='--') + plt.plot(ccdoff[0],ccdoff[1],'ko',markersize=3,alpha=0.3) + plt.xlabel(r'$\Delta$ RA (arcsec)',fontsize=12) + plt.ylabel(r'$\Delta$ Dec (arcsec)',fontsize=12) + ax2=plt.subplot(122) + plt.grid(color='grey',ls='--') + plt.hist(ccdoff[0],bins=np.arange(-1,1, bin),histtype="step",color="r",label=r'$\Delta$RA (arcsec)') + plt.hist(ccdoff[1],bins=np.arange(-1,1, bin),histtype="step",color="b",label=r'$\Delta$Dec (arcsec)') + plt.legend() + a=str(float(ccdoff[2])) + b=str(float(ccdoff[4])) + c=str(float(ccdoff[3])) + d=str(float(ccdoff[5])) + plt.text(-0.95,45,r'$\mu$='+a[0:6]+r', $\sigma$='+b[0:5]+' (arcsec)',color='red') + plt.text(-0.95,35,r'$\mu$='+c[0:6]+r', $\sigma$='+d[0:5]+' (arcsec)',color='blue') + plt.xlabel('coord_diff (arcsec)',fontsize=12) + plt.savefig(self.dm.pc_radecoff,dpi=300) + + + def prepare(self, dm): + self.dm = dm + + def run(self, img_list, wht_list, flg_list, ccd_ids, path_gaia, search_radius, pm_correct=True, wcs_refine=False, plot=False): + print('preparing files for position calibration....') + self.join_data(img_list, wht_list, flg_list) + print('################## run sextractor ###################') + p = Pool() + prod_x = partial(self.run_sextractor) + result = p.map(prod_x, ccd_ids) + p.close() + p.join() + print('################## sextractor done ###################') + print('############### combine sextractor catalog ###############') + self.combine_catalog(ccd_ids) + print('############### get reference catalog ###############3') + self.get_refcat(img_list, path_gaia=path_gaia, search_radius=search_radius, silent=True,pm_correct=True) + print('############### run scamp ##################') + self.run_scamp(wcs_refine=wcs_refine) + print('################ scamp done #################') + print('Checking astrometry quality....') + ccdoff = self.check_astrometry(img_list) + print('################ updating headers.... #############') + self.write_headers(ccd_ids) + print('#### Position calibration process done ####') + + if plot: + make_plot(ccdoff) + + def cleanup(self): + pass + diff --git a/csst/msc/data_manager.py~ b/csst/msc/data_manager.py~ new file mode 100644 index 0000000000000000000000000000000000000000..c0880aeb89f874cac4c1ae0009c83b1c22f24cbd --- /dev/null +++ b/csst/msc/data_manager.py~ @@ -0,0 +1,289 @@ +import glob +import re + +from astropy.io import fits + +from .backbone import VER_SIMS, CCD_FILTER_MAPPING + + +class CsstMscDataManager: + """ this class defines the file format of the input & output of CSST MSC L1 pipeline + + C3: + MSC_MS_210525220000_100000020_06_raw.fits + MSC_CRS_210525220000_100000020_06_raw.fits + MSC_210525120000_0000020_06.cat + + C5.1: + CSST_MSC_MS_SCI_20270810081950_20270810082220_100000100_06_L0_1.fits + CSST_MSC_MS_CRS_20270810081950_20270810082220_100000100_06_L0_1.fits + MSC_10000100_chip_06_filt_y.cat + MSC_10000100_chip_06_filt_y.log + + """ + + def __init__(self, ver_sim="C5.1", dir_l0="", dir_l1="", dir_pcref="", path_aux="", force_all_ccds=False): + assert ver_sim in VER_SIMS + + self.dir_l0 = dir_l0 + self.dir_l1 = dir_l1 + self.dir_pcref = dir_pcref + self.path_aux = path_aux + self.ver_sim = ver_sim + + fps_img = self.glob_image(dir_l0, ver_sim=ver_sim) + fps_cat = self.glob_cat(dir_l0, ver_sim=ver_sim) + + if force_all_ccds: + assert len(fps_img) == len(CCD_ID_LIST) + else: + assert len(fps_img) > 0 + + if ver_sim == "C3": + # get info + # print(re.split(r"[_.]", fps[0])) + self._instrument, self._survey, \ + self._exp_start, self._exp_id, \ + _ccd_id, self._l0_suffix, _ext = re.split(r"[_.]", fps_img[0]) + self._cat_id = re.split(r"[_.]", fps_cat[0])[1] + + self._exp_start = int(self._exp_start) + self._exp_id = int(self._exp_id) + + # available CCD IDs + self.available_ccd_ids = [int(re.split(r"[_.]", fp)[4]) for fp in fps_img] + self.available_ccd_ids.sort() + + elif ver_sim == "C5.1": + # get info + # print(re.split(r"[_.]", fps[0])) + self._telescope, self._instrument, self._survey, self._imagetype, \ + self._exp_start, self._exp_stop, self._exp_id, \ + _ccd_id, self._l0_suffix, self._version, _ext = re.split(r"[_.]", fps_img[0]) + self._cat_id = re.split(r"[_.]", fps_cat[0])[1] + + self._exp_start = int(self._exp_start) + self._exp_stop = int(self._exp_stop) + self._exp_id = int(self._exp_id) + + # available CCD IDs + self.available_ccd_ids = [int(re.split(r"[_.]", fp)[7]) for fp in fps_img] + self.available_ccd_ids.sort() + + @staticmethod + def glob_image(dir_l0, ver_sim="C5"): + """ glob files in L0 data directory """ + if ver_sim == "C3": + pattern = os.path.join(dir_l0, "MSC_MS_*_raw.fits") + elif ver_sim == "C5.1": + pattern = os.path.join(dir_l0, "CSST_MSC_MS_SCI_*.fits") + fps = glob.glob(pattern) + fps = [os.path.basename(fp) for fp in fps] + fps.sort() + + print("@DM.glob_dir: {} files found with pattern: {}".format(len(fps), pattern)) + return fps + + @staticmethod + def glob_cat(dir_l0, ver_sim="C5"): + """ glob input catalogs in L0 data directory """ + if ver_sim == "C3": + pattern = os.path.join(dir_l0, "MSC_*.cat") + elif ver_sim == "C5.1": + pattern = os.path.join(dir_l0, "MSC_*.cat") + fps = glob.glob(pattern) + fps = [os.path.basename(fp) for fp in fps] + fps.sort() + + print("@DM.glob_dir: {} files found with pattern: {}".format(len(fps), pattern)) + return fps + + def l0_cat(self, ccd_id=6): + """ the L0 cat file path""" + if self.ver_sim == "C3": + fn = "{}_{}_{:07d}_{:02d}.cat".format( + self._instrument, self._cat_id, self._exp_id - 100000000, ccd_id) + elif self.ver_sim == "C5.1": + fn = "{}_{}_chip_{:02d}_filt_{}.cat".format( + self._instrument, self._exp_id - 90000000, ccd_id, CCD_FILTER_MAPPING[ccd_id]) + return os.path.join(self.dir_l0, fn) + + def l0_log(self, ccd_id=6): + """ L0 log file path """ + if self.ver_sim == "C5.1": + fn = "{}_{}_chip_{:02d}_filt_{}.log".format( + self._instrument, self._exp_id - 90000000, ccd_id, CCD_FILTER_MAPPING[ccd_id]) + return os.path.join(self.dir_l0, fn) + + def l0_sci(self, ccd_id=6): + """ L0 image file path """ + if self.ver_sim == "C3": + fn = "{}_{}_{}_{}_{:02d}_raw.fits".format( + self._instrument, self._survey, self._exp_start, self._exp_id, ccd_id) + elif self.ver_sim == "C5.1": + fn = "{}_{}_{}_SCI_{}_{}_{}_{:02d}_L0_1.fits".format( + self._telescope, self._instrument, self._survey, + self._exp_start, self._exp_stop, self._exp_id, ccd_id) + return os.path.join(self.dir_l0, fn) + + def l0_crs(self, ccd_id=6): + """ L0 cosmic ray file path """ + if self.ver_sim == "C3": + fn = "{}_CRS_{}_{}_{:02d}_raw.fits".format( + self._instrument, self._exp_start, self._exp_id, ccd_id) + elif self.ver_sim == "C5.1": + fn = "{}_{}_{}_CRS_{}_{}_{}_{:02d}_L0_1.fits".format( + self._telescope, self._instrument, self._survey, + self._exp_start, self._exp_stop, self._exp_id, ccd_id) + return os.path.join(self.dir_l0, fn) + + def l1_sci(self, ccd_id=6, suffix="img_whead", ext="fits"): + """ generate L1 file path + + Parameters + ---------- + ccd_id: + CCD ID + suffix: + {"img", "wht", "flg", "img_L1", "wht_L1", "flg_L1", "whead"} + ext: + {"fits", "cat", "rcat"} + + Returns + ------- + L1 file path + + """ + if self.ver_sim == "C3": + fn = "{}_{}_{}_{}_{:02d}_{}.{}".format( + self._instrument, self._survey, + self._exp_start, self._exp_id, ccd_id, suffix, ext) + elif self.ver_sim == "C5.1": + fn = "{}_{}_{}_SCI_{}_{}_{}_{:02d}_{}.{}".format( + self._telescope, self._instrument, self._survey, + self._exp_start, self._exp_stop, self._exp_id, ccd_id, suffix, ext) + return os.path.join(self.dir_l1, fn) + + def pc_combined_file(self, suffix="img", ext="fits"): + """ combined images + + Parameters + ---------- + suffix: + {"img", "wht", "flg", "cat", "head"} + ext: + {"fits", "head", "ahead" } + + Returns + ------- + combined image path + + """ + if self.ver_sim == "C3": + fn = "{}_{}_{}_{}_{}.{}".format( + self._instrument, self._survey, + self._exp_start, self._exp_id, suffix, ext) + elif self.ver_sim == "C5.1": + fn = "{}_{}_{}_SCI_{}_{}_{}_{}.{}".format( + self._telescope, self._instrument, self._survey, + self._exp_start, self._exp_stop, self._exp_id, suffix, ext) + return os.path.join(self.dir_l1, fn) + + @property + def pc_ref_cat(self): + """ reference catalog """ + return os.path.join(self.dir_l1, "ref.cat") + + @property + def pc_radecoff(self): + """ plot coordinate diff between obs and ref """ + return os.path.join(self.dir_l1, "radec_off.png") + + @property + def pc_check_fits(self): + """ check fits """ + return os.path.join(self.dir_l1, "check.fits") + + @property + def pc_combined_head_fits(self): + """ combined head (converted to) fits """ + return os.path.join(self.dir_l1, "combined_cat.head.fits") + + @property + def pc_scamp_coord(self): + """ SCAMP coord """ + return os.path.join(self.dir_l1, "scamp_coord.txt") + + def get_ccd_ids(self, ccd_ids=None): + """ """ + if ccd_ids is None: + # default ccd_ids + ccd_ids = self.available_ccd_ids + else: + try: + # assert ccd_ids is a subset of available ccd_ids + assert set(ccd_ids).issubset(set(self.available_ccd_ids)) + except AssertionError as ae: + print("@DM: available CCD IDs are ", self.available_ccd_ids) + print("@DM: target CCD IDs are ", ccd_ids) + raise ae + return ccd_ids + + def get_bias(self, ccd_id=6): + fp = glob.glob(self.path_aux.format("CLB", ccd_id))[0] + return fits.getdata(fp) + + def get_dark(self, ccd_id=6): + fp = glob.glob(self.path_aux.format("CLD", ccd_id))[0] + return fits.getdata(fp) + + def get_flat(self, ccd_id=6): + fp = glob.glob(self.path_aux.format("CLF", ccd_id))[0] + return fits.getdata(fp) + + +if __name__ == "__main__": + # test C3 + import os + from csst.msc.data_manager import CsstMscDataManager + + dm = CsstMscDataManager( + ver_sim="C3", dir_l0="/data/L1Pipeline/msc/MSC_0000020", dir_l1="/data/L1Pipeline/msc/work") + print("----- L0 images -----") + print(dm.l0_sci(ccd_id=6)) + print(os.path.exists(dm.l0_sci(ccd_id=6))) + print("----- L0 crs -----") + print(dm.l0_crs(ccd_id=6)) + print(os.path.exists(dm.l0_sci(ccd_id=8))) + print("----- L0 input cat -----") + print(dm.l0_cat(8)) + print(os.path.exists(dm.l0_cat(ccd_id=8))) + print("----- available ccd_ids -----") + print(dm.available_ccd_ids) + print("----- L1 images -----") + print(dm.l1_sci(25, "img", "fits")) + + # test C5.1 + import os + from csst.msc.data_manager import CsstMscDataManager + from csst.msc.backbone import CCD_ID_LIST + + dm = CsstMscDataManager( + ver_sim="C5.1", dir_l0="/data/sim_data/MSC_0000100", dir_l1="/home/user/L1Pipeline/msc/work") + print("----- available ccd_ids -----") + print(dm.available_ccd_ids) + for ccd_id in dm.available_ccd_ids[:2]: + print("----- L0 images -----") + print(dm.l0_sci(ccd_id=ccd_id)) + print(os.path.exists(dm.l0_sci(ccd_id=ccd_id))) + print("----- L0 crs -----") + print(dm.l0_crs(ccd_id=ccd_id)) + print(os.path.exists(dm.l0_sci(ccd_id=ccd_id))) + print("----- L0 input cat -----") + print(dm.l0_cat(ccd_id=ccd_id)) + print(os.path.exists(dm.l0_cat(ccd_id=ccd_id))) + print("----- L0 input log -----") + print(dm.l0_log(ccd_id=ccd_id)) + print(os.path.exists(dm.l0_log(ccd_id=ccd_id))) + print("----- L1 images -----") + print(dm.l1_sci(ccd_id, "img", "fits")) diff --git a/csst/msc/pipeline.py b/csst/msc/pipeline.py index 01fd85a7f306e83439b19ae4b3efbc08bf1977b2..d19dd8ae1eecae8bbefb46d4832a9ea096f3f40a 100644 --- a/csst/msc/pipeline.py +++ b/csst/msc/pipeline.py @@ -107,8 +107,24 @@ def do_one_exposure(ver_sim="C5.1", dir_l0="", dir_l1="", dir_pcref="", path_aux # Step 2. Calibrate Position pcProc = CsstProcMscPositionCalibration() - pcProc.run(img_list, wht_list, flg_list, fn_list, dir_pcref, dir_l1, 2.0) - pcProc.cleanup(img_list, dir_l1) + pcProc.prepare(dm) + if (img_list): + pcProc.run(img_list, wht_list, flg_list, ccd_ids, dir_pcref, 2.0, pm_correct=True, wcs_refine=True, plot=True) + else: + for this_ccd_id in ccd_ids: + fp_img = dm.l1_sci(ccd_id=this_ccd_id, suffix="img",ext="fits") + fp_wht = dm.l1_sci(ccd_id=this_ccd_id, suffix="wht",ext="fits") + fp_flg = dm.l1_sci(ccd_id=this_ccd_id, suffix="flg",ext="fits") + img = CsstMscImgData.read(fp_img) + wht = CsstMscImgData.read(fp_wht) + flg = CsstMscImgData.read(fp_flg) + img_list.append(img) + wht_list.append(wht) + flg_list.append(flg) + pcProc.run(img_list, wht_list, flg_list, ccd_ids, dir_pcref, 2.0, pm_correct=True, wcs_refine=True, plot=True) + + pcProc.cleanup() + """ pcProc = CsstProcMscPositionCalibration() pcProc.prepare(dm) @@ -120,13 +136,13 @@ def do_one_exposure(ver_sim="C5.1", dir_l0="", dir_l1="", dir_pcref="", path_aux img_list, dir_l1 """ - # Step 3. Calibrate Flux - fcProc = CsstProcFluxCalibration() - # fcProc.prepare() - fcProc.run( - fn_list, img_list, wht_list, flg_list, wcsdir=dir_l1, L1dir=dir_l1, workdir=dir_l1, refdir=dir_l0, - addhead=True, morehead=False, plot=False, nodel=False, update=False, upcat=True) - fcProc.cleanup(fn_list, dir_l1) + ## Step 3. Calibrate Flux + #fcProc = CsstProcFluxCalibration() + ## fcProc.prepare() + #fcProc.run( + #fn_list, img_list, wht_list, flg_list, wcsdir=dir_l1, L1dir=dir_l1, workdir=dir_l1, refdir=dir_l0, + #addhead=True, morehead=False, plot=False, nodel=False, update=False, upcat=True) + #fcProc.cleanup(fn_list, dir_l1) """ fcProc = CsstProcFluxCalibration() diff --git a/csst/msc/pipeline.py~ b/csst/msc/pipeline.py~ new file mode 100644 index 0000000000000000000000000000000000000000..b57f8a4c2ab22d287c83700447436187234b0e4a --- /dev/null +++ b/csst/msc/pipeline.py~ @@ -0,0 +1,187 @@ +import glob +import os + +from csst.msc.backbone import CCD_ID_LIST, VER_SIMS +from csst.msc.calib_flux import CsstProcFluxCalibration +from csst.msc.calib_pos import CsstProcMscPositionCalibration +from csst.msc.data import CsstMscImgData +from csst.msc.inst_corr import CsstMscInstrumentProc +from csst.msc.phot import CsstMscPhotometryProc +from csst.msc.data_manager import CsstMscDataManager + + +CONFIG_DANDELION = dict( + # test and working directory + dir_l0="/home/csstpipeline/L1Pipeline/msc/MSC_0000020", + dir_l1="/home/csstpipeline/L1Pipeline/msc/work/", + # on Dandelion + path_aux="/home/csstpipeline/L1Pipeline/msc/ref/MSC_{}_*_{:02d}_combine.fits", + # gaia catalog directory (for position calibration) + dir_pcref="/home/csstpipeline/L1Pipeline/msc/gaia_dr3/", + # version of simulation data + ver_sim="C3", + # only 18 cores available in cloud machine from PMO + n_jobs=18, + # shut down backend multithreading + backend_multithreading=False +) + + +CONFIG_PMO = dict( + # test and working directory + dir_l0="/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C5/NGP_AstrometryON_shearOFF/MSC_0000100", + dir_l1="/home/user/L1Pipeline/msc/work/", + # on PMO + path_aux="/data/sim_data/MSC_0000100/ref/MSC_{}_*_{:02d}_combine.fits", + # gaia catalog directory (for position calibration) + dir_pcref="/home/user/L1Pipeline/msc/gaia_dr3/", + # version of simulation data + ver_sim="C5.1", + # only 18 cores available in cloud machine from PMO + n_jobs=18, + # shut down backend multithreading + backend_multithreading=False +) + + +def do_one_exposure(ver_sim="C5.1", dir_l0="", dir_l1="", dir_pcref="", path_aux="", ccd_ids=None, + n_jobs=18, backend_multithreading=False): + + # currently C3 and C5.1 are tested + try: + assert ver_sim in VER_SIMS + except AssertionError as ae: + print("Available options for ver_sim: ", VER_SIMS) + raise ae + + # shut down backend multi-threading + if not backend_multithreading: + # prohibit multi-threading in backend + os.environ["MKL_NUM_THREADS"] = '1' + os.environ["NUMEXPR_NUM_THREADS"] = '1' + os.environ["OMP_NUM_THREADS"] = '1' + + # set paths + dm = CsstMscDataManager(ver_sim=ver_sim, dir_l0=dir_l0, dir_l1=dir_l1, dir_pcref=dir_pcref, path_aux=path_aux) + # request for ccd_ids + ccd_ids = dm.get_ccd_ids(ccd_ids) + + # Step 1. Correct instrumental effect + os.chdir(dir_l1) + + img_list = [] + wht_list = [] + flg_list = [] + fn_list = [] + for this_ccd_id in ccd_ids: + print("processing CCD {}".format(this_ccd_id)) + fp_raw = dm.l0_sci(ccd_id=this_ccd_id) + + # read data with CsstMscImgData.read + raw = CsstMscImgData.read(fp_raw) + + # in the future, use get_* functions grab + bias = dm.get_bias(this_ccd_id) + dark = dm.get_dark(this_ccd_id) + flat = dm.get_flat(this_ccd_id) + + # initialize Instrument Processor + instProc = CsstMscInstrumentProc() + instProc.prepare(n_jobs=n_jobs) + img, wht, flg = instProc.run(raw, bias, dark, flat, ver_sim) + instProc.cleanup() + fp_img = img[0].header["FILENAME"] + '.fits' + + # save img, wht, flg to somewhere + img.writeto(dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="fits"), overwrite=True) + wht.writeto(dm.l1_sci(ccd_id=this_ccd_id, suffix="wht", ext="fits"), overwrite=True) + flg.writeto(dm.l1_sci(ccd_id=this_ccd_id, suffix="flg", ext="fits"), overwrite=True) + # save header + img[1].header.tofile(dm.l1_sci(ccd_id=this_ccd_id, suffix="img", ext="head"), overwrite=True) + + # append img, wht, flg list + img_list.append(img) + wht_list.append(wht) + flg_list.append(flg) + fn_list.append(fp_img) + + # Step 2. Calibrate Position + pcProc = CsstProcMscPositionCalibration() + pcProc.prepare(dm) + if (img_list): + pcProc.run(img_list, wht_list, flg_list, ccd_ids, dir_pcref, 2.0, pm_correct=True, wcs_refine=True, plot=True) + else: + for this_ccd_id in ccd_ids: + fp_img = dm.l1_sci(ccd_id=this_ccd_id, suffix="img",ext="fits") + fp_wht = dm.l1_sci(ccd_id=this_ccd_id, suffix="wht",ext="fits") + fp_flg = dm.l1_sci(ccd_id=this_ccd_id, suffix="flg",ext="fits") + img = CsstMscImgData.read(fp_img) + wht = CsstMscImgData.read(fp_wht) + flg = CsstMscImgData.read(fp_flg) + img_list.append(img) + wht_list.append(wht) + flg_list.append(flg) + pcProc.run(img_list, wht_list, flg_list, ccd_ids, dir_pcref, 2.0, pm_correct=True, wcs_refine=True, plot=True) + + + + pcProc.cleanup() + + """ + pcProc = CsstProcMscPositionCalibration() + pcProc.prepare(dm) + pcProc.run(2.0) + pcProc.cleanup() + + get these parameters from dm.l1_sci(*): + img_list, wht_list, flg_list, fn_list, dir_pcref, dir_l1 + img_list, dir_l1 + """ + + ## Step 3. Calibrate Flux + #fcProc = CsstProcFluxCalibration() + ## fcProc.prepare() + #fcProc.run( + #fn_list, img_list, wht_list, flg_list, wcsdir=dir_l1, L1dir=dir_l1, workdir=dir_l1, refdir=dir_l0, + #addhead=True, morehead=False, plot=False, nodel=False, update=False, upcat=True) + #fcProc.cleanup(fn_list, dir_l1) + + """ + fcProc = CsstProcFluxCalibration() + fcProc.prepare(dm) + fcProc.run(addhead=True, morehead=False, plot=False, nodel=False, update=False, upcat=True) + fcProc.cleanup() + + get these parameters from dm.l1_sci(*): + fn_list, img_list, wht_list, flg_list, wcsdir=dir_l1, L1dir=dir_l1, workdir=dir_l1, refdir=dir_l0, + fn_list, dir_l1 + """ + + # Step 4. Photometry + # ptProc = CsstMscPhotometryProc() + # ptProc.prepare() + # ptProc.run(fn_list, out_dir=dir_l1, n_jobs=n_jobs) + # ptProc.cleanup() + + return + + +if __name__ == "__main__": + # identify where you are + HOSTNAME = os.uname()[1] + # you have to run this pipeline in some well-defined servers + assert HOSTNAME in ["ubuntu", "Dandelion"] + # get config parameters + if HOSTNAME == "ubuntu": + config = CONFIG_PMO + elif HOSTNAME == "Dandelion": + config = CONFIG_DANDELION + else: + raise ValueError("HOSTNAME {} not known!".format(HOSTNAME)) + + # process this exposure + do_one_exposure(**config) + + for k, v in config.items(): + eval("{}=config[\"{}\"]".format(k, k)) + diff --git a/csst/msc/pos_calib_config/default2.scamp b/csst/msc/pos_calib_config/csst.scamp similarity index 97% rename from csst/msc/pos_calib_config/default2.scamp rename to csst/msc/pos_calib_config/csst.scamp index 74b59777ccea4283be16b7a768f688d5c2dece81..c3df80227ca2683350b82916875ec63abdf6f856 100644 --- a/csst/msc/pos_calib_config/default2.scamp +++ b/csst/msc/pos_calib_config/csst.scamp @@ -56,12 +56,12 @@ FIXFOCALPLANE_NMIN 1 # Min number of dets for FIX_FOCALPLANE #---------------------------- Cross-identification ---------------------------- -CROSSID_RADIUS 2.0 # Cross-id initial radius (arcsec) +CROSSID_RADIUS 1.0 # Cross-id initial radius (arcsec) #---------------------------- Astrometric solution ---------------------------- SOLVE_ASTROM Y # Compute astrometric solution (Y/N) ? -PROJECTION_TYPE TAN # SAME, TPV or TAN +PROJECTION_TYPE TPV # SAME, TPV or TAN ASTRINSTRU_KEY FILTER,QRUNID # FITS keyword(s) defining the astrom STABILITY_TYPE EXPOSURE # EXPOSURE, PRE-DISTORTED or INSTRUMENT CENTROID_KEYS XWIN_IMAGE,YWIN_IMAGE # Cat. parameters for centroiding @@ -69,8 +69,8 @@ CENTROIDERR_KEYS ERRAWIN_IMAGE,ERRBWIN_IMAGE,ERRTHETAWIN_IMAGE # Cat. params for centroid err ellipse DISTORT_KEYS XWIN_IMAGE,YWIN_IMAGE # Cat. parameters or FITS keywords DISTORT_GROUPS 1,1 # Polynom group for each context key -DISTORT_DEGREES 2 # Polynom degree for each group -FOCDISTORT_DEGREE 1 # Polynom degree for focal plane coords +DISTORT_DEGREES 3 # Polynom degree for each group +FOCDISTORT_DEGREE 3 # Polynom degree for focal plane coords ASTREF_WEIGHT 1000000.0 # Relative weight of ref.astrom.ca ASTRACCURACY_TYPE SIGMA-PIXEL # SIGMA-PIXEL, SIGMA-ARCSEC, # or TURBULENCE-ARCSEC diff --git a/csst/msc/pos_calib_config/csst.scamp~ b/csst/msc/pos_calib_config/csst.scamp~ new file mode 100644 index 0000000000000000000000000000000000000000..c3df80227ca2683350b82916875ec63abdf6f856 --- /dev/null +++ b/csst/msc/pos_calib_config/csst.scamp~ @@ -0,0 +1,142 @@ +# Default configuration file for SCAMP 2.0.4 +# EB 2016-11-04 +# + +#----------------------------- Field grouping --------------------------------- + +FGROUP_RADIUS 2.0 # Max dist (deg) between field groups + +#---------------------------- Reference catalogs ------------------------------ + +REF_SERVER cocat1.u-strasbg.fr # Internet addresses of catalog servers +REF_PORT 80 # Ports to connect to catalog servers +CDSCLIENT_EXEC aclient_cgi # CDSclient executable +ASTREF_CATALOG FILE # NONE, FILE, USNO-A1,USNO-A2,USNO-B1, + # GSC-1.3,GSC-2.2,GSC-2.3, + # TYCHO-2, UCAC-1,UCAC-2,UCAC-3,UCAC-4, + # NOMAD-1, PPMX, CMC-14, 2MASS, DENIS-3, + # SDSS-R3,SDSS-R5,SDSS-R6,SDSS-R7, + # SDSS-R8, SDSS-R9 +ASTREF_BAND DEFAULT # Photom. band for astr.ref.magnitudes + # or DEFAULT, BLUEST, or REDDEST +ASTREFCAT_NAME ref.cat # Local astrometric reference catalogs +ASTREFCENT_KEYS X_WORLD,Y_WORLD # Local ref.cat. centroid parameters +ASTREFERR_KEYS ERRA_WORLD, ERRB_WORLD , ERRTHETA_WORLD + # Local ref.cat. err. ellipse params +ASTREFMAG_KEY MAG # Local ref.cat. magnitude parameter +ASTREFMAGERR_KEY MAGERR # Local ref.cat. mag. error parameter +ASTREFOBSDATE_KEY OBSDATE # Local ref.cat. obs. date parameter +ASTREFMAG_LIMITS -99,99 # Select magnitude range in ASTREF_BAND +SAVE_REFCATALOG Y # Save ref catalogs in FITS-LDAC format? +REFOUT_CATPATH ./ # Save path for reference catalogs + +#--------------------------- Merged output catalogs --------------------------- + +MERGEDOUTCAT_TYPE FITS_LDAC # NONE, ASCII_HEAD, ASCII, FITS_LDAC +MERGEDOUTCAT_NAME merged.cat # Merged output catalog filename + +#--------------------------- Full output catalogs --------------------------- + +FULLOUTCAT_TYPE FITS_LDAC # NONE, ASCII_HEAD, ASCII, FITS_LDAC +FULLOUTCAT_NAME full.cat # Full output catalog filename + +#----------------------------- Pattern matching ------------------------------- + +MATCH Y # Do pattern-matching (Y/N) ? +MATCH_NMAX 0 # Max.number of detections for MATCHing + # (0=auto) +PIXSCALE_MAXERR 2.0 # Max scale-factor uncertainty +POSANGLE_MAXERR 2.0 # Max position-angle uncertainty (deg) +POSITION_MAXERR 2.0 # Max positional uncertainty (arcmin) +MATCH_RESOL 0 # Matching resolution (arcsec); 0=auto +MATCH_FLIPPED Y # Allow matching with flipped axes? +MOSAIC_TYPE FIX_FOCALPLANE # UNCHANGED, SAME_CRVAL, SHARE_PROJAXIS, + # FIX_FOCALPLANE or LOOSE +FIXFOCALPLANE_NMIN 1 # Min number of dets for FIX_FOCALPLANE + +#---------------------------- Cross-identification ---------------------------- + +CROSSID_RADIUS 1.0 # Cross-id initial radius (arcsec) + +#---------------------------- Astrometric solution ---------------------------- + +SOLVE_ASTROM Y # Compute astrometric solution (Y/N) ? +PROJECTION_TYPE TPV # SAME, TPV or TAN +ASTRINSTRU_KEY FILTER,QRUNID # FITS keyword(s) defining the astrom +STABILITY_TYPE EXPOSURE # EXPOSURE, PRE-DISTORTED or INSTRUMENT +CENTROID_KEYS XWIN_IMAGE,YWIN_IMAGE # Cat. parameters for centroiding +CENTROIDERR_KEYS ERRAWIN_IMAGE,ERRBWIN_IMAGE,ERRTHETAWIN_IMAGE + # Cat. params for centroid err ellipse +DISTORT_KEYS XWIN_IMAGE,YWIN_IMAGE # Cat. parameters or FITS keywords +DISTORT_GROUPS 1,1 # Polynom group for each context key +DISTORT_DEGREES 3 # Polynom degree for each group +FOCDISTORT_DEGREE 3 # Polynom degree for focal plane coords +ASTREF_WEIGHT 1000000.0 # Relative weight of ref.astrom.ca +ASTRACCURACY_TYPE SIGMA-PIXEL # SIGMA-PIXEL, SIGMA-ARCSEC, + # or TURBULENCE-ARCSEC +ASTRACCURACY_KEY ASTRACCU # FITS keyword for ASTR_ACCURACY param. +ASTR_ACCURACY 0.00000001 # Astrom. uncertainty floor parameter +ASTRCLIP_NSIGMA 3 # Astrom. clipping threshold in sigmas +COMPUTE_PARALLAXES N # Compute trigonom. parallaxes (Y/N)? +COMPUTE_PROPERMOTIONS N # Compute proper motions (Y/N)? +CORRECT_COLOURSHIFTS N # Correct for colour shifts (Y/N)? +INCLUDE_ASTREFCATALOG N # Include ref.cat in prop.motions (Y/N)? +ASTR_FLAGSMASK 0x00fc # Astrometry rejection mask on SEx FLAGS +ASTR_IMAFLAGSMASK 0x0 # Astrometry rejection mask on IMAFLAGS + +#---------------------------- Photometric solution ---------------------------- + +SOLVE_PHOTOM Y # Compute photometric solution (Y/N) ? +MAGZERO_OUT 0.0 # Magnitude zero-point(s) in output +MAGZERO_INTERR 0.01 # Internal mag.zero-point accuracy +MAGZERO_REFERR 0.03 # Photom.field mag.zero-point accuracy +PHOTINSTRU_KEY FILTER # FITS keyword(s) defining the photom. +MAGZERO_KEY PHOT_C # FITS keyword for the mag zero-point +EXPOTIME_KEY EXPTIME # FITS keyword for the exposure time (s) +AIRMASS_KEY AIRMASS # FITS keyword for the airmass +EXTINCT_KEY PHOT_K # FITS keyword for the extinction coeff +PHOTOMFLAG_KEY PHOTFLAG # FITS keyword for the photometry flag +PHOTFLUX_KEY FLUX_AUTO # Catalog param. for the flux measurement +PHOTFLUXERR_KEY FLUXERR_AUTO # Catalog parameter for the flux error +PHOTCLIP_NSIGMA 3.0 # Photom.clipping threshold in sigmas +PHOT_ACCURACY 1e-3 # Photometric uncertainty floor (frac.) +PHOT_FLAGSMASK 0x00fc # Photometry rejection mask on SEx FLAGS +PHOT_IMAFLAGSMASK 0x0 # Photometry rejection mask on IMAFLAGS + +#------------------------------- Check-plots ---------------------------------- + +CHECKPLOT_CKEY SCAMPCOL # FITS keyword for PLPLOT field colour +CHECKPLOT_DEV PS # NULL, XWIN, TK, PS, PSC, XFIG, PNG, + # JPEG, AQT, PDF or SVG +CHECKPLOT_RES 0 # Check-plot resolution (0 = default) +CHECKPLOT_ANTIALIAS Y # Anti-aliasing using convert (Y/N) ? +CHECKPLOT_TYPE FGROUPS,DISTORTION,ASTR_INTERROR2D,ASTR_INTERROR1D,ASTR_REFERROR2D,ASTR_REFERROR1D,ASTR_CHI2,PHOT_ERROR +CHECKPLOT_NAME fgroups,distort,astr_interror2d,astr_interror1d,astr_referror2d,astr_referror1d,astr_chi2,psphot_error # Check-plot filename(s) + +#------------------------------- Check-images --------------------------------- + +CHECKIMAGE_TYPE AS_PAIR # NONE, AS_PAIR, AS_REFPAIR, or AS_XCORR +CHECKIMAGE_NAME check.fits # Check-image filename(s) + +#------------------------------ Miscellaneous --------------------------------- + +SN_THRESHOLDS 10.0,100.0 # S/N thresholds (in sigmas) for all and + # high-SN sample +FWHM_THRESHOLDS 0.0,100.0 # FWHM thresholds (in pixels) for sources +ELLIPTICITY_MAX 0.5 # Max. source ellipticity +FLAGS_MASK 0x00f0 # Global rejection mask on SEx FLAGS +WEIGHTFLAGS_MASK 0x00ff # Global rejec. mask on SEx FLAGS_WEIGHT +IMAFLAGS_MASK 0x0 # Global rejec. mask on SEx IMAFLAGS_ISO +AHEADER_GLOBAL scamp.abcdhead # Filename of the global INPUT header +AHEADER_SUFFIX .ahead # Filename extension for additional + # INPUT headers +HEADER_SUFFIX .head # Filename extension for OUTPUT headers +HEADER_TYPE NORMAL # NORMAL or FOCAL_PLANE +VERBOSE_TYPE NORMAL # QUIET, NORMAL, LOG or FULL +WRITE_XML N # Write XML file (Y/N)? +XML_NAME scamp.xml # Filename for XML output +XSL_URL file:///usr/local/share/scamp/scamp.xsl + # Filename for XSL style-sheet +NTHREADS 0 # Number of simultaneous threads for + # the SMP version of SCAMP + # 0 = automatic