#i3meds.py import meds import os import numpy as np from . import utils from . import structs DEFAULT_OPTIONS = dict( starting_params_only=False, rescale_weight=False, first_exposure=1, max_nexposures=None, use_nearest_pix_mask=False, reject_outliers=True, replace_psf_path=None, psfex_rerun_version=None, psf_input_method='psfex', ) class I3MultiExposureInputs(object): def __init__(self): self.exposures = [] self.meta = {} def __getitem__(self, n): return self.exposures[n] def __len__(self): return len(self.exposures) def __iter__(self): return iter(self.exposures) def all(self, name): return [info[name] for info in self.exposures] def append(self, exposure): self.exposures.append(exposure) class I3MEDS(meds.MEDS): def __init__(self, *args, **kwargs): super(I3MEDS, self).__init__(*args, **kwargs) self.saved_filename = args[0] #I found that loading the PSFs from file #is taking about 15-20% of the run time on #edison. So I'm caching them self.psfex_cache={} #This is only relevant for GREAT-DES #In the Great-DES scheme we just look up the index #of the psf to use in the truth table self.truth_filename = self.saved_filename.replace(".meds.", ".truth.").replace(".noisefree","") self.psf_id_col = None self.tilename = utils.find_tilename(self.saved_filename) def get_psfex_psf(self, iobj, iexp, options): desdata_old=self.get_meta()['DESDATA'][0] desdata_new=os.environ['DESDATA'] psf_path = self.get_source_path(iobj, iexp).replace( '.fits.fz','_psfcat.psf').replace( desdata_old,desdata_new) # fiddle with the PSF path, in one of two ways, # a replacement or a specific new dir if options.replace_psf_path: psf_path = psf_path.replace(options.replace_psf_path.split(',', maxsplit=1)) #specific new dir for re-run psfs if options.psfex_rerun_version: psf_path = psf_path.replace("OPS", "EXTRA") psf_path = psf_path.replace("red/DECam", "psfex-rerun/{}/DECam".format(options.psfex_rerun_version)) #Load the PSFEx data import galsim.des try: psfex_i = galsim.des.DES_PSFEx(psf_path) except: print "PSF path missing", psf_path return None #Get the image array return self.extract_psfex_psf(psfex_i, iobj, iexp, options) def get_bundled_psfex_psf(self, iobj, iexp, options): #Find the exposure name source_path = self.get_source_path(iobj, iexp) source = os.path.split(source_path)[1].split('.')[0] #get the HDU name corresponding to that name. #Might need to tweak that naming hdu_name = "psf_"+options.psfex_rerun_version+"_"+source+"_psfcat" #Maybe we have it in the cache already if hdu_name in self.psfex_cache: psfex_i = self.psfex_cache[hdu_name] else: #Turn the HDU into a PSFEx object #PSFEx expects a pyfits HDU not fitsio. #This is insane. I know this. import galsim.des hdu = galsim.pyfits.open(self._filename)[hdu_name] psfex_i = galsim.des.DES_PSFEx(hdu) self.psfex_cache[hdu_name] = psfex_i #Get the image array return self.extract_psfex_psf(psfex_i, iobj, iexp, options) def extract_psfex_psf(self, psfex_i, iobj, iexp, options): psf_size=(options.stamp_size+options.padding)*options.upsampling # locations in PSF data at which to interpolate orig_col = self['orig_col'][iobj][iexp] orig_row = self['orig_row'][iobj][iexp] psf = utils.getPSFExarray(psfex_i, orig_col, orig_row, psf_size, psf_size, options.upsampling) return psf def get_fake_psf(self, iobj, iexp, options): logger.warning("Warning: Using false PSF") psf = utils.small_round_psf_image(psf_size) return psf def get_great_des_psf(self, iobj, iexp, options): import pyfits #only load the table once per process if self.psf_id_col is None: self.psf_id_col = pyfits.getdata(self.truth_filename)['id_psf'] #find and load the corresponding PSF image, #again from a fixed filename psf_id = self.psf_id_col[iobj] psf_filename = os.path.split(self.truth_filename)[0]+ "/nbc.psf.hires.nsub%d.fits"%options.upsampling if not os.path.exists(psf_filename): raise ValueError("Could not find PSF file for specified upsampling: %s"%psf_filename) #load psf, just once since single exposure psf = pyfits.getdata(psf_filename, psf_id) return psf def get_psf(self, iobj, iexp, options, psf_input_method): """ @brief get a list of numpy arrays of PSF images for each exposure @options im3shape options struct @use_psfex if yes, use the PSFEx files, otherwise use dummy PSF @iobj object id in the meds file @selection list of indices of cutouts to use @return list of numpy arrays with images of PSFs """ # get the exposure images list, if there is a coadd (n_exposures>1) then skip it # n_cutout - number of available cutouts # n_image - number of single exposure images, or 1 if n_cutout==1 # index_first_image - 0 if n_cutout==1 and 1 if n_cutout!=1 psf_size=(options.stamp_size+options.padding)*options.upsampling # decide whath PSF to use if psf_input_method is None: psf_input_method='none' psf_input_method = psf_input_method.lower() psf_function = { 'psfex': self.get_psfex_psf, 'none': self.get_fake_psf, 'bundled-psfex': self.get_bundled_psfex_psf, 'great_des': self.get_great_des_psf, }.get(psf_input_method) if psf_function is None: raise ValueError('"%s" unknown PSF input method' % psf_input_method) psf = psf_function(iobj, iexp, options) return psf def get_i3transform(self, iobj, exposure_id, stamp_size): # Build the transform from sky to image coordinates transform = structs.Transform_struct() transform.ra0 = 0.0 # these are not used in this case transform.dec0 = 0.0 transform.cosdec0 = 0.0 transform.x0 = self['cutout_col'][iobj][exposure_id]+0.5 transform.y0 = self['cutout_row'][iobj][exposure_id]+0.5 dudcol=self['dudcol'][iobj][exposure_id] dudrow=self['dudrow'][iobj][exposure_id] dvdcol=self['dvdcol'][iobj][exposure_id] dvdrow=self['dvdrow'][iobj][exposure_id] B = np.array([ [dudcol,dudrow],[dvdcol,dvdrow]]) A = np.linalg.inv(B) transform.A[0][0] = A[0,0] transform.A[0][1] = A[0,1] transform.A[1][0] = A[1,0] transform.A[1][1] = A[1,1] #print 'transform det = ', np.linalg.det(B) #print 'A = ', A return transform def i3_save_source_paths(self, filename): """ Save the list of source paths into a catalog with indices""" outfile = open(filename, 'w') for i,path in enumerate(self.i3_source_path_list): outfile.write("%d %s\n"% (i,path.strip())) outfile.close() def get_i3transform_inverse_matrix(self, iobj, exposure_id, stamp_size): transform = self.get_i3transform(iobj, exposure_id, stamp_size) A = transform.A A = [[A[0][0],A[0][1]],[A[1][0],A[1][1]]] return np.linalg.inv(A) # This next function is taken from wcsutils.py in the ucl_des_shear repository def convert_g_image2sky(self, iobj, exposure_id, stamp_size, g1image, g2image): """Return the ellipticity (g1sky, g2sky) in sky coordinates corresponding to the input (g1image, g2image). Uses the ellipticity convention |g| = (a-b)/(a+b). Currently only works for scalar input g1image, g2image, but can be called in list comprehensions. """ Aimage2sky = self.get_i3transform_inverse_matrix(iobj, exposure_id, stamp_size) e1image, e2image = utils.convert_e_linear_to_quadratic(g1image, g2image) # Build the ellipticity matrix Ei = np.array(((1. + e1image, e2image), (e2image, 1. - e1image))) # Perform the transformation Es = np.dot(np.dot(Aimage2sky, Ei), Aimage2sky.T) # Extract the ellipticities and convert back to the linear |g|=(a-b)/(a+b) ellips Estrace = Es.trace() e1sky, e2sky = (Es[0, 0] - Es[1, 1]) / Estrace, 2. * Es[0, 1] / Estrace g1sky, g2sky = utils.convert_e_quadratic_to_linear(e1sky, e2sky) return g1sky, g2sky def convert_g_sky2image(self, iobj, exposure_id, stamp_size, g1sky, g2sky): """Return the ellipticity (g1image, g2image) in image coordinates corresponding to the input (g1sky, g2sky). Uses the ellipticity convention |g| = (a-b)/(a+b). Currently only works for scalar input g1sky, g2sky, but can be called in list comprehensions. """ transform = self.get_i3transform(iobj, exposure_id, stamp_size) Asky2image = np.array(transform.A) e1sky, e2sky = utils.convert_e_linear_to_quadratic(g1sky, g2sky) # Build the ellipticity matrix Es = np.array(((1. + e1sky, e2sky), (e2sky, 1. - e1sky))) # Perform the transformation Ei = np.dot(np.dot(Asky2image, Es), Asky2image.T) # Extract the ellipticities and convert back to the linear |g|=(a-b)/(a+b) ellips Eitrace = Ei.trace() e1image, e2image = (Ei[0, 0] - Ei[1, 1]) / Eitrace, 2. * Ei[0, 1] / Eitrace g1image, g2image = utils.convert_e_quadratic_to_linear(e1image, e2image) return g1image, g2image def select_exposures(self, iobj, options, first, last): #Use all exposures if not using flags if not options.use_image_flags: return range(first,last) #Return a list of indices of exposures to select. #Get the image flag for each exposure flags = [ self.get_source_info(iobj, i)['image_flags'] for i in xrange(first, last) ] #Return a list of indices where the image flag is zero selected_exposures = [i+1 for (i,f) in enumerate(flags) if f==0] return selected_exposures def get_ccd_info(self, iobj, iexp): try: source_path = self.get_source_path(iobj, iexp) except: return 0, "N/A" source_base = os.path.basename(source_path).split('.')[0] decam, exposure_name, ccd = source_base.split("_") ccd = int(ccd) return ccd, exposure_name def run_uberseg(self, iobj, iexp, weight): ID = self.get_number(iobj) seg = self.interpolate_coadd_seg(iobj,iexp) uber_mask = utils.uberseg(ID,weight,seg) return uber_mask def get_im3shape_weight(self, iobj, iexp, image, options): weight = self.get_cweight_cutout(iobj, iexp) #blank out the weight mask where other objects #are assigned if options.uberseg: uber_mask=self.run_uberseg(iobj,iexp,weight) else: uber_mask=None #Rescale weights (i.e. check they are correctly normalized) if options.rescale_weight: utils.rescale_weight_edge(image, weight, options) # set the negative weights to zero (those are flags) weight_max = weight.max() weight[weightpixel transform stamp_size = image.shape[0] transform = self.get_i3transform(iobj, iexp, stamp_size) #Auxiliary diagnostic data ccd, exposure_name = self.get_ccd_info(iobj, iexp) orig_row = self['orig_row'][iobj][iexp] orig_col = self['orig_col'][iobj][iexp] info = { 'image' : image, 'weight' : weight, 'psf' : psf, 'uber_mask' : uber_mask, 'stamp_size' : stamp_size, 'orig_row' : orig_row, 'orig_col' : orig_col, 'exposure_name':exposure_name, 'ccd': ccd, 'transform':transform, 'iexp':iexp, 'iobj':iobj, } return info def get_im3shape_inputs(self, iobj, options): """ @brief analyse a multi-exposure galaxy @param options Im3shape options structs @param psf_input_method choose from {psfex,single,none} @param iobj id of the object in the MEDS file """ #Get the number of exposures from the MEDS file n_image = self.get_cat()['ncutout'][iobj] ID = self.get_number(iobj) coadd_objects_id = self['id'][iobj] if n_image==1: first=0 last=1 else: first=1 last=n_image #Choose which exposures we want and #get all the info for them selection = self.select_exposures(iobj, options, first, last) exposures = [self.get_exposure_info(iobj, iexp, options) for iexp in selection] #Get a few bits and bobs of metadata inputs = I3MultiExposureInputs() inputs.meta['ID'] = coadd_objects_id inputs.meta['iobj'] = iobj inputs.meta['row_id'] = ID inputs.meta['tilename'] = self.tilename ra, dec = self.get_radec(iobj) if ra is not None: inputs.meta['ra'] = ra inputs.meta['dec'] = dec #Choose which exposures to use and save them selection_bitmask = 0 for e,exposure in zip(selection,exposures): if exposure is not None: selection_bitmask |= (1<