updated combined_file

csst/msc/calib_pos.py~

+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
+

csst/msc/data_manager.py~

+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"))

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()

csst/msc/pipeline.py~

+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))
+

csst/msc/pos_calib_config/default2.scamp → 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

csst/msc/pos_calib_config/csst.scamp~

+# 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