update flux calibration (ac59766f) · Commits · Bo Zhang / csst

csst/msc/calib_flux.py

+83 −65

Original line number	Diff line number	Diff line
		@@ -14,8 +14,9 @@ from astropy.wcs.utils import proj_plane_pixel_scales

		from .. import PACKAGE_PATH
		from ..core.processor import CsstProcessor
		from naming import CsstMscNamingRules as nr

		# Edited on Jun. 17, 2016
		# Edited on Jun. 17, 2022
		# add color term
		__author__ = 'ZZM'

		@@ -171,14 +172,15 @@ class CsstProcFluxCalibration(CsstProcessor):

		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)
		def read_inputcat(self, dm, ccdid, outcat='refcat.fits', 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)
		cat = dm.l0_cat(ccd_id)
		data = table.Table.read(cat, format='ascii')
		data.write(outcat, format='fits', overwrite=True)

		@@ -237,7 +239,7 @@ class CsstProcFluxCalibration(CsstProcessor):
		pathsplit = os.path.split(filename)
		if workdir == '':
		workdir = pathsplit[0]
		headname = os.path.join(workdir, pathsplit[1] + '.whead.fits')
		headname = os.path.join(workdir, pathsplit[1] + '_whead.fits')
		wheader = head[i].header
		prihdu = fits.PrimaryHDU(header=wheader)
		prihdu.writeto(headname, overwrite=True)
		@@ -269,22 +271,32 @@ class CsstProcFluxCalibration(CsstProcessor):

		return cat

		def prepare(self, image, wcsdir, workdir, usewcsresult=False, newcat=False):
		def prepare(self, dm, ccdid, 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')

		# 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')
		#
		workdir = dm.l1_dir
		wcscat0 = dm.pc_combined_cat
		wcscat1 = dm.l1_sci(ccd_id, suffix="img", ext="acat")
		cat = dm.l1_sci(ccd_id, suffix="img", ext="acat")
		ref = dm.l1_sci(ccd_id, suffix="img", ext="rcat")
		wcshead2 = dm.l1_sci(ccd_id, suffix="img_whead", ext="fits")

		wcshead0 = dm.pc_combined_head_fits
		wcshead0txt = dm.pc_combined_head
		wcshead1 = dm.l1_sci(ccd_id, suffix="img.acat.head", ext="fits")

		# cali_ref='GAIA' #calibration reference data

		@@ -489,21 +501,22 @@ class CsstProcFluxCalibration(CsstProcessor):
		# 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):
		def calib(self, dm, ccdid, image, imgdata, whtdata, flgdata, addhead=False,
		morehead=True, plot=False, update=False, upcat=True):
		# wcsdir='./', L1dir='./', workdir='./', refdir='',

		print('calibration for:', image)
		L1dir = dm.dir_l1()
		workdir = dm.dir_l1()
		refdir = dm.dir_l0()
		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):
		# print((image+' has been calibrated.'))
		# return
		#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')

		psname = os.path.join(workdir, imname + '_calib.png')
		if plot and os.path.isfile(psname) and (not update):
		@@ -518,7 +531,7 @@ class CsstProcFluxCalibration(CsstProcessor):

		cali_ref = 'GAIA' # calibration reference data
		# Get the photometric catalog, astrometry head
		cat, newcat, ref, header = self.prepare(image, wcsdir, workdir, newcat=False)
		cat, newcat, ref, header = self.prepare(dm, ccdid, newcat=False)

		# get obs. Information from header
		# k=list(header.keys())
		@@ -527,8 +540,8 @@ class CsstProcFluxCalibration(CsstProcessor):
		# !#airmass=header['AIRMASS']
		# exptime=header['EXPTIME']
		exptime = 150.0
		naxis1 = 9216 # header['NAXIS1']
		naxis2 = 9232 # header['NAXIS2']
		naxis1 = header['NAXIS1']
		naxis2 = 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)
		@@ -582,8 +595,7 @@ class CsstProcFluxCalibration(CsstProcessor):

		# 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
		ps = self.read_inputcat(dm, ccdid, outcat=ref)

		if np.size(ps) < 1:
		wcscheck = open(ckf, "a")
		@@ -716,10 +728,10 @@ class CsstProcFluxCalibration(CsstProcessor):
		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')]
		if morehead or (not addhead):
		image_L1_head = dm.l1_sci(ccd_id=ccdid, suffix="img_L1_head", ext=".fits")
		prihdu = fits.PrimaryHDU(header=header)
		prihdu.writeto(image_head_file, overwrite=True)
		prihdu.writeto(image_L1_head, overwrite=True)
		# chead = fits.open(image_head_file,mode='update')
		# imh=chead[0].header
		# imh.extend(header,unique=True,update=True)
		@@ -732,20 +744,22 @@ class CsstProcFluxCalibration(CsstProcessor):

		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')
		#imwht = image.replace('_img', '_wht')
		#imflg = image.replace('_img', '_flg')
		image_L1 = dm.l1_sci(ccd_id=ccdid, suffix="img_L1", ext=".fits")
		imwht_L1 = dm.l1_sci(ccd_id=ccdid, suffix="wht_L1", ext=".fits")
		imflg_L1 = dm.l1_sci(ccd_id=ccdid, suffix="flg_L1", ext=".fits")

		for img, im in zip([image_L1, imwht_L1, imflg_L1], [imgdata, whtdata, flgdata]):
		#newimage = os.path.join(L1dir, os.path.split(img)[-1])
		#newimage = newimage.replace('.fits', '_L1.fits')
		newimage = img
		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

		qc1.close()

		# plot
		@@ -781,42 +795,46 @@ class CsstProcFluxCalibration(CsstProcessor):

		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):
		def run(self, dm, img_list=[], wht_list=[], flg_list=[], wcsdir='./', L1dir='./', workdir='./', refdir='',
		addhead=True, plot=False, update=False, upcat=True):

		if len(fn_list) == 0:
		ccdids = dm.available_ccd_ids
		if len(ccdids) == 0:
		print('Flux calibration: No input images in img_list!')
		return
		# time1=time.time()
		print('\n\n############### run flux calibration ###############')
		for i in range(len(fn_list)):
		image = workdir + fn_list[i]
		if len(img_list) == len(fn_list):
		for i,ccdid in enumerate(ccdids):
		image = dm.l1_sci(ccd_id=ccdid, suffix="img", ext=".fits")
		imwht = dm.l1_sci(ccd_id=ccdid, suffix="wht", ext=".fits")
		imflg = dm.l1_sci(ccd_id=ccdid, suffix="flg", ext=".fits")

		if len(img_list) == len(ccdids):
		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')
		whtdata = fits.open(imwht, mode='readonly')
		flgdata = fits.open(imflg, mode='readonly')

		# print(('flux calibration: '+image))
		if not os.path.isfile(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)
		self.calib(dm,ccdid,image, imgdata, whtdata, flgdata, addhead=addhead, morehead=morehead, plot=plot, update=update)
		# wcsdir=wcsdir, L1dir=L1dir, workdir=workdir, refdir=refdir,
		# time2=time.time()
		print('\n############### flux calibration done #############\n')

		def cleanup(self, fn_list, workdir, nodel=False):
		def cleanup(self, dm, workdir, nodel=False):
		ccdids = dm.available_ccd_ids
		# 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')
		for ccdid in ccdids:
		cat = dm.l1_sci(ccd_id=ccdid, suffix="img", ext=".acat")
		ref = dm.l1_sci(ccd_id=ccdid, suffix="img", ext=".rcat")
		whead = dm.l1_sci(ccd_id=ccdid, suffix="img_whead", ext=".fits")
		if not nodel:
		try:
		os.remove(cat)