reformated flux calibration module (1f8e6306) · Commits · Bo Zhang / csst

csst/msc/calib_flux.py

+631 −641

Original line number	Original line	Diff line number	Diff line
	#!/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 os
			import time
			from subprocess import Popen

	import shutil
	import numpy as np		import numpy as np
			from astropy import table
	from astropy import units as u		from astropy import units as u
	from astropy.coordinates import SkyCoord		from astropy.coordinates import SkyCoord
			from astropy.io import fits
	from astropy.stats import sigma_clip		from astropy.stats import sigma_clip
	from astropy.stats import sigma_clipped_stats		from astropy.stats import sigma_clipped_stats
	from astropy.io import fits
	from astropy.wcs import WCS		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		from astropy.wcs.utils import proj_plane_pixel_scales
	import time

	from .. import PACKAGE_PATH		from .. import PACKAGE_PATH
	from ..core.processor import CsstProcessor		from ..core.processor import CsstProcessor
	#---------------------------------------------------------------------------
			# Edited on Jun. 17, 2016
			# add color term
			__author__ = 'ZZM'


	class CsstProcFluxCalibration(CsstProcessor):		class CsstProcFluxCalibration(CsstProcessor):

	def __init__(self):		def __init__(self, **kwargs):
	pass		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		# get ps1 color-corrected magnitude for specific filter
	@@ -71,7 +66,6 @@ class CsstProcFluxCalibration(CsstProcessor):

	return mags, magerr, ra, dec, brmedian		return mags, magerr, ra, dec, brmedian

	#----------------------------------------------------------------------
	def read_ps1cat(self, ra, dec, outcat, path='/line12/Pan-STARRS/chunks-qz-star-v2/', silent=True):		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		# Read a piece of the PS1 calibration catalog containing ra,dec
	from healpy import ang2pix		from healpy import ang2pix
	@@ -99,7 +93,7 @@ class CsstProcFluxCalibration(CsstProcessor):
	ps = table.vstack(ps)		ps = table.vstack(ps)
	if outcat: ps.write(outcat, format='fits', overwrite=True)		if outcat: ps.write(outcat, format='fits', overwrite=True)
	return ps		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		# function: Read a piece of the GAIA calibration catalog containing ra,dec
	from healpy import ang2pix		from healpy import ang2pix
	@@ -114,7 +108,15 @@ class CsstProcFluxCalibration(CsstProcessor):
	if npix > 8:		if npix > 8:
	print('too many healpix files:', npix, ' check image wcs!')		print('too many healpix files:', npix, ' check image wcs!')
	return		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', '?')])		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)		ps = table.Table(dtype=dt)
	for i in pix:		for i in pix:
	fname = path + 'chunk-' + '%5.5d' % i + '.fits'		fname = path + 'chunk-' + '%5.5d' % i + '.fits'
	@@ -125,7 +127,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	ps = table.vstack(ps)		ps = table.vstack(ps)
	if outcat: ps.write(outcat, format='fits', overwrite=True)		if outcat: ps.write(outcat, format='fits', overwrite=True)
	return ps		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		# get Gaia color-corrected magnitude for specific filter
	@@ -169,7 +170,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	refc = SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='fk5')		refc = SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='fk5')

	return mags, magerr, gimedian, refc		return mags, magerr, gimedian, refc
	#########################################

	def read_inputcat(self, image, outcat='refcat.fits', refdir='', silent=True):		def read_inputcat(self, image, outcat='refcat.fits', refdir='', silent=True):
	# imname=os.path.split(image)[-1]		# imname=os.path.split(image)[-1]
	@@ -183,7 +183,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	data.write(outcat, format='fits', overwrite=True)		data.write(outcat, format='fits', overwrite=True)

	return data		return data
	###########################

	def input_mags(self, cat, usepixcood=False):		def input_mags(self, cat, usepixcood=False):
	goodid = (cat['mag'] < 20) & (cat['mag'] > 16)		goodid = (cat['mag'] < 20) & (cat['mag'] > 16)
	@@ -207,7 +206,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	brmedian = 0.0		brmedian = 0.0

	return mags, magerr, brmedian, refc		return mags, magerr, brmedian, refc
	#########################################

	def rewrite_sex_cat(self, cat, workdir=''):		def rewrite_sex_cat(self, cat, workdir=''):
	data = fits.open(cat)		data = fits.open(cat)
	@@ -226,7 +224,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	hdu.writeto(subfile, overwrite=True)		hdu.writeto(subfile, overwrite=True)
	index.append(filename)		index.append(filename)
	return index		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:		if len(im_index) == 0:
	@@ -244,7 +241,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	wheader = head[i].header		wheader = head[i].header
	prihdu = fits.PrimaryHDU(header=wheader)		prihdu = fits.PrimaryHDU(header=wheader)
	prihdu.writeto(headname, overwrite=True)		prihdu.writeto(headname, overwrite=True)
	##########################

	def run_sextractor(self, image, cat):		def run_sextractor(self, image, cat):
	image0 = image		image0 = image
	@@ -272,7 +268,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	p.wait()		p.wait()

	return cat		return cat
	##########################

	def prepare(self, image, wcsdir, workdir, usewcsresult=False, newcat=False):		def prepare(self, image, wcsdir, workdir, usewcsresult=False, newcat=False):
	# make calibration files: Obs. cat (SExtractor), PSF model (PSFex)		# make calibration files: Obs. cat (SExtractor), PSF model (PSFex)
	@@ -310,7 +305,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	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		# ebv=0.0
	@@ -325,7 +319,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	# if tid.sum()==1: ebv=tiles['EBV'][tid].data[0]		# if tid.sum()==1: ebv=tiles['EBV'][tid].data[0]
	# #header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')		# #header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')
	# return ebv		# 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		# E(B-V): 3.995, 3.214, 2.165, 1.592, 1.211, 1.064 for ugrizY decals
	@@ -342,7 +335,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	flim = noise * snr		flim = noise * snr
	mlim = -2.5 * np.log10(flim) + zpt - k * ebv		mlim = -2.5 * np.log10(flim) + zpt - k * ebv
	return mlim		return mlim
	#########################

	def getfwhm(self, fwhmsex, ellip, obsme, obsflags):		def getfwhm(self, fwhmsex, ellip, obsme, obsflags):
	nanid = np.isnan(fwhmsex)		nanid = np.isnan(fwhmsex)
	@@ -368,7 +360,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	# print id,fwhmid.data[~fwhmid.mask][id]		# print id,fwhmid.data[~fwhmid.mask][id]
	# print fwhmid, fwhmsex		# print fwhmid, fwhmsex
	return fwhm		return fwhm
	#########################

	def makedatedir(self, path):		def makedatedir(self, path):
	index = path.rfind('/')		index = path.rfind('/')
	@@ -403,7 +394,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	f1.close()		f1.close()
	f.close()		f.close()
	return wcshead		return wcshead
	##########################

	def combine_head(self, image, wcshead1='', wcshead2='', prime=False):		def combine_head(self, image, wcshead1='', wcshead2='', prime=False):
	# combine image head and wcs head keywords		# combine image head and wcs head keywords
	@@ -433,7 +423,6 @@ class CsstProcFluxCalibration(CsstProcessor):
	h0['CTYPE2'] = 'DEC--TPV'		h0['CTYPE2'] = 'DEC--TPV'

	return h0		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		# matching obs and ref catalogs
	@@ -468,7 +457,8 @@ class CsstProcFluxCalibration(CsstProcessor):
	obs_uidlim = obs_uid[uidlim]		obs_uidlim = obs_uid[uidlim]
	ref_uidlim = ref_uid[uidlim]		ref_uidlim = ref_uid[uidlim]

	ccdraoff=np.median((obsc[obs_uidlim].ra- refc[ref_uidlim].ra).arcsecnp.cos(obsc[obs_uidlim].dec.degnp.pi/180))		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)		ccddecoff = np.median((obsc[obs_uidlim].dec - refc[ref_uidlim].dec).arcsec)
	# if fwhmsex.size >1:		# if fwhmsex.size >1:
	# fwhm=np.median(fwhmsex[obs_uidlim])		# fwhm=np.median(fwhmsex[obs_uidlim])
	@@ -498,9 +488,9 @@ class CsstProcFluxCalibration(CsstProcessor):
	# print(coeff0,cstd,csize,cobsm,crefm,ccdraoff,ccddecoff)		# print(coeff0,cstd,csize,cobsm,crefm,ccdraoff,ccddecoff)
	# print('crefm.min,max=',refm.min(),refm.max())		# print('crefm.min,max=',refm.min(),refm.max())
	return coeff0, cstd, csize, cobsm, crefm, ccdraoff, ccddecoff, fwhm		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, 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('calibration for:', image)
	if not os.path.exists(L1dir):		if not os.path.exists(L1dir):
	@@ -664,11 +654,12 @@ class CsstProcFluxCalibration(CsstProcessor):
	vernum = 'FluxCalib_v1.0'		vernum = 'FluxCalib_v1.0'
	vernum_com = 'version of calibration code'		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']		# 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)		# tmpa=round(coeff0[0],4)
	#########################
	##set header keywords: flux calibration information##		# set header keywords: flux calibration information##
	header.set('cali_ref', cali_ref, cali_ref_com)		header.set('cali_ref', cali_ref, cali_ref_com)
	header.set('COMMENT', '=' * 66, before='cali_ref')		header.set('COMMENT', '=' * 66, before='cali_ref')
	header.set('COMMENT', 'Flux calibration information', before='cali_ref')		header.set('COMMENT', 'Flux calibration information', before='cali_ref')
	@@ -700,7 +691,7 @@ class CsstProcFluxCalibration(CsstProcessor):
	# header.set('colt_par',colortpar,colortpar_com)		# header.set('colt_par',colortpar,colortpar_com)
	# header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')		# header.set('ebv',round(ebv,4),'E(B-V) from SFD1998')

	###Calculate and set SKY & magnitude limiting#####################		# Calculate and set SKY & magnitude limiting
	if not ('SKYRMS' in list(header.keys())):		if not ('SKYRMS' in list(header.keys())):
	# imdata=fits.getdata(image, 0)		# imdata=fits.getdata(image, 0)
	imdata = imgdata[1].data		imdata = imgdata[1].data
	@@ -719,7 +710,7 @@ class CsstProcFluxCalibration(CsstProcessor):
	mlim = self.getmlim(fwhm=fwhm, avsky=avsky, rdnoise=5.0, zpt=ccdzp, ebv=ebv, filter=band)		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'		mlim_com = 'magnitude limiting of 5-sigma galaxy detection'
	header.set('mlim', round(mlim, 2), mlim_com)		header.set('mlim', round(mlim, 2), mlim_com)
	#######set signals of calibration progress################		# set signals of calibration progress
	opetime = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())		opetime = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())
	header.set('FLUX_S', 0, 'flux calibration status')		header.set('FLUX_S', 0, 'flux calibration status')
	header.set('FLUX_V', '1.3', vernum_com)		header.set('FLUX_V', '1.3', vernum_com)
	@@ -757,7 +748,7 @@ class CsstProcFluxCalibration(CsstProcessor):
	# close the QC1 list file		# close the QC1 list file
	qc1.close()		qc1.close()

	#################plot######################		# plot
	if plot:		if plot:
	print('plot calibration chart ...')		print('plot calibration chart ...')
	import matplotlib		import matplotlib
	@@ -789,9 +780,9 @@ class CsstProcFluxCalibration(CsstProcessor):
	plt.savefig(psname)		plt.savefig(psname)

	return coeff, std, match		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, 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:		if len(fn_list) == 0:
	print('Flux calibration: No input images in img_list!')		print('Flux calibration: No input images in img_list!')
	@@ -814,11 +805,11 @@ class CsstProcFluxCalibration(CsstProcessor):
	if not os.path.isfile(image):		if not os.path.isfile(image):
	print(('cannot find the file:' + image))		print(('cannot find the file:' + image))
	else:		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(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()		# time2=time.time()
	print('\n############### flux calibration done #############\n')		print('\n############### flux calibration done #############\n')

	######################################
	def cleanup(self, fn_list, workdir, nodel=False):		def cleanup(self, fn_list, workdir, nodel=False):
	# clean up environment		# clean up environment
	for image in fn_list:		for image in fn_list:
	@@ -829,14 +820,13 @@ class CsstProcFluxCalibration(CsstProcessor):
	if not nodel:		if not nodel:
	try:		try:
	os.remove(cat)		os.remove(cat)
	except(FileNotFoundError):		except FileNotFoundError:
	print()		print()
	try:		try:
	os.remove(ref)		os.remove(ref)
	except(FileNotFoundError):		except FileNotFoundError:
	print()		print()
	try:		try:
	os.remove(whead)		os.remove(whead)
	except(FileNotFoundError):		except FileNotFoundError:
	print()		print()