From 1f8e630633e43adc3ddee5e310a05024b873c02a Mon Sep 17 00:00:00 2001
From: Bo Zhang <bozhang@nao.cas.cn>
Date: Wed, 11 May 2022 18:03:35 +0800
Subject: [PATCH] reformated flux calibration module
---
csst/msc/calib_flux.py | 1272 ++++++++++++++++++++--------------------
1 file changed, 631 insertions(+), 641 deletions(-)
diff --git a/csst/msc/calib_flux.py b/csst/msc/calib_flux.py
index 578c9d7..d1b0b9b 100644
--- a/csst/msc/calib_flux.py
+++ b/csst/msc/calib_flux.py
@@ -1,318 +1,312 @@
-#!/usr/bin/env python3.8
-#!/usr/bin/python3.8
-#Written by ZZM
-#Edited on Jun. 17, 2016
-#add color term
-__author__ = 'ZZM'
-
-import sys
import os
+import time
+from subprocess import Popen
-import shutil
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.stats import sigma_clip
from astropy.stats import sigma_clipped_stats
-from astropy.io import fits
from astropy.wcs import WCS
-from astropy import table
-from multiprocessing import Pool
-from subprocess import Popen, PIPE
-from astropy.wcs.utils import proj_plane_pixel_scales
-import time
+from astropy.wcs.utils import proj_plane_pixel_scales
from .. import PACKAGE_PATH
from ..core.processor import CsstProcessor
-#---------------------------------------------------------------------------
+
+# Edited on Jun. 17, 2016
+# add color term
+__author__ = 'ZZM'
+
+
class CsstProcFluxCalibration(CsstProcessor):
- def __init__(self):
- pass
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
- def ps1_mags(self,cat,band,obs_region):
+ def ps1_mags(self, cat, band, obs_region):
# get ps1 color-corrected magnitude for specific filter
- #ref:https://desi.lbl.gov/trac/wiki/DecamLegacy/Reductions
- color_gi = cat['MEDIAN'][:,0] - cat['MEDIAN'][:,2]
- coeff=np.array([0.0,0.0,0.0,0.0])
+ # ref:https://desi.lbl.gov/trac/wiki/DecamLegacy/Reductions
+ color_gi = cat['MEDIAN'][:, 0] - cat['MEDIAN'][:, 2]
+ coeff = np.array([0.0, 0.0, 0.0, 0.0])
if band == 'g':
ifilt = 0
- magcut=21.5
- coeff=np.array([0.00225738,-0.02543153,0.10575001,0.01077462])
+ magcut = 21.5
+ coeff = np.array([0.00225738, -0.02543153, 0.10575001, 0.01077462])
elif band == 'bokr' or band == 'r':
ifilt = 1
- magcut=20.0
- coeff=np.array([-0.00966711,0.02808938,-0.07650502,0.0024329])
+ magcut = 20.0
+ coeff = np.array([-0.00966711, 0.02808938, -0.07650502, 0.0024329])
elif band == 'z':
- magcut=20.0
+ magcut = 20.0
ifilt = 3
- coeff=np.array([-0.0094723,0.03040996,-0.08774486,0.01029133])
- #coeff=coeff=np.array([0.0,0.0,0.0,0.0])
+ coeff = np.array([-0.0094723, 0.03040996, -0.08774486, 0.01029133])
+ # coeff=coeff=np.array([0.0,0.0,0.0,0.0])
else:
- print(('no filter:',band))
+ print(('no filter:', band))
return
- colorterm= coeff[0]*color_gi**3 + coeff[1]*color_gi**2 + coeff[2]*color_gi + coeff[3]
- mags=cat['MEDIAN'][:,ifilt] + colorterm
-
- #select good phot_mag
- goodid= (cat['NMAG_OK'][:,0] >1) & (cat['NMAG_OK'][:,2] >1) \
- & (cat['NMAG_OK'][:,ifilt] >1) & (cat['STDEV'][:,ifilt] <0.1)\
- & (color_gi <2.7) & (color_gi >0.4) \
- & (cat['MEDIAN'][:,ifilt] <magcut) & (cat['MEDIAN'][:,ifilt] >16)\
- & (cat['RA']>obs_region[0]) & (cat['RA']<obs_region[1])\
- & (cat['DEC']>obs_region[2]) & (cat['DEC']<obs_region[3]) #2016.06.03
- print('starPS1.sum()=',goodid.sum())
- if goodid.sum() < 20: goodid=list(range(mags.size))
+ colorterm = coeff[0] * color_gi ** 3 + coeff[1] * color_gi ** 2 + coeff[2] * color_gi + coeff[3]
+ mags = cat['MEDIAN'][:, ifilt] + colorterm
+
+ # select good phot_mag
+ goodid = (cat['NMAG_OK'][:, 0] > 1) & (cat['NMAG_OK'][:, 2] > 1) \
+ & (cat['NMAG_OK'][:, ifilt] > 1) & (cat['STDEV'][:, ifilt] < 0.1) \
+ & (color_gi < 2.7) & (color_gi > 0.4) \
+ & (cat['MEDIAN'][:, ifilt] < magcut) & (cat['MEDIAN'][:, ifilt] > 16) \
+ & (cat['RA'] > obs_region[0]) & (cat['RA'] < obs_region[1]) \
+ & (cat['DEC'] > obs_region[2]) & (cat['DEC'] < obs_region[3]) # 2016.06.03
+ print('starPS1.sum()=', goodid.sum())
+ if goodid.sum() < 20: goodid = list(range(mags.size))
mags = mags[goodid]
- magerr=cat['STDEV'][:,ifilt][goodid]
+ magerr = cat['STDEV'][:, ifilt][goodid]
ra = cat['RA'][goodid]
dec = cat['DEC'][goodid]
- brmedian=np.median(color_gi[goodid])
-
- return mags,magerr,ra,dec,brmedian
-
- #----------------------------------------------------------------------
- def read_ps1cat(self,ra,dec,outcat,path='/line12/Pan-STARRS/chunks-qz-star-v2/',silent=True):
- #Read a piece of the PS1 calibration catalog containing ra,dec
+ brmedian = np.median(color_gi[goodid])
+
+ return mags, magerr, ra, dec, brmedian
+
+ def read_ps1cat(self, ra, dec, outcat, path='/line12/Pan-STARRS/chunks-qz-star-v2/', silent=True):
+ # Read a piece of the PS1 calibration catalog containing ra,dec
from healpy import ang2pix
- c=SkyCoord(ra,dec, unit=(u.deg, u.deg))
- if isinstance(ra[0],str):
- c=SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
- phi = c.ra.deg/(180./np.pi)
- theta = (90.-c.dec.deg)/(180/np.pi)
- ipring=ang2pix(32,theta,phi)
- pix=np.unique(ipring)
- npix=pix.size
- if npix >8:
- print('too many healpix files:',npix,' check image wcs!')
+ c = SkyCoord(ra, dec, unit=(u.deg, u.deg))
+ if isinstance(ra[0], str):
+ c = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
+ phi = c.ra.deg / (180. / np.pi)
+ theta = (90. - c.dec.deg) / (180 / np.pi)
+ ipring = ang2pix(32, theta, phi)
+ pix = np.unique(ipring)
+ npix = pix.size
+ if npix > 8:
+ print('too many healpix files:', npix, ' check image wcs!')
return
- dt=np.dtype([('OBJ_ID', '>i8'), ('RA', '>f8'), ('DEC', '>f8'), ('NMAG_OK', '>i2', (5,)),\
- ('STDEV', '>f4', (5,)), ('MEAN', '>f4', (5,)), ('MEDIAN', '>f4', (5,)), \
- ('MEAN_AP', '>f4', (5,)), ('MEDIAN_AP', '>f4', (5,))])
- ps=table.Table(dtype=dt)
+ dt = np.dtype([('OBJ_ID', '>i8'), ('RA', '>f8'), ('DEC', '>f8'), ('NMAG_OK', '>i2', (5,)), \
+ ('STDEV', '>f4', (5,)), ('MEAN', '>f4', (5,)), ('MEDIAN', '>f4', (5,)), \
+ ('MEAN_AP', '>f4', (5,)), ('MEDIAN_AP', '>f4', (5,))])
+ ps = table.Table(dtype=dt)
for i in pix:
- fname = path+'ps1-'+'%5.5d'%i+'.fits'
+ fname = path + 'ps1-' + '%5.5d' % i + '.fits'
if not silent: print(('Reading ', fname))
if os.path.isfile(fname):
- d=table.Table.read(fname)
- ps=[ps,d]
- ps=table.vstack(ps)
- if outcat: ps.write(outcat,format='fits',overwrite=True)
+ d = table.Table.read(fname)
+ ps = [ps, d]
+ ps = table.vstack(ps)
+ if outcat: ps.write(outcat, format='fits', overwrite=True)
return ps
- #-----------------------------------------------------------------------
- def read_gaiacat(self,ra,dec,outcat,path='./',silent=True):
+
+ def read_gaiacat(self, ra, dec, outcat, path='./', silent=True):
# function: Read a piece of the GAIA calibration catalog containing ra,dec
from healpy import ang2pix
- c=SkyCoord(ra,dec, unit=(u.deg, u.deg))
- if isinstance(ra[0],str):
- c=SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
- phi = c.ra.deg/(180./np.pi)
- theta = (90.-c.dec.deg)/(180/np.pi)
- ipring=ang2pix(32,theta,phi)
- pix=np.unique(ipring)
- npix=pix.size
- if npix >8:
- print('too many healpix files:',npix,' check image wcs!')
+ c = SkyCoord(ra, dec, unit=(u.deg, u.deg))
+ if isinstance(ra[0], str):
+ c = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
+ phi = c.ra.deg / (180. / np.pi)
+ theta = (90. - c.dec.deg) / (180 / np.pi)
+ ipring = ang2pix(32, theta, phi)
+ pix = np.unique(ipring)
+ npix = pix.size
+ if npix > 8:
+ print('too many healpix files:', npix, ' check image wcs!')
return
- dt=np.dtype([('source_id', '>i8'), ('ra', '>f8'), ('dec', '>f8'), ('ra_error', '>f4'), ('dec_error', '>f4'), ('phot_g_mean_mag', '>f4'), ('phot_g_mean_flux_over_error', '>f4'), ('phot_g_n_obs', '>i2'), ('phot_bp_mean_mag', '>f4'), ('phot_bp_mean_flux_over_error', '>f4'), ('phot_bp_n_obs', '>i2'), ('phot_rp_mean_mag', '>f4'), ('phot_rp_mean_flux_over_error', '>f4'), ('phot_rp_n_obs', '>i2'), ('astrometric_weight_al', '>f4'), ('astrometric_n_obs_al', '>i2'), ('astrometric_n_good_obs_al', '>i2'), ('astrometric_excess_noise', '>f4'), ('astrometric_excess_noise_sig', '>f4'), ('duplicated_source', '?'), ('ref_epoch', '>f4'), ('parallax', '>f4'), ('parallax_error', '>f4'), ('pmra', '>f4'), ('pmra_error', '>f4'), ('pmdec', '>f4'), ('pmdec_error', '>f4'), ('phot_variable_flag', '?')])
- ps=table.Table(dtype=dt)
+ dt = np.dtype([('source_id', '>i8'), ('ra', '>f8'), ('dec', '>f8'), ('ra_error', '>f4'), ('dec_error', '>f4'),
+ ('phot_g_mean_mag', '>f4'), ('phot_g_mean_flux_over_error', '>f4'), ('phot_g_n_obs', '>i2'),
+ ('phot_bp_mean_mag', '>f4'), ('phot_bp_mean_flux_over_error', '>f4'), ('phot_bp_n_obs', '>i2'),
+ ('phot_rp_mean_mag', '>f4'), ('phot_rp_mean_flux_over_error', '>f4'), ('phot_rp_n_obs', '>i2'),
+ ('astrometric_weight_al', '>f4'), ('astrometric_n_obs_al', '>i2'),
+ ('astrometric_n_good_obs_al', '>i2'), ('astrometric_excess_noise', '>f4'),
+ ('astrometric_excess_noise_sig', '>f4'), ('duplicated_source', '?'), ('ref_epoch', '>f4'),
+ ('parallax', '>f4'), ('parallax_error', '>f4'), ('pmra', '>f4'), ('pmra_error', '>f4'),
+ ('pmdec', '>f4'), ('pmdec_error', '>f4'), ('phot_variable_flag', '?')])
+ ps = table.Table(dtype=dt)
for i in pix:
- fname = path+'chunk-'+'%5.5d'%i+'.fits'
+ fname = path + 'chunk-' + '%5.5d' % i + '.fits'
if not silent: print(('Reading ', fname))
if os.path.isfile(fname):
- d=table.Table.read(fname)
- ps=[ps,d]
- ps=table.vstack(ps)
- if outcat: ps.write(outcat,format='fits',overwrite=True)
+ d = table.Table.read(fname)
+ ps = [ps, d]
+ ps = table.vstack(ps)
+ if outcat: ps.write(outcat, format='fits', overwrite=True)
return ps
- ##########################
- def gaia_mags(self,cat,band,obs_region):
+ def gaia_mags(self, cat, band, obs_region):
# get Gaia color-corrected magnitude for specific filter
- #color transp. parameters
- filterlist=['N','u','g','r','i','z','y']
- #ffpar=np.array([[-0.12372707,0.87415984,2.68152378,0.46688053],\
+ # color transp. parameters
+ filterlist = ['N', 'u', 'g', 'r', 'i', 'z', 'y']
+ # ffpar=np.array([[-0.12372707,0.87415984,2.68152378,0.46688053],\
# [ 0.0715824,0.42018862,1.58311364,0.30184685],\
# [-0.02193249,0.22376644,0.51403047,-0.04977904],\
# [ 0.0850583, 0.0045054, -0.17416451, 0.02696644],\
# [ 0.03638692, 0.0513301, -0.58839035, 0.02328468],\
# [-0.14492398, 0.26053385, -0.82220974, -0.0115362 ],\
# [-0.27416685, 0.36465315, -0.89341154, -0.02812699]])
- ffpar=np.array([[-1.92165818e+00, 7.88689593e+00, -1.60028226e+00, 1.53590717e+00],\
- [-2.00544954e+00, 7.17304891e+00, -4.92617803e+00, 3.05285639e+00],\
- [ 2.36560800e-02, 1.78153510e-01, 4.17186700e-01, 2.06523985e+00],\
- [ 8.04794200e-02, -6.17626400e-02, -1.81644150e-01, 1.89622143e+00],\
- [ 2.01861000e-03, 1.22630040e-01, -6.96895960e-01, 9.65536332e-01],\
- [-4.75608000e-02, 1.86804250e-01, -9.25372450e-01, 2.67927720e+00],\
- [-7.93956600e-02, 2.28310630e-01, -1.02675552e+00, 1.76769826e+00]])
- coeff=ffpar[filterlist.index(band),:]
- #print('filterlist.index(band):', filterlist.index(band))
- color_br = cat['phot_bp_mean_mag'] - cat['phot_rp_mean_mag'] ##- 0.3161466496 #VEGAMAG -> AB
- colorterm= coeff[0]*color_br**3 + coeff[1]*color_br**2 + coeff[2]*color_br + coeff[3]
- gaia_G=cat['phot_g_mean_mag'] ## +0.1001113078 #VEGAMAG -> AB
- mags=gaia_G + colorterm
- magcut=19.0
-
- #select good phot_mag
- goodid= (cat['phot_g_n_obs'] >3) & (cat['phot_bp_n_obs'] >1) & (cat['phot_rp_n_obs'] >1)\
- & (color_br <2.5) & (color_br >0.5) \
- & (gaia_G <magcut) & (gaia_G >16)\
- & (cat['ra']>obs_region[0]) & (cat['ra']<obs_region[1])\
- & (cat['dec']>obs_region[2]) & (cat['dec']<obs_region[3])
- print('starGAIA.sum()=',goodid.sum())
- if goodid.sum() < 20: goodid=list(range(mags.size))
+ ffpar = np.array([[-1.92165818e+00, 7.88689593e+00, -1.60028226e+00, 1.53590717e+00], \
+ [-2.00544954e+00, 7.17304891e+00, -4.92617803e+00, 3.05285639e+00], \
+ [2.36560800e-02, 1.78153510e-01, 4.17186700e-01, 2.06523985e+00], \
+ [8.04794200e-02, -6.17626400e-02, -1.81644150e-01, 1.89622143e+00], \
+ [2.01861000e-03, 1.22630040e-01, -6.96895960e-01, 9.65536332e-01], \
+ [-4.75608000e-02, 1.86804250e-01, -9.25372450e-01, 2.67927720e+00], \
+ [-7.93956600e-02, 2.28310630e-01, -1.02675552e+00, 1.76769826e+00]])
+ coeff = ffpar[filterlist.index(band), :]
+ # print('filterlist.index(band):', filterlist.index(band))
+ color_br = cat['phot_bp_mean_mag'] - cat['phot_rp_mean_mag'] ##- 0.3161466496 #VEGAMAG -> AB
+ colorterm = coeff[0] * color_br ** 3 + coeff[1] * color_br ** 2 + coeff[2] * color_br + coeff[3]
+ gaia_G = cat['phot_g_mean_mag'] ## +0.1001113078 #VEGAMAG -> AB
+ mags = gaia_G + colorterm
+ magcut = 19.0
+
+ # select good phot_mag
+ goodid = (cat['phot_g_n_obs'] > 3) & (cat['phot_bp_n_obs'] > 1) & (cat['phot_rp_n_obs'] > 1) \
+ & (color_br < 2.5) & (color_br > 0.5) \
+ & (gaia_G < magcut) & (gaia_G > 16) \
+ & (cat['ra'] > obs_region[0]) & (cat['ra'] < obs_region[1]) \
+ & (cat['dec'] > obs_region[2]) & (cat['dec'] < obs_region[3])
+ print('starGAIA.sum()=', goodid.sum())
+ if goodid.sum() < 20: goodid = list(range(mags.size))
mags = mags[goodid]
- magerr= mags * 0.0
+ magerr = mags * 0.0
ra = cat['ra'][goodid]
dec = cat['dec'][goodid]
- gimedian=np.median(color_br[goodid])
- refc=SkyCoord(ra,dec, unit=(u.deg, u.deg),frame='fk5')
-
- return mags,magerr,gimedian,refc
- #########################################
-
- def read_inputcat(self,image,outcat='refcat.fits',refdir='',silent=True):
- #imname=os.path.split(image)[-1]
+ gimedian = np.median(color_br[goodid])
+ refc = SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='fk5')
+
+ return mags, magerr, gimedian, refc
+
+ def read_inputcat(self, image, outcat='refcat.fits', refdir='', silent=True):
+ # imname=os.path.split(image)[-1]
if refdir == '':
- refdir=os.path.split(image)[0]
- #catname='MSC_'+imname[7:]+'.cat'
- catid=image[image.rfind('_')-10:image.rfind('_')]
- catname='MSC_210525120000_'+catid+'.cat'
- cat=os.path.join(refdir,catname)
- data=table.Table.read(cat,format='ascii')
- data.write(outcat,format='fits',overwrite=True)
-
+ refdir = os.path.split(image)[0]
+ # catname='MSC_'+imname[7:]+'.cat'
+ catid = image[image.rfind('_') - 10:image.rfind('_')]
+ catname = 'MSC_210525120000_' + catid + '.cat'
+ cat = os.path.join(refdir, catname)
+ data = table.Table.read(cat, format='ascii')
+ data.write(outcat, format='fits', overwrite=True)
+
return data
- ###########################
-
- def input_mags(self,cat,usepixcood=False):
- goodid=(cat['mag']<20) & (cat['mag']>16)
- n=0
- while(n<4):
- goodid=(cat['mag']<20+n) & (cat['mag']>16)
- n=n+0.5
- if goodid.sum()>20:
- n=n+5
- mags=cat['mag'][goodid]
- magerr= mags * 0.0
+
+ def input_mags(self, cat, usepixcood=False):
+ goodid = (cat['mag'] < 20) & (cat['mag'] > 16)
+ n = 0
+ while (n < 4):
+ goodid = (cat['mag'] < 20 + n) & (cat['mag'] > 16)
+ n = n + 0.5
+ if goodid.sum() > 20:
+ n = n + 5
+ mags = cat['mag'][goodid]
+ magerr = mags * 0.0
if usepixcood:
- x=cat['xImage'][goodid]
- y=cat['yImage'][goodid]
- z=np.zeros_like(x)
- refc = SkyCoord(x=x,y=y,z=z,unit=u.pix,representation_type='cartesian')
+ x = cat['xImage'][goodid]
+ y = cat['yImage'][goodid]
+ z = np.zeros_like(x)
+ refc = SkyCoord(x=x, y=y, z=z, unit=u.pix, representation_type='cartesian')
else:
ra = cat['ra'][goodid]
dec = cat['dec'][goodid]
- refc=SkyCoord(ra,dec, unit=(u.deg, u.deg),frame='fk5')
- brmedian=0.0
-
- return mags,magerr,brmedian,refc
- #########################################
-
- def rewrite_sex_cat(self,cat,workdir=''):
- data=fits.open(cat)
+ refc = SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='fk5')
+ brmedian = 0.0
+
+ return mags, magerr, brmedian, refc
+
+ def rewrite_sex_cat(self, cat, workdir=''):
+ data = fits.open(cat)
if workdir == '':
- workdir=os.path.split(cat)[0]
- index=[]
- if (np.size(data)%3 == 0) & (np.size(data) !=0):
- for i in range(int(np.size(data)/3)):
- #print(i)
- hdr=data[3*i+1].data
- filename=str(hdr)[str(hdr).rfind('FITSFILE')+11:
- str(hdr).rfind(".fits' / File name of the analyse")]
- subfile=os.path.join(workdir,filename+'.acat')
+ workdir = os.path.split(cat)[0]
+ index = []
+ if (np.size(data) % 3 == 0) & (np.size(data) != 0):
+ for i in range(int(np.size(data) / 3)):
+ # print(i)
+ hdr = data[3 * i + 1].data
+ filename = str(hdr)[str(hdr).rfind('FITSFILE') + 11:
+ str(hdr).rfind(".fits' / File name of the analyse")]
+ subfile = os.path.join(workdir, filename + '.acat')
if not os.path.exists(subfile):
- hdu=fits.HDUList([data[0],data[3*i+1],data[3*i+2]])
- hdu.writeto(subfile,overwrite=True)
+ hdu = fits.HDUList([data[0], data[3 * i + 1], data[3 * i + 2]])
+ hdu.writeto(subfile, overwrite=True)
index.append(filename)
return index
- ##########################
- def split_wcs_head(self,wcshead,im_index=[],workdir=''):
+ def split_wcs_head(self, wcshead, im_index=[], workdir=''):
if len(im_index) == 0:
return
- head=fits.open(wcshead,ignore_missing_simple=True)
+ head = fits.open(wcshead, ignore_missing_simple=True)
if np.size(head) != np.size(im_index):
return
-
- for i,filename in enumerate(im_index):
- pathsplit=os.path.split(filename)
+
+ for i, filename in enumerate(im_index):
+ pathsplit = os.path.split(filename)
if workdir == '':
workdir = pathsplit[0]
- headname=os.path.join(workdir,pathsplit[1] + '.whead.fits')
- wheader=head[i].header
+ headname = os.path.join(workdir, pathsplit[1] + '.whead.fits')
+ wheader = head[i].header
prihdu = fits.PrimaryHDU(header=wheader)
- prihdu.writeto(headname,overwrite=True)
- ##########################
-
- def run_sextractor(self,image,cat):
- image0=image
- header0=fits.getheader(image, 0)
- header=fits.getheader(image, 1)
-
+ prihdu.writeto(headname, overwrite=True)
+
+ def run_sextractor(self, image, cat):
+ image0 = image
+ header0 = fits.getheader(image, 0)
+ header = fits.getheader(image, 1)
+
CONFIG_PATH = PACKAGE_PATH + '/msc/flux_calib_config/'
- #sexfile='/home/zhouzm/data/csst/code/cali_default.sex'
+ # sexfile='/home/zhouzm/data/csst/code/cali_default.sex'
sexfile = CONFIG_PATH + 'cali_default.sex'
sex_comd3 = ' -PARAMETERS_NAME ' + CONFIG_PATH + 'cali_default.param' + ' -FILTER_NAME ' + CONFIG_PATH + 'cali_default.conv' + ' -STARNNW_NAME ' + CONFIG_PATH + 'cali_default.nnw'
- exptime=header0['EXPTIME']
- gain=header['GAIN1']*exptime
- #saturate=header['SATURATE']/exptime
- sex_zp='0.0'
- imwht=image.replace('_img','_wht')
- #apture='4,9,15,21,27,36,39,44'
- apture='10'
- sex='sex '+ image+' -c '+ sexfile+' -CATALOG_NAME '+cat+\
- ' -PHOT_APERTURES '+apture+' -MAG_ZEROPOINT '+sex_zp+\
- ' -WEIGHT_IMAGE '+imwht +' -GAIN '+str(gain) +\
- ' -GAIN_KEY '+ 'abcdefg' + sex_comd3
- #sex WFST_000001_r_210103150723_9.fits -c cali_default.sex -CATALOG_NAME ./WFST_000001_r_210103150723_9_calibphot.fits -PHOT_APERTURES 36 -MAG_ZEROPOINT 0.0 -WEIGHT_IMAGE WFST_000001_r_210103150723_9.wht.fits -GAIN 2.047026888573 -GAIN_KEY abcdefg
+ exptime = header0['EXPTIME']
+ gain = header['GAIN1'] * exptime
+ # saturate=header['SATURATE']/exptime
+ sex_zp = '0.0'
+ imwht = image.replace('_img', '_wht')
+ # apture='4,9,15,21,27,36,39,44'
+ apture = '10'
+ sex = 'sex ' + image + ' -c ' + sexfile + ' -CATALOG_NAME ' + cat + \
+ ' -PHOT_APERTURES ' + apture + ' -MAG_ZEROPOINT ' + sex_zp + \
+ ' -WEIGHT_IMAGE ' + imwht + ' -GAIN ' + str(gain) + \
+ ' -GAIN_KEY ' + 'abcdefg' + sex_comd3
+ # sex WFST_000001_r_210103150723_9.fits -c cali_default.sex -CATALOG_NAME ./WFST_000001_r_210103150723_9_calibphot.fits -PHOT_APERTURES 36 -MAG_ZEROPOINT 0.0 -WEIGHT_IMAGE WFST_000001_r_210103150723_9.wht.fits -GAIN 2.047026888573 -GAIN_KEY abcdefg
print(sex)
- p=Popen(sex,shell=True)
+ p = Popen(sex, shell=True)
p.wait()
-
+
return cat
- ##########################
- def prepare(self,image,wcsdir,workdir,usewcsresult=False,newcat=False):
- #make calibration files: Obs. cat (SExtractor), PSF model (PSFex)
- #sex_zp=str(25.5+2.5*np.log10(exptime))
+ def prepare(self, image, wcsdir, workdir, usewcsresult=False, newcat=False):
+ # make calibration files: Obs. cat (SExtractor), PSF model (PSFex)
+ # sex_zp=str(25.5+2.5*np.log10(exptime))
if not os.path.exists(workdir):
os.mkdir(workdir)
- imag_file_ename=os.path.split(image)[-1]
- imname=imag_file_ename[0:imag_file_ename.find('.fits')]
- imname0=imag_file_ename[0:imag_file_ename.find('_img')-3]
- wcscat0=os.path.join(wcsdir,imname0+'.acat.fits')
- wcscat1=os.path.join(wcsdir,imname+'.acat')
- wcshead0=os.path.join(wcsdir,imname0+'.acat.head.fits')
- wcshead0txt=wcshead0[:wcshead0.rfind('.fits')]
- wcshead1=os.path.join(wcsdir,imname+'.acat.head.fits')
- cat=os.path.join(workdir,imname+'.acat')
- ref=os.path.join(workdir,imname+'.rcat')
- wcshead2=os.path.join(workdir,imname+'.whead.fits')
-
- #cali_ref='GAIA' #calibration reference data
-
- #run SExtractor
- im_index=[]
+ imag_file_ename = os.path.split(image)[-1]
+ imname = imag_file_ename[0:imag_file_ename.find('.fits')]
+ imname0 = imag_file_ename[0:imag_file_ename.find('_img') - 3]
+ wcscat0 = os.path.join(wcsdir, imname0 + '.acat.fits')
+ wcscat1 = os.path.join(wcsdir, imname + '.acat')
+ wcshead0 = os.path.join(wcsdir, imname0 + '.acat.head.fits')
+ wcshead0txt = wcshead0[:wcshead0.rfind('.fits')]
+ wcshead1 = os.path.join(wcsdir, imname + '.acat.head.fits')
+ cat = os.path.join(workdir, imname + '.acat')
+ ref = os.path.join(workdir, imname + '.rcat')
+ wcshead2 = os.path.join(workdir, imname + '.whead.fits')
+
+ # cali_ref='GAIA' #calibration reference data
+
+ # run SExtractor
+ im_index = []
if (os.path.isfile(cat)) & (newcat):
os.remove(cat)
if (not os.path.isfile(wcscat1)) or ((os.path.isfile(wcscat0)) & (usewcsresult)):
- im_index=self.rewrite_sex_cat(wcscat0,workdir)
+ im_index = self.rewrite_sex_cat(wcscat0, workdir)
if os.path.isfile(wcshead0):
- self.split_wcs_head(wcshead0,im_index,workdir)
+ self.split_wcs_head(wcshead0, im_index, workdir)
elif os.path.isfile(wcshead0txt):
- wcshead0=self.rewrite_wcs_head(wcshead0txt)
- self.split_wcs_head(wcshead0,im_index,workdir)
+ wcshead0 = self.rewrite_wcs_head(wcshead0txt)
+ self.split_wcs_head(wcshead0, im_index, workdir)
if (not os.path.isfile(wcscat1)):
- self.run_sextractor(image,wcscat1)
+ self.run_sextractor(image, wcscat1)
+
+ header = self.combine_head(image, wcshead1, wcshead2, prime=False)
- header=self.combine_head(image,wcshead1,wcshead2,prime=False)
-
- return wcscat1,cat,ref,header
- #-----------------------------------------------------------------------
+ return wcscat1, cat, ref, header
- #def getebv(image):
+ # def getebv(image):
# ebv=0.0
# header=fits.getheader(image, 0)
# if image[-3:]=='.fz':
@@ -325,518 +319,514 @@ class CsstProcFluxCalibration(CsstProcessor):
# if tid.sum()==1: ebv=tiles['EBV'][tid].data[0]
# #header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')
# return ebv
- #########################
- def getmlim(self,fwhm=0.15,avsky=7.0,rdnoise=5.0,zpt=25.8,ebv=0.0,filter='g'):
+ def getmlim(self, fwhm=0.15, avsky=7.0, rdnoise=5.0, zpt=25.8, ebv=0.0, filter='g'):
# E(B-V): 3.995, 3.214, 2.165, 1.592, 1.211, 1.064 for ugrizY decals
- #https://github.com/dstndstn/tractor/blob/master/tractor/sfd.py#L10
- #if filter=='g':k=3.214; pixsize=0.33
- k=0.0
- pixsize=0.075
- snr=5.0
- #nea=(((fwhm/2.35)**2*4*np.pi)**(1/1.15)+(8.91*(0.45/pixsize)**2)**(1/1.15))**1.15
- #noise=np.abs(avsky)*nea*exptime+7.3**2*nea
- #flim=0.5/exptime*(snr**2+np.sqrt(snr**4+4*snr**2*noise))
-
- noise=np.sqrt(avsky*np.pi*(fwhm/2.35)**2+rdnoise**2)
- flim=noise*snr
- mlim=-2.5*np.log10(flim)+zpt-k*ebv
+ # https://github.com/dstndstn/tractor/blob/master/tractor/sfd.py#L10
+ # if filter=='g':k=3.214; pixsize=0.33
+ k = 0.0
+ pixsize = 0.075
+ snr = 5.0
+ # nea=(((fwhm/2.35)**2*4*np.pi)**(1/1.15)+(8.91*(0.45/pixsize)**2)**(1/1.15))**1.15
+ # noise=np.abs(avsky)*nea*exptime+7.3**2*nea
+ # flim=0.5/exptime*(snr**2+np.sqrt(snr**4+4*snr**2*noise))
+
+ noise = np.sqrt(avsky * np.pi * (fwhm / 2.35) ** 2 + rdnoise ** 2)
+ flim = noise * snr
+ mlim = -2.5 * np.log10(flim) + zpt - k * ebv
return mlim
- #########################
- def getfwhm(self,fwhmsex,ellip,obsme,obsflags):
- nanid=np.isnan(fwhmsex)
+ def getfwhm(self, fwhmsex, ellip, obsme, obsflags):
+ nanid = np.isnan(fwhmsex)
if nanid.sum() > 0:
- fwhmsex[np.isnan(fwhmsex)]=0.0
- starid1= (ellip <=0.15) & (obsme <=0.15) &(obsme >0)& (obsflags<1) &(fwhmsex>0)
- starid2= fwhmsex[starid1] < 5.0 #temp. value, 20220314
- starid=starid1.nonzero()[0][starid2]
- if starid2.sum()<10:
- starid=list(range(fwhmsex.size))
-
- fwhmsub=fwhmsex[starid]
- fwhmid=sigma_clip(fwhmsub,sigma=2.5)
- f1,f2=np.histogram(fwhmid.data[~fwhmid.mask])
- fbin=np.where(f1==f1.max())
- f2l=f2[np.max([fbin[0][0]-1,0])]
- f2h=f2[np.min([fbin[0][0]+1,10])]
- id=np.where((fwhmid.data[~fwhmid.mask]>f2l) & (fwhmid.data[~fwhmid.mask]<f2h))
- fwhm=np.median(fwhmid.data[~fwhmid.mask][id])
+ fwhmsex[np.isnan(fwhmsex)] = 0.0
+ starid1 = (ellip <= 0.15) & (obsme <= 0.15) & (obsme > 0) & (obsflags < 1) & (fwhmsex > 0)
+ starid2 = fwhmsex[starid1] < 5.0 # temp. value, 20220314
+ starid = starid1.nonzero()[0][starid2]
+ if starid2.sum() < 10:
+ starid = list(range(fwhmsex.size))
+
+ fwhmsub = fwhmsex[starid]
+ fwhmid = sigma_clip(fwhmsub, sigma=2.5)
+ f1, f2 = np.histogram(fwhmid.data[~fwhmid.mask])
+ fbin = np.where(f1 == f1.max())
+ f2l = f2[np.max([fbin[0][0] - 1, 0])]
+ f2h = f2[np.min([fbin[0][0] + 1, 10])]
+ id = np.where((fwhmid.data[~fwhmid.mask] > f2l) & (fwhmid.data[~fwhmid.mask] < f2h))
+ fwhm = np.median(fwhmid.data[~fwhmid.mask][id])
if np.isnan(fwhm):
- fwhm=np.mean(fwhmsex)
- #print 'fwhm=',fwhm,f1,f2,fbin, f2l,f2h
- #print id,fwhmid.data[~fwhmid.mask][id]
- #print fwhmid, fwhmsex
- return fwhm
- #########################
-
- def makedatedir(self,path):
- index=path.rfind('/')
- rootpath=path[0:path.rfind('/',0,index)]
- walk=os.walk(rootpath)
- datedir=walk.next()[1]
+ fwhm = np.mean(fwhmsex)
+ # print 'fwhm=',fwhm,f1,f2,fbin, f2l,f2h
+ # print id,fwhmid.data[~fwhmid.mask][id]
+ # print fwhmid, fwhmsex
+ return fwhm
+
+ def makedatedir(self, path):
+ index = path.rfind('/')
+ rootpath = path[0:path.rfind('/', 0, index)]
+ walk = os.walk(rootpath)
+ datedir = walk.next()[1]
for i in datedir:
if not os.path.exists(i):
- os.mkdir(i)
- print(('makedir ', i))
- ##########################
+ os.mkdir(i)
+ print(('makedir ', i))
+ ##########################
- def rewrite_wcs_head(self,head):
- #rewrite the WCS head from Scamp to the standard fits header.
- wcshead=head+'.fits'
+ def rewrite_wcs_head(self, head):
+ # rewrite the WCS head from Scamp to the standard fits header.
+ 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))
+ 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',''))
+ 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 combine_head(self,image,wcshead1='',wcshead2='',prime=False):
- #combine image head and wcs head keywords
- inst_head_file=image[0:image.find('.fits')]+'.head'
+ def combine_head(self, image, wcshead1='', wcshead2='', prime=False):
+ # combine image head and wcs head keywords
+ inst_head_file = image[0:image.find('.fits')] + '.head'
if os.path.isfile(inst_head_file):
- h0=fits.getheader(inst_head_file,ignore_missing_simple=True)
+ h0 = fits.getheader(inst_head_file, ignore_missing_simple=True)
else:
- h0=fits.getheader(image,1,ignore_missing_simple=True)
+ h0 = fits.getheader(image, 1, ignore_missing_simple=True)
if os.path.isfile(wcshead1):
- hw=fits.getheader(wcshead1,ignore_missing_simple=True)
+ hw = fits.getheader(wcshead1, ignore_missing_simple=True)
elif os.path.isfile(wcshead2):
- hw=fits.getheader(wcshead2,ignore_missing_simple=True)
+ hw = fits.getheader(wcshead2, ignore_missing_simple=True)
if prime:
- hprime=fits.getheader(image,0,ignore_missing_simple=True)
- h0.extend(hprime,unique=True, update=True)
+ hprime = fits.getheader(image, 0, ignore_missing_simple=True)
+ h0.extend(hprime, unique=True, update=True)
try:
- h0.extend(hw,unique=True, update=True)
+ h0.extend(hw, unique=True, update=True)
except:
print(image, ': no wcs header file')
-
- #update the keywords CTYPE1 & CTYPE2 in header
+
+ # update the keywords CTYPE1 & CTYPE2 in header
if ('CTYPE1' in list(h0.keys())) & ('PV1_0' in list(h0.keys())):
- h0['CTYPE1'] = 'RA---TPV'
- h0['CTYPE2'] = 'DEC--TPV'
-
+ h0['CTYPE1'] = 'RA---TPV'
+ h0['CTYPE2'] = 'DEC--TPV'
+
return h0
- ##########################
- def match_calib(self,obsc,refc,obsm,refm,obsme,refme,obsflags,fwhmsex=np.array([])):
+ def match_calib(self, obsc, refc, obsm, refm, obsme, refme, obsflags, fwhmsex=np.array([])):
# matching obs and ref catalogs
'''
coeff0=self.match_calib(obsc,refc,obsm,refm)
obsc,refc: obs and ref coordinates
obsm,refm: obs and ref magnitudes
'''
- coeff0=cstd=cobsm=crefm=-1
- ccdraoff=ccddecoff=-1
- csize=0
- fwhm=-1
- if obsm.size >0 and refm.size>0:
- obsid=(obsme < 0.2) #0.2->0.1 #20170518
- obsc=obsc[obsid]
- obsm=obsm[obsid]
- obsme=obsme[obsid]
- obsflags=obsflags[obsid]
- idx, d2d,d3d=obsc.match_to_catalog_sky(refc)
- ref_uid=np.unique(idx)
- obs_uid=np.full_like(ref_uid,-1)
- tmpj=-1
+ coeff0 = cstd = cobsm = crefm = -1
+ ccdraoff = ccddecoff = -1
+ csize = 0
+ fwhm = -1
+ if obsm.size > 0 and refm.size > 0:
+ obsid = (obsme < 0.2) # 0.2->0.1 #20170518
+ obsc = obsc[obsid]
+ obsm = obsm[obsid]
+ obsme = obsme[obsid]
+ obsflags = obsflags[obsid]
+ idx, d2d, d3d = obsc.match_to_catalog_sky(refc)
+ ref_uid = np.unique(idx)
+ obs_uid = np.full_like(ref_uid, -1)
+ tmpj = -1
for i in ref_uid:
- tmpj=tmpj+1
- iid=(idx == i)
+ 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 <0.5 # set match radius=0.5 arcsec
- if uidlim.sum()>0:
- #csize=uidlim.sum()
- obs_uidlim=obs_uid[uidlim]
- ref_uidlim=ref_uid[uidlim]
-
- ccdraoff=np.median((obsc[obs_uidlim].ra- refc[ref_uidlim].ra).arcsec*np.cos(obsc[obs_uidlim].dec.deg*np.pi/180))
- ccddecoff=np.median((obsc[obs_uidlim].dec- refc[ref_uidlim].dec).arcsec)
- #if fwhmsex.size >1:
+ obs_uid[tmpj] = iid.nonzero()[0][iiid.nonzero()[0]][0]
+
+ uidlim = d2d[obs_uid].arcsecond < 0.5 # set match radius=0.5 arcsec
+ if uidlim.sum() > 0:
+ # csize=uidlim.sum()
+ obs_uidlim = obs_uid[uidlim]
+ ref_uidlim = ref_uid[uidlim]
+
+ ccdraoff = np.median(
+ (obsc[obs_uidlim].ra - refc[ref_uidlim].ra).arcsec * np.cos(obsc[obs_uidlim].dec.deg * np.pi / 180))
+ ccddecoff = np.median((obsc[obs_uidlim].dec - refc[ref_uidlim].dec).arcsec)
+ # if fwhmsex.size >1:
# fwhm=np.median(fwhmsex[obs_uidlim])
- #matched magnitude
-
- #select good photometric value
- goodphot= (obsme[obs_uidlim] < 0.1) & (obsflags[obs_uidlim] < 2 )
- csize=goodphot.sum()
- if goodphot.sum()<=5:
- goodphot= (obsme[obs_uidlim] < 0.2) & (obsflags[obs_uidlim] < 5 )
- if goodphot.sum()>=1:
- obsm=obsm[obs_uidlim][goodphot]
- refm=refm[ref_uidlim][goodphot]
- obsme=obsme[obs_uidlim][goodphot]
- refme=refme[ref_uidlim][goodphot]
- combme=np.sqrt(obsme**2+refme**2)
- clip=sigma_clip(refm-obsm,sigma=2.5) #sigma=2.5
- coeff0=np.median(clip.data[~clip.mask])
- #wht=1./combme
- #coeff0=np.average(refm-obsm,weights=wht)
- cstd=np.std(clip.data[~clip.mask])
- cobsm=obsm[~clip.mask]
- crefm=refm[~clip.mask]
- #if fwhmsex.size >1:
+ # matched magnitude
+
+ # select good photometric value
+ goodphot = (obsme[obs_uidlim] < 0.1) & (obsflags[obs_uidlim] < 2)
+ csize = goodphot.sum()
+ if goodphot.sum() <= 5:
+ goodphot = (obsme[obs_uidlim] < 0.2) & (obsflags[obs_uidlim] < 5)
+ if goodphot.sum() >= 1:
+ obsm = obsm[obs_uidlim][goodphot]
+ refm = refm[ref_uidlim][goodphot]
+ obsme = obsme[obs_uidlim][goodphot]
+ refme = refme[ref_uidlim][goodphot]
+ combme = np.sqrt(obsme ** 2 + refme ** 2)
+ clip = sigma_clip(refm - obsm, sigma=2.5) # sigma=2.5
+ coeff0 = np.median(clip.data[~clip.mask])
+ # wht=1./combme
+ # coeff0=np.average(refm-obsm,weights=wht)
+ cstd = np.std(clip.data[~clip.mask])
+ cobsm = obsm[~clip.mask]
+ crefm = refm[~clip.mask]
+ # if fwhmsex.size >1:
# fwhm=np.median(fwhmsex[obsid][obs_uidlim][goodphot])
- #print 'fwhm=',fwhm*0.26
- #print(coeff0,cstd,csize,cobsm,crefm,ccdraoff,ccddecoff)
- #print('crefm.min,max=',refm.min(),refm.max())
- return coeff0,cstd,csize,cobsm,crefm,ccdraoff,ccddecoff,fwhm
- ########################
-
- def calib(self,image,imgdata,whtdata,flgdata,wcsdir='./',L1dir='./',workdir='./',refdir='',addhead=False,morehead=True,plot=False,nodel=True,update=False,upcat=True):
-
- print('calibration for:',image)
+ # print 'fwhm=',fwhm*0.26
+ # print(coeff0,cstd,csize,cobsm,crefm,ccdraoff,ccddecoff)
+ # print('crefm.min,max=',refm.min(),refm.max())
+ return coeff0, cstd, csize, cobsm, crefm, ccdraoff, ccddecoff, fwhm
+
+ def calib(self, image, imgdata, whtdata, flgdata, wcsdir='./', L1dir='./', workdir='./', refdir='', addhead=False,
+ morehead=True, plot=False, nodel=True, update=False, upcat=True):
+
+ print('calibration for:', image)
if not os.path.exists(L1dir):
os.mkdir(L1dir)
if not os.path.exists(workdir):
os.mkdir(workdir)
-
- imag_file_ename=os.path.split(image)[-1]
- imname=imag_file_ename[0:imag_file_ename.rfind('.')]
- image_head_file=os.path.join(workdir,imname+'.head')
- #if os.path.exists(image_head_file) and (not update):
+
+ imag_file_ename = os.path.split(image)[-1]
+ imname = imag_file_ename[0:imag_file_ename.rfind('.')]
+ image_head_file = os.path.join(workdir, imname + '.head')
+ # if os.path.exists(image_head_file) and (not update):
# print((image+' has been calibrated.'))
# return
-
- psname =os.path.join(workdir,imname+'_calib.png')
+
+ psname = os.path.join(workdir, imname + '_calib.png')
if plot and os.path.isfile(psname) and (not update):
- print((psname,' is existed, pass'))
+ print((psname, ' is existed, pass'))
return
-
- ckf=os.path.join(workdir,'checkwcs.ls')
- ckim=os.path.join(workdir,'checkim.ls')
-
- #get astrometry head
- #header=self.combine_head(image)
-
- cali_ref='GAIA' #calibration reference data
- #Get the photometric catalog, astrometry head
- cat,newcat,ref,header=self.prepare(image,wcsdir,workdir,newcat=False)
-
- #get obs. Information from header
- #k=list(header.keys())
- w=WCS(header)
- band=header['FILTER'][0] #['CCDLABEL']
- #!#airmass=header['AIRMASS']
- #exptime=header['EXPTIME']
- exptime=150.0
- naxis1=9216 #header['NAXIS1']
- naxis2=9232 #header['NAXIS2']
- #gain=header['GAIN1']
- pixsize=np.mean(proj_plane_pixel_scales(w))*3600.0 #pixsize=0.33
- imr0,imd0=w.all_pix2world(naxis1/2.,naxis2/2.,0)
-
- #obs catalog: coord & magnitude
- catdata=fits.open(cat)
- sexcat=catdata[2].data
- #sexcat=table.Table.read(cat,2)
+
+ ckf = os.path.join(workdir, 'checkwcs.ls')
+ ckim = os.path.join(workdir, 'checkim.ls')
+
+ # get astrometry head
+ # header=self.combine_head(image)
+
+ cali_ref = 'GAIA' # calibration reference data
+ # Get the photometric catalog, astrometry head
+ cat, newcat, ref, header = self.prepare(image, wcsdir, workdir, newcat=False)
+
+ # get obs. Information from header
+ # k=list(header.keys())
+ w = WCS(header)
+ band = header['FILTER'][0] # ['CCDLABEL']
+ # !#airmass=header['AIRMASS']
+ # exptime=header['EXPTIME']
+ exptime = 150.0
+ naxis1 = 9216 # header['NAXIS1']
+ naxis2 = 9232 # header['NAXIS2']
+ # gain=header['GAIN1']
+ pixsize = np.mean(proj_plane_pixel_scales(w)) * 3600.0 # pixsize=0.33
+ imr0, imd0 = w.all_pix2world(naxis1 / 2., naxis2 / 2., 0)
+
+ # obs catalog: coord & magnitude
+ catdata = fits.open(cat)
+ sexcat = catdata[2].data
+ # sexcat=table.Table.read(cat,2)
if upcat:
- x=sexcat['XWIN_IMAGE']
- y=sexcat['YWIN_IMAGE']
- r,d=w.all_pix2world(x,y,0)
- sexcat['ALPHAWIN_J2000']=r
- sexcat['DELTAWIN_J2000']=d
- catdata.writeto(newcat,overwrite=True)
+ x = sexcat['XWIN_IMAGE']
+ y = sexcat['YWIN_IMAGE']
+ r, d = w.all_pix2world(x, y, 0)
+ sexcat['ALPHAWIN_J2000'] = r
+ sexcat['DELTAWIN_J2000'] = d
+ catdata.writeto(newcat, overwrite=True)
catdata.close()
-
- select_stars=(sexcat['MAG_APER']-sexcat['MAG_AUTO']<0.05)
- if select_stars.sum()>0:
- sexcat=sexcat[select_stars]
- r=sexcat['ALPHAWIN_J2000'] #ALPHA_J2000->ALPHAWIN_J2000
- d=sexcat['DELTAWIN_J2000']
- #back=sexcat['BACKGROUND']
- fwhmsex=sexcat['FWHM_IMAGE']
- obsflags=sexcat['FLAGS']
- cstar=sexcat['CLASS_STAR']
- ellip=sexcat['ELLIPTICITY']
- obsm=sexcat['MAG_APER'] #aperture radii=13pix
- obsme=sexcat['MAGERR_APER']
- ccdnstar=obsm.size #ccdnstar; total number of stars detected on a CCD
+
+ select_stars = (sexcat['MAG_APER'] - sexcat['MAG_AUTO'] < 0.05)
+ if select_stars.sum() > 0:
+ sexcat = sexcat[select_stars]
+ r = sexcat['ALPHAWIN_J2000'] # ALPHA_J2000->ALPHAWIN_J2000
+ d = sexcat['DELTAWIN_J2000']
+ # back=sexcat['BACKGROUND']
+ fwhmsex = sexcat['FWHM_IMAGE']
+ obsflags = sexcat['FLAGS']
+ cstar = sexcat['CLASS_STAR']
+ ellip = sexcat['ELLIPTICITY']
+ obsm = sexcat['MAG_APER'] # aperture radii=13pix
+ obsme = sexcat['MAGERR_APER']
+ ccdnstar = obsm.size # ccdnstar; total number of stars detected on a CCD
if cstar.size <= 0:
- #ckim=image[image.rfind('/')+1:image.rfind('.')][0:5]+'checkim.ls'
- imcheck= open(ckim,"a")
+ # ckim=image[image.rfind('/')+1:image.rfind('.')][0:5]+'checkim.ls'
+ imcheck = open(ckim, "a")
print('line501: image needs update. Pass.')
print(image, file=imcheck)
imcheck.close()
return
-
- #obs coordinate
- obsc = SkyCoord(r,d, unit=(u.deg, u.deg),frame='fk5')
- obsrmin=r.min()
- obsrmax=r.max()
- obsdmin=d.min()
- obsdmax=d.max()
- obs_region=[obsrmin,obsrmax,obsdmin,obsdmax]
- #measure image seeing:
- seeing=-1
- fwhm=-1
- fwhm=self.getfwhm(fwhmsex,ellip,obsme,obsflags)
- seeing=fwhm*pixsize
-
- #get reference catalog
- #ps=self.read_gaiacat(r,d,outcat=ref,silent=True)
- ps=self.read_inputcat(image,outcat=ref,refdir=refdir)
- #remove tmp file: cat & ref
-
- if np.size(ps) <1:
- wcscheck= open(ckf,"a")
+
+ # obs coordinate
+ obsc = SkyCoord(r, d, unit=(u.deg, u.deg), frame='fk5')
+ obsrmin = r.min()
+ obsrmax = r.max()
+ obsdmin = d.min()
+ obsdmax = d.max()
+ obs_region = [obsrmin, obsrmax, obsdmin, obsdmax]
+ # measure image seeing:
+ seeing = -1
+ fwhm = -1
+ fwhm = self.getfwhm(fwhmsex, ellip, obsme, obsflags)
+ seeing = fwhm * pixsize
+
+ # get reference catalog
+ # ps=self.read_gaiacat(r,d,outcat=ref,silent=True)
+ ps = self.read_inputcat(image, outcat=ref, refdir=refdir)
+ # remove tmp file: cat & ref
+
+ if np.size(ps) < 1:
+ wcscheck = open(ckf, "a")
print(image, file=wcscheck)
wcscheck.close()
return
else:
- #refm,refme,gimedian,refc=self.gaia_mags(ps,band,obs_region)
- refm,refme,gimedian,refc=self.input_mags(ps,usepixcood=False)
-
- #4: different doors
- im1= (obsc.ra.deg >= imr0) & (obsc.dec.deg <= imd0)
- im2= (obsc.ra.deg < imr0) & (obsc.dec.deg <= imd0)
- im3= (obsc.ra.deg >= imr0) & (obsc.dec.deg > imd0)
- im4= (obsc.ra.deg < imr0) & (obsc.dec.deg > imd0)
- ref1= (refc.ra.deg >= imr0) & (refc.dec.deg <= imd0)
- ref2= (refc.ra.deg < imr0) & (refc.dec.deg <= imd0)
- ref3= (refc.ra.deg >= imr0) & (refc.dec.deg > imd0)
- ref4= (refc.ra.deg < imr0) & (refc.dec.deg > imd0)
- #for all chips
- coeff0=self.match_calib(obsc,refc,obsm,refm,obsme,refme,obsflags,fwhmsex)
+ # refm,refme,gimedian,refc=self.gaia_mags(ps,band,obs_region)
+ refm, refme, gimedian, refc = self.input_mags(ps, usepixcood=False)
+
+ # 4: different doors
+ im1 = (obsc.ra.deg >= imr0) & (obsc.dec.deg <= imd0)
+ im2 = (obsc.ra.deg < imr0) & (obsc.dec.deg <= imd0)
+ im3 = (obsc.ra.deg >= imr0) & (obsc.dec.deg > imd0)
+ im4 = (obsc.ra.deg < imr0) & (obsc.dec.deg > imd0)
+ ref1 = (refc.ra.deg >= imr0) & (refc.dec.deg <= imd0)
+ ref2 = (refc.ra.deg < imr0) & (refc.dec.deg <= imd0)
+ ref3 = (refc.ra.deg >= imr0) & (refc.dec.deg > imd0)
+ ref4 = (refc.ra.deg < imr0) & (refc.dec.deg > imd0)
+ # for all chips
+ coeff0 = self.match_calib(obsc, refc, obsm, refm, obsme, refme, obsflags, fwhmsex)
if coeff0[2] <= 0:
- imcheck= open(ckim,"a")
+ imcheck = open(ckim, "a")
print('line560: image needs update. Pass.')
print(image, file=imcheck)
imcheck.close()
- #return
-
- coeff1=self.match_calib(obsc[im1],refc[ref1],obsm[im1],refm[ref1],obsme[im1],refme[ref1],obsflags[im1])
- coeff2=self.match_calib(obsc[im2],refc[ref2],obsm[im2],refm[ref2],obsme[im2],refme[ref2],obsflags[im2])
- coeff3=self.match_calib(obsc[im3],refc[ref3],obsm[im3],refm[ref3],obsme[im3],refme[ref3],obsflags[im3])
- coeff4=self.match_calib(obsc[im4],refc[ref4],obsm[im4],refm[ref4],obsme[im4],refme[ref4],obsflags[im4])
- coeff=[coeff0[0],coeff1[0],coeff2[0],coeff3[0],coeff4[0]]
- std=[coeff0[1],coeff1[1],coeff2[1],coeff3[1],coeff4[1]]
- match=[coeff0[2],coeff1[2],coeff2[2],coeff3[2],coeff4[2]]
- ccdraoff=coeff0[5]
- ccddecoff=coeff0[6]
-
- colortpar='0.0,0.0,0.0,0.0'
- aperture_radii=10
- aper_r_com='(pixels) photo-aperture radius'
- colortpar_com='par.in [br^3,br^2,br^1,br^0]'
-
- ccdzp_com='zero point for CCD'
- #ccdzpa_com='zero point for CCD ampA'
- #ccdzpb_com='zero point for CCD ampB'
- #ccdzpc_com='zero point for CCD ampC'
- #ccdzpd_com='zero point for CCD ampD'
- ccdphrms_com='zpt rms of the matched objects in CCD'
- #ccdphrms_coma='zpt rms of the matched objects in CCD ampA'
- #ccdphrms_comb='zpt rms of the matched objects in CCD ampB'
- #ccdphrms_comc='zpt rms of the matched objects in CCD ampC'
- #ccdphrms_comd='zpt rms of the matched objects in CCD ampD'
-
- fwhm_com='FWHM in pixel'
- seeing_com='seeing in arcsec'
- ccdraoff_com='median positional offset from GAIA, in arcsec'
- ccddecoff_com='median positional offset from GAIA, in arcsec'
- #transparency_com='median transparency.'#
- ccdnstar_com='total number of stars detected on a CCD'
- ccdnmatch_com='total number of matched stars in 2 arcsec'
- #ccdnmatcha_com='number of matched stars in CCD ampA'
- #ccdnmatchb_com='number of matched stars in CCD ampB'
- #ccdnmatchc_com='number of matched stars in CCD ampC'
- #ccdnmatchd_com='number of matched stars in CCD ampD'
- ccdmdncol_com='median (BP-RP)_GAIA of matched stars in CCD'
+ # return
+
+ coeff1 = self.match_calib(obsc[im1], refc[ref1], obsm[im1], refm[ref1], obsme[im1], refme[ref1], obsflags[im1])
+ coeff2 = self.match_calib(obsc[im2], refc[ref2], obsm[im2], refm[ref2], obsme[im2], refme[ref2], obsflags[im2])
+ coeff3 = self.match_calib(obsc[im3], refc[ref3], obsm[im3], refm[ref3], obsme[im3], refme[ref3], obsflags[im3])
+ coeff4 = self.match_calib(obsc[im4], refc[ref4], obsm[im4], refm[ref4], obsme[im4], refme[ref4], obsflags[im4])
+ coeff = [coeff0[0], coeff1[0], coeff2[0], coeff3[0], coeff4[0]]
+ std = [coeff0[1], coeff1[1], coeff2[1], coeff3[1], coeff4[1]]
+ match = [coeff0[2], coeff1[2], coeff2[2], coeff3[2], coeff4[2]]
+ ccdraoff = coeff0[5]
+ ccddecoff = coeff0[6]
+
+ colortpar = '0.0,0.0,0.0,0.0'
+ aperture_radii = 10
+ aper_r_com = '(pixels) photo-aperture radius'
+ colortpar_com = 'par.in [br^3,br^2,br^1,br^0]'
+
+ ccdzp_com = 'zero point for CCD'
+ # ccdzpa_com='zero point for CCD ampA'
+ # ccdzpb_com='zero point for CCD ampB'
+ # ccdzpc_com='zero point for CCD ampC'
+ # ccdzpd_com='zero point for CCD ampD'
+ ccdphrms_com = 'zpt rms of the matched objects in CCD'
+ # ccdphrms_coma='zpt rms of the matched objects in CCD ampA'
+ # ccdphrms_comb='zpt rms of the matched objects in CCD ampB'
+ # ccdphrms_comc='zpt rms of the matched objects in CCD ampC'
+ # ccdphrms_comd='zpt rms of the matched objects in CCD ampD'
+
+ fwhm_com = 'FWHM in pixel'
+ seeing_com = 'seeing in arcsec'
+ ccdraoff_com = 'median positional offset from GAIA, in arcsec'
+ ccddecoff_com = 'median positional offset from GAIA, in arcsec'
+ # transparency_com='median transparency.'#
+ ccdnstar_com = 'total number of stars detected on a CCD'
+ ccdnmatch_com = 'total number of matched stars in 2 arcsec'
+ # ccdnmatcha_com='number of matched stars in CCD ampA'
+ # ccdnmatchb_com='number of matched stars in CCD ampB'
+ # ccdnmatchc_com='number of matched stars in CCD ampC'
+ # ccdnmatchd_com='number of matched stars in CCD ampD'
+ ccdmdncol_com = 'median (BP-RP)_GAIA of matched stars in CCD'
cali_ref_com = 'the reference database for calibration'
- vernum='FluxCalib_v1.0'
- vernum_com='version of calibration code'
-
- #keys=['cali_ref','ccdzp','ccdzpa','ccdzpb','ccdzpc','ccdzpd','ccdphoff','ccdphrms', 'phrmsA','phrmsB','phrmsC','phrmsD','aper_r','fwhm','seeing','raoff','decoff','trans','ccdnstar','nmatch','nmatcha','nmatchb','nmatchc','nmatchd', 'mdncol','colt_par','ebv','EXTNAME','cali_v']
-
- #tmpa=round(coeff0[0],4)
- #########################
- ##set header keywords: flux calibration information##
- header.set('cali_ref',cali_ref,cali_ref_com)
- header.set('COMMENT','='*66,before='cali_ref')
- header.set('COMMENT','Flux calibration information',before='cali_ref')
- header.set('COMMENT','='*66,before='cali_ref')
- header.set('ccdzp',round(coeff0[0],4),ccdzp_com)
- #header.set('ccdzpa',float('%.4f' % coeff1[0]),ccdzpa_com)
- #header.set('ccdzpb',float('%.4f' % coeff2[0]),ccdzpb_com)
- #header.set('ccdzpc',float('%.4f' % coeff3[0]),ccdzpc_com)
- #header.set('ccdzpd',float('%.4f' % coeff4[0]),ccdzpd_com)
- header.set('ccdphrms',round(coeff0[1],4),ccdphrms_com)
- #header.set('phrmsA',float('%.4f' % coeff1[1]),ccdphrms_coma)
- #header.set('phrmsB',float('%.4f' % coeff2[1]),ccdphrms_comb)
- #header.set('phrmsC',float('%.4f' % coeff3[1]),ccdphrms_comc)
- #header.set('phrmsD',float('%.4f' % coeff4[1]),ccdphrms_comd)
-
- header.set('aper_r',round(aperture_radii,4),aper_r_com)
- header.set('fwhm',float('%.4f' % fwhm),fwhm_com)
- #header.set('seeing',round(seeing,4),seeing_com)
- header.set('raoff',round(ccdraoff,4),ccdraoff_com)
- header.set('decoff',round(ccddecoff,4),ccddecoff_com)
- #header.set('trans',round(transparency,4),transparency_com)
- header.set('ccdnstar',round(ccdnstar),ccdnstar_com)
- header.set('nmatch',round(coeff0[2]),ccdnmatch_com)
- #header.set('nmatcha',round(coeff1[2]),ccdnmatcha_com)
- #header.set('nmatchb',round(coeff2[2]),ccdnmatchb_com)
- #header.set('nmatchc',round(coeff3[2]),ccdnmatchc_com)
- #header.set('nmatchd',round(coeff4[2]),ccdnmatchd_com)
- header.set('mdncol',round(gimedian,4),ccdmdncol_com)
- #header.set('colt_par',colortpar,colortpar_com)
- #header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')
-
- ###Calculate and set SKY & magnitude limiting#####################
- if not('SKYRMS' in list(header.keys())):
- #imdata=fits.getdata(image, 0)
- imdata=imgdata[1].data
- skystat=sigma_clipped_stats(imdata[500:1500,500:1500], sigma=3.)
- ccdskyrms=skystat[2]#rms/pixel of the sky in CCD
- ccdsky_com='(e-/s per pixel)'
- ccdskyrms_com='rms/pixel of the sky in unit of e-/s'
- ccdsky=skystat[1]
- header.set('sky',float('%.4f' % ccdsky),ccdsky_com)
- header.set('skyrms',float('%.4f' % ccdskyrms),ccdskyrms_com)
-
- avsky=header['SKYRMS']
- #ebv=getebv(image)
- ebv=0.0
- ccdzp=coeff0[0]
- mlim=self.getmlim(fwhm=fwhm,avsky=avsky,rdnoise=5.0,zpt=ccdzp,ebv=ebv,filter=band)
- mlim_com='magnitude limiting of 5-sigma galaxy detection'
- header.set('mlim',round(mlim,2),mlim_com)
- #######set signals of calibration progress################
- opetime=time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())
- header.set('FLUX_S',0,'flux calibration status')
- header.set('FLUX_V','1.3',vernum_com)
- header.set('FLUX_TOL',opetime,'flux calibration operation time')
-
+ vernum = 'FluxCalib_v1.0'
+ vernum_com = 'version of calibration code'
+
+ # keys=['cali_ref','ccdzp','ccdzpa','ccdzpb','ccdzpc','ccdzpd','ccdphoff','ccdphrms', 'phrmsA','phrmsB',
+ # 'phrmsC','phrmsD','aper_r','fwhm','seeing','raoff','decoff','trans','ccdnstar','nmatch','nmatcha',
+ # 'nmatchb','nmatchc','nmatchd', 'mdncol','colt_par','ebv','EXTNAME','cali_v']
+ # tmpa=round(coeff0[0],4)
+
+ # set header keywords: flux calibration information##
+ header.set('cali_ref', cali_ref, cali_ref_com)
+ header.set('COMMENT', '=' * 66, before='cali_ref')
+ header.set('COMMENT', 'Flux calibration information', before='cali_ref')
+ header.set('COMMENT', '=' * 66, before='cali_ref')
+ header.set('ccdzp', round(coeff0[0], 4), ccdzp_com)
+ # header.set('ccdzpa',float('%.4f' % coeff1[0]),ccdzpa_com)
+ # header.set('ccdzpb',float('%.4f' % coeff2[0]),ccdzpb_com)
+ # header.set('ccdzpc',float('%.4f' % coeff3[0]),ccdzpc_com)
+ # header.set('ccdzpd',float('%.4f' % coeff4[0]),ccdzpd_com)
+ header.set('ccdphrms', round(coeff0[1], 4), ccdphrms_com)
+ # header.set('phrmsA',float('%.4f' % coeff1[1]),ccdphrms_coma)
+ # header.set('phrmsB',float('%.4f' % coeff2[1]),ccdphrms_comb)
+ # header.set('phrmsC',float('%.4f' % coeff3[1]),ccdphrms_comc)
+ # header.set('phrmsD',float('%.4f' % coeff4[1]),ccdphrms_comd)
+
+ header.set('aper_r', round(aperture_radii, 4), aper_r_com)
+ header.set('fwhm', float('%.4f' % fwhm), fwhm_com)
+ # header.set('seeing',round(seeing,4),seeing_com)
+ header.set('raoff', round(ccdraoff, 4), ccdraoff_com)
+ header.set('decoff', round(ccddecoff, 4), ccddecoff_com)
+ # header.set('trans',round(transparency,4),transparency_com)
+ header.set('ccdnstar', round(ccdnstar), ccdnstar_com)
+ header.set('nmatch', round(coeff0[2]), ccdnmatch_com)
+ # header.set('nmatcha',round(coeff1[2]),ccdnmatcha_com)
+ # header.set('nmatchb',round(coeff2[2]),ccdnmatchb_com)
+ # header.set('nmatchc',round(coeff3[2]),ccdnmatchc_com)
+ # header.set('nmatchd',round(coeff4[2]),ccdnmatchd_com)
+ header.set('mdncol', round(gimedian, 4), ccdmdncol_com)
+ # header.set('colt_par',colortpar,colortpar_com)
+ # header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')
+
+ # Calculate and set SKY & magnitude limiting
+ if not ('SKYRMS' in list(header.keys())):
+ # imdata=fits.getdata(image, 0)
+ imdata = imgdata[1].data
+ skystat = sigma_clipped_stats(imdata[500:1500, 500:1500], sigma=3.)
+ ccdskyrms = skystat[2] # rms/pixel of the sky in CCD
+ ccdsky_com = '(e-/s per pixel)'
+ ccdskyrms_com = 'rms/pixel of the sky in unit of e-/s'
+ ccdsky = skystat[1]
+ header.set('sky', float('%.4f' % ccdsky), ccdsky_com)
+ header.set('skyrms', float('%.4f' % ccdskyrms), ccdskyrms_com)
+
+ avsky = header['SKYRMS']
+ # ebv=getebv(image)
+ ebv = 0.0
+ ccdzp = coeff0[0]
+ mlim = self.getmlim(fwhm=fwhm, avsky=avsky, rdnoise=5.0, zpt=ccdzp, ebv=ebv, filter=band)
+ mlim_com = 'magnitude limiting of 5-sigma galaxy detection'
+ header.set('mlim', round(mlim, 2), mlim_com)
+ # set signals of calibration progress
+ opetime = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())
+ header.set('FLUX_S', 0, 'flux calibration status')
+ header.set('FLUX_V', '1.3', vernum_com)
+ header.set('FLUX_TOL', opetime, 'flux calibration operation time')
+
if morehead:
- #catroot=cat[:cat.rfind('phot.fits')]
+ # catroot=cat[:cat.rfind('phot.fits')]
prihdu = fits.PrimaryHDU(header=header)
- prihdu.writeto(image_head_file,overwrite=True)
- #chead = fits.open(image_head_file,mode='update')
- #imh=chead[0].header
- #imh.extend(header,unique=True,update=True)
- #chead.flush()
- #chead.close()
-
- #crate and open QC1 list file to save the calibrated image names
- qc1list=os.path.join(workdir,'file_list.tmp') #the name of QC1 list
- qc1= open(qc1list,"a") #open the file
-
+ prihdu.writeto(image_head_file, overwrite=True)
+ # chead = fits.open(image_head_file,mode='update')
+ # imh=chead[0].header
+ # imh.extend(header,unique=True,update=True)
+ # chead.flush()
+ # chead.close()
+
+ # crate and open QC1 list file to save the calibrated image names
+ qc1list = os.path.join(workdir, 'file_list.tmp') # the name of QC1 list
+ qc1 = open(qc1list, "a") # open the file
+
if addhead:
- #print(('add head for ',image))
- imwht=image.replace('_img','_wht')
- imflg=image.replace('_img','_flg')
- for img,im in zip([image,imwht,imflg],[imgdata,whtdata,flgdata]):
- newimage=os.path.join(L1dir,os.path.split(img)[-1])
- newimage=newimage.replace('.fits','_L1.fits')
- #im = fits.open(img,mode='readonly')
- h0 = fits.PrimaryHDU(data=im[0].data,header=im[0].header)
- h1 = fits.ImageHDU(data=im[1].data,header=header)
- hdulist=fits.HDUList([h0,h1])
- hdulist.writeto(newimage,overwrite=True)
- #im.writeto(newimage, overwrite=True)
- #print the calibrated image name to QC1 list file
+ # print(('add head for ',image))
+ imwht = image.replace('_img', '_wht')
+ imflg = image.replace('_img', '_flg')
+ for img, im in zip([image, imwht, imflg], [imgdata, whtdata, flgdata]):
+ newimage = os.path.join(L1dir, os.path.split(img)[-1])
+ newimage = newimage.replace('.fits', '_L1.fits')
+ # im = fits.open(img,mode='readonly')
+ h0 = fits.PrimaryHDU(data=im[0].data, header=im[0].header)
+ h1 = fits.ImageHDU(data=im[1].data, header=header)
+ hdulist = fits.HDUList([h0, h1])
+ hdulist.writeto(newimage, overwrite=True)
+ # im.writeto(newimage, overwrite=True)
+ # print the calibrated image name to QC1 list file
print(newimage, file=qc1)
- #close the QC1 list file
+ # close the QC1 list file
qc1.close()
-
- #################plot######################
+
+ # plot
if plot:
print('plot calibration chart ...')
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
- fig = plt.figure(figsize = (6,5), dpi=120)
- plt.subplots_adjust(left= 0.2,bottom=0.15)
- plt.ylim(coeff0[0]-5*coeff0[1]-0.1,coeff0[0]+5*coeff0[1]+0.1)
- plt.xlim(np.min(coeff0[4])-0.2,np.max(coeff0[4])+0.2)
+ fig = plt.figure(figsize=(6, 5), dpi=120)
+ plt.subplots_adjust(left=0.2, bottom=0.15)
+ plt.ylim(coeff0[0] - 5 * coeff0[1] - 0.1, coeff0[0] + 5 * coeff0[1] + 0.1)
+ plt.xlim(np.min(coeff0[4]) - 0.2, np.max(coeff0[4]) + 0.2)
plt.xlabel(r'$\mathrm{MAG_{Gaia}}$')
plt.ylabel(r'$\mathrm{MAG_{Gaia}-MAG_{inst}}$')
- plt.plot(coeff1[4],-coeff1[3]+coeff1[4],'r.')
- plt.plot(coeff2[4],-coeff2[3]+coeff2[4],'y.')
- plt.plot(coeff3[4],-coeff3[3]+coeff3[4],'b.')
- plt.plot(coeff4[4],-coeff4[3]+coeff4[4],'g.')
- label0='zpt=%6.3f, rms:%6.3f, matched:%i' %(coeff0[0:3])
- label1='zpta=%6.3f, rms:%6.3f, matched:%i' %(coeff1[0:3])
- label2='zptb=%6.3f, rms:%6.3f, matched:%i' %(coeff2[0:3])
- label3='zptc=%6.3f, rms:%6.3f, matched:%i' %(coeff3[0:3])
- label4='zptd=%6.3f, rms:%6.3f, matched:%i' %(coeff4[0:3])
- xtmin=np.min(coeff0[4])-0.2
- xtmax=np.max(coeff0[4])+0.2
- l1,=plt.plot([xtmin,xtmax],np.array([0,0])+coeff0[0],'k-',lw=1,label=label0)
- l2,=plt.plot([xtmin,xtmax],np.array([0,0])+coeff1[0],'r-',lw=1,label=label1)
- l3,=plt.plot([xtmin,xtmax],np.array([0,0])+coeff2[0],'y-',lw=1,label=label2)
- l4,=plt.plot([xtmin,xtmax],np.array([0,0])+coeff3[0],'b-',lw=1,label=label3)
- l5,=plt.plot([xtmin,xtmax],np.array([0,0])+coeff4[0],'g-',lw=1,label=label4)
- plt.legend(handles=[l1,l2,l3,l4,l5],frameon=False,loc=2,prop={'size': 10})
+ plt.plot(coeff1[4], -coeff1[3] + coeff1[4], 'r.')
+ plt.plot(coeff2[4], -coeff2[3] + coeff2[4], 'y.')
+ plt.plot(coeff3[4], -coeff3[3] + coeff3[4], 'b.')
+ plt.plot(coeff4[4], -coeff4[3] + coeff4[4], 'g.')
+ label0 = 'zpt=%6.3f, rms:%6.3f, matched:%i' % (coeff0[0:3])
+ label1 = 'zpta=%6.3f, rms:%6.3f, matched:%i' % (coeff1[0:3])
+ label2 = 'zptb=%6.3f, rms:%6.3f, matched:%i' % (coeff2[0:3])
+ label3 = 'zptc=%6.3f, rms:%6.3f, matched:%i' % (coeff3[0:3])
+ label4 = 'zptd=%6.3f, rms:%6.3f, matched:%i' % (coeff4[0:3])
+ xtmin = np.min(coeff0[4]) - 0.2
+ xtmax = np.max(coeff0[4]) + 0.2
+ l1, = plt.plot([xtmin, xtmax], np.array([0, 0]) + coeff0[0], 'k-', lw=1, label=label0)
+ l2, = plt.plot([xtmin, xtmax], np.array([0, 0]) + coeff1[0], 'r-', lw=1, label=label1)
+ l3, = plt.plot([xtmin, xtmax], np.array([0, 0]) + coeff2[0], 'y-', lw=1, label=label2)
+ l4, = plt.plot([xtmin, xtmax], np.array([0, 0]) + coeff3[0], 'b-', lw=1, label=label3)
+ l5, = plt.plot([xtmin, xtmax], np.array([0, 0]) + coeff4[0], 'g-', lw=1, label=label4)
+ plt.legend(handles=[l1, l2, l3, l4, l5], frameon=False, loc=2, prop={'size': 10})
plt.savefig(psname)
-
- return coeff,std,match
-#---------------------------------------------------------------------------
- def run(self,fn_list,img_list=[], wht_list=[], flg_list=[],wcsdir='./',L1dir='./',workdir='./',refdir='',addhead=True,morehead=False,plot=False,nodel=True,update=False,upcat=True):
+ return coeff, std, match
+
+ def run(self, fn_list, img_list=[], wht_list=[], flg_list=[], wcsdir='./', L1dir='./', workdir='./', refdir='',
+ addhead=True, morehead=False, plot=False, nodel=True, update=False, upcat=True):
if len(fn_list) == 0:
print('Flux calibration: No input images in img_list!')
return
- #time1=time.time()
+ # time1=time.time()
print('\n\n############### run flux calibration ###############')
for i in range(len(fn_list)):
- image=workdir + fn_list[i]
+ image = workdir + fn_list[i]
if len(img_list) == len(fn_list):
- imgdata=img_list[i]
- whtdata=wht_list[i]
- flgdata=flg_list[i]
+ imgdata = img_list[i]
+ whtdata = wht_list[i]
+ flgdata = flg_list[i]
else:
print("Image data is not defined, will be read from files")
- imgdata=fits.open(image,mode='readonly')
- whtdata=fits.open(image.replace('_img','_wht'),mode='readonly')
- flgdata=fits.open(image.replace('_img','_flg'),mode='readonly')
+ imgdata = fits.open(image, mode='readonly')
+ whtdata = fits.open(image.replace('_img', '_wht'), mode='readonly')
+ flgdata = fits.open(image.replace('_img', '_flg'), mode='readonly')
- #print(('flux calibration: '+image))
+ # print(('flux calibration: '+image))
if not os.path.isfile(image):
- print(('cannot find the file:'+image))
+ print(('cannot find the file:' + image))
else:
- self.calib(image,imgdata,whtdata,flgdata,wcsdir=wcsdir,L1dir=L1dir,workdir=workdir,refdir=refdir,addhead=addhead,morehead=morehead,plot=plot,nodel=nodel,update=update)
- #time2=time.time()
+ self.calib(image, imgdata, whtdata, flgdata, wcsdir=wcsdir, L1dir=L1dir, workdir=workdir, refdir=refdir,
+ addhead=addhead, morehead=morehead, plot=plot, nodel=nodel, update=update)
+ # time2=time.time()
print('\n############### flux calibration done #############\n')
-######################################
def cleanup(self, fn_list, workdir, nodel=False):
# clean up environment
for image in fn_list:
- cat=os.path.join(workdir,image[:-5]+'.acat')
- ref=os.path.join(workdir,image[:-5]+'.rcat')
- whead=os.path.join(workdir,image[:-5]+'.whead.fits')
- #psname =os.path.join(workdir,image+'_calib.png')
+ cat = os.path.join(workdir, image[:-5] + '.acat')
+ ref = os.path.join(workdir, image[:-5] + '.rcat')
+ whead = os.path.join(workdir, image[:-5] + '.whead.fits')
+ # psname =os.path.join(workdir,image+'_calib.png')
if not nodel:
- try:
+ try:
os.remove(cat)
- except(FileNotFoundError):
+ except FileNotFoundError:
print()
try:
os.remove(ref)
- except(FileNotFoundError):
+ except FileNotFoundError:
print()
try:
os.remove(whead)
- except(FileNotFoundError):
+ except FileNotFoundError:
print()
-
--
GitLab