import py3shape import meds from inspect import currentframe, getframeinfo from .analyze import get_psf_params, compute_rgpp_rp, convert_psf_ellipticities_to_wcs from .common import MAX_N_EXPOSURE from . import utils import sys import os import numpy as np import argparse import logging import py3shape.output import time FLAG_ERROR = 1 FLAG_OK = 0 FLAG_INCOMPLETE = 2 DES_PIXEL_SCALE = 0.27 logging.basicConfig(format="%(message)s", level=logging.INFO, stream=sys.stdout) logger = logging.getLogger("des_meds") psf_id_col = None def saveResidualImages(id,image_data): import pyfits import pylab import numpy image = numpy.concatenate(image_data[0]) weight = numpy.concatenate(image_data[1]) psf = numpy.concatenate(image_data[2]) bestfit = image_data[3] # this assumes that im3shape got the scale right, anyway doesn't matter for shape scale1 = (bestfit * weight).sum() scale2 = (image * weight).sum() residual = bestfit/scale1 - image/scale2 masked = (image * weight) if len(image_data[0]) == 1: nx = 1 ny = 6 else: nx = 6 ny = 1 pylab.figure() pylab.subplot(nx,ny,1) pylab.imshow(image,interpolation='none') pylab.xticks([len(image) -1],['%d' % len(image)]) pylab.yticks([len(image) -1],['%d' % len(image)]) pylab.title('orig') pylab.subplot(nx,ny,2) pylab.imshow(bestfit,interpolation='none') pylab.xticks([]) pylab.yticks([]) pylab.title('im3shape') pylab.subplot(nx,ny,3) pylab.imshow(masked,interpolation='none') pylab.xticks([]) pylab.yticks([]) pylab.title('masked') pylab.subplot(nx,ny,4) pylab.imshow(residual,interpolation='none') pylab.xticks([]) pylab.yticks([]) pylab.title('residual') pylab.subplot(nx,ny,5) pylab.imshow(weight,interpolation='none') pylab.xticks([]) pylab.yticks([]) pylab.title('weight') pylab.subplot(nx,ny,6) pylab.imshow(psf,interpolation='none') pylab.xticks([]) pylab.yticks([]) pylab.title('PSF') filename_fig = 'images.%d.png' % id pylab.savefig(filename_fig,bbox_inches='tight',dpi=100) pylab.close() logger.info('saved %s' % filename_fig) #Stolen from Barney Rowe. #For now I'm not going to use this but instead just use a small round gaussian def getPSFExarray(psfex, x_image, y_image, nsidex=32, nsidey=32, upsampling=1, offset=None): """Return an image of the PSFEx model of the PSF as a NumPy array. Arguments --------- psfex A galsim.des.PSFEx instance opened using, for example, `psfex = galsim.des.DES_PSFEx(psfex_file_name)`. x_image Floating point x position on image [pixels] y_image Floating point y position on image [pixels] nsidex Size of PSF image along x [pixels] nsidey Size of PSF image along y [pixels] upsampling Upsampling (see Zuntz et al 2013) Returns a NumPy array with shape (nsidey, nsidex) - note the reversal of y and x to match the NumPy internal [y, x] style array ordering. This is to ensure that `pyfits.writeto()` using the ouput array creates FITS-compliant output. """ import galsim image = galsim.ImageD(nsidex, nsidey) psf = psfex.getPSF(galsim.PositionD(x_image, y_image)) psf.draw(image, scale=1.0/upsampling, offset=offset) return image.array def getShapeletPSFarray(des_shapelet, x_image, y_image, nsidex=32, nsidey=32, upsampling=1, offset=None): """Return an image of the PSFEx model of the PSF as a NumPy array. Arguments --------- des_shapelet A galsim.des.DES_Shapelet instance opened using, for example, `des_shapelet = galsim.des.DES_Shapelet(shapelet_file_name)`. Usually stored as *_fitpsf.fits files. x_image Floating point x position on image [pixels] y_image Floating point y position on image [pixels] nsidex Size of PSF image along x [pixels] nsidey Size of PSF image along y [pixels] upsampling Upsampling (see Zuntz et al 2013) Returns a NumPy array with shape (nsidey, nsidex) - note the reversal of y and x to match the NumPy internal [y, x] style array ordering. This is to ensure that `pyfits.writeto()` using the ouput array creates FITS-compliant output. """ import galsim image = galsim.ImageD(nsidex, nsidey) psf = des_shapelet.getPSF(galsim.PositionD(x_image, y_image), pixel_scale=1.) psf.draw(image, scale=1. / upsampling, offset=offset) return image.array def uberseg(object_id,weight,seg): #Object id in seg map should be #First check that expected object is in seg map if object_id not in seg: print 'ID not in seg...' raise ValueError #First get all indices of all seg map pixels which contain an object i.e. are not equal to zero obj_inds = np.where(seg!=0) #Then loop through pixels in seg map, check which obj ind it is closest to. #If the closest obj ind does not correspond to the target, set this pixel in the weight map to zero. uber_weight = np.ones_like(weight) for i,row in enumerate(seg): for j, element in enumerate(row): obj_dists = (i-obj_inds[0])**2 + (j-obj_inds[1])**2 ind_min=np.argmin(obj_dists) if seg[obj_inds[0][ind_min],obj_inds[1][ind_min]] != object_id: uber_weight[i,j] = 0. weight*=uber_weight return uber_weight class I3MEDS(meds.MEDS): def __init__(self, *args, **kwargs): self.replace_psf_path=kwargs.pop("replace_psf_path", None) self.psfex_rerun_version=kwargs.pop("psfex_rerun_version", None) super(I3MEDS, self).__init__(*args, **kwargs) self.i3_source_path_list = self.get_image_info()['image_path'].tolist() self.saved_filename = args[0] self.mask_dilation = 0 def get_PSFs(self, options, psf_input_method, iobj, max_nexposures=None): """ @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 @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 n_cutouts = self.get_cat()['ncutout'][iobj] if n_cutouts == 1: n_image = 1 index_first_image = 0 else: n_image = n_cutouts-1 index_first_image = 1 if max_nexposures is not None: n_image = min(n_image, max_nexposures) # get the number of exposures, skip the coadd # get the PSF size psf_size=(options.stamp_size+options.padding)*options.upsampling # decide whath PSF to use if psf_input_method.lower() == 'psfex': desdata_old=self.get_meta()['DESDATA'][0] desdata_new=os.environ['DESDATA'] import galsim.des # find the PSFEx file using filenames psf_paths = [ self.get_source_path(iobj, i).replace('.fits.fz','_psfcat.psf').replace(desdata_old,desdata_new) for i in xrange(index_first_image,n_image+index_first_image) ] if self.replace_psf_path is not None: psf_paths = [p.replace(*self.replace_psf_path) for p in psf_paths] #.replace('/astro/astronfs03/workarea/mjarvis/','/global/project/projectdirs/des/wl/desdata/wlpipe/end-to-end/') if self.psfex_rerun_version: new_psf_paths = [] for p in psf_paths: p = p.replace("OPS", "EXTRA") p = p.replace("red/DECam", "psfex-rerun/{}/DECam".format(self.psfex_rerun_version)) new_psf_paths.append(p) psf_paths = new_psf_paths # create PSFEx images orig_cols = self['orig_col'][iobj][index_first_image:] orig_rows = self['orig_row'][iobj][index_first_image:] #Loop through PSFs psfs = [] for i in xrange(n_image): try: psfex_i = galsim.des.DES_PSFEx(psf_paths[i]) psf = getPSFExarray(psfex_i, orig_cols[i], orig_rows[i], psf_size, psf_size, options.upsampling) except: psf = None psfs.append(psf) elif psf_input_method.lower() == 'shapelet': try: shapelet_path = os.environ['SHAPELET_PATH'] except: logger.error("Error: Set environment variable SHAPELET_PATH before running im3shape. Otherwise don't know there to look for shaplet PSF files.") # Construct shapelet PSF filenames psf_paths = [] for i in xrange(index_first_image,n_image+index_first_image): exposure, ccd = m.get_source_info(10,3)[0].split('/')[-1].split('.')[0].split('_')[1:] shapelet_filename = '%s/DECam_%s/DECam_%s_%s_psf.fits' % \ (shapelet_path, str(exposure), str(exposure), str(ccd)) psf_paths.append(shapelet_filename) # create Shapelets objects shapelets = [ galsim.des.DES_Shapelet(psf_path) for psf_path in psf_paths ] # create Shapelets images orig_cols = self['orig_col'][iobj][index_first_image:] orig_rows = self['orig_row'][iobj][index_first_image:] psfs = [ getShapeletPSFarray(shapelets[i], orig_cols[i], orig_rows[i], psf_size, psf_size, options.upsampling) for i in xrange(n_image) ] elif psf_input_method == None or psf_input_method.lower()=='none': logger.warning("Warning: Using false PSF") psfs = [py3shape.utils.small_round_psf_image(psf_size) for i in xrange(index_first_image,n_image+index_first_image)] elif psf_input_method.lower() == 'single': import pyfits psf_single = pyfits.getdata(args.filename_psf_image) psfs = [psf_single for i in xrange(index_first_image,n_image+index_first_image) ] elif psf_input_method.lower() == 'single_moffat': import pyfits # psf_single = pyfits.getdata(args.filename_psf_image) psf_beta, psf_fwhm, psf_e1, psf_e2 = np.loadtxt(args.filename_psf_single_moffat) import galsim psf_image = galsim.ImageD(psf_size, psf_size) psf = galsim.Moffat(beta=psf_beta, fwhm=psf_fwhm) psf.applyShear(e1=psf_e1, e2=psf_e2) pix = galsim.Pixel(DES_PIXEL_SCALE) psfpix = galsim.Convolve([psf,pix]) psfpix.draw(psf_image, scale=DES_PIXEL_SCALE / float(options.upsampling)) psfs = [psf_image.array for i in xrange(index_first_image,n_image+index_first_image) ] elif psf_input_method.lower() == 'great_des': import pyfits #In the Great-DES scheme we just look up the index #of the psf to use in the truth table truth_filename = self.saved_filename.replace(".meds.", ".truth.").replace(".noisefree","") #only load the table once per process global psf_id_col if psf_id_col is None: psf_id_col = pyfits.getdata(truth_filename)['id_psf'] #find and load the corresponding PSF image, #again from a fixed filename psf_id = psf_id_col[iobj] logger.debug('using psf_id = %d'% psf_id) psf_filename = os.path.split(truth_filename)[0]+ "/nbc.psf.hires.fits" #load psf, just once since single exposure psfs = [pyfits.getdata(psf_filename, psf_id)] logger.debug('psf sum = %f'% psfs[0].sum()) else: logger.error('"%s" unknown PSF input method' % psf_input_method) raise NameError('"%s" unknown PSF input method' % psf_input_method) return psfs def get_i3transform(self, iobj, exposure_id, stamp_size): # Build the transform from sky to image coordinates transform = py3shape.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_object_source_path_indices(self, iobj, index_first_image, n_image): # Get the path to each of the exposures that we have just run on. # This is useful output metadata source_path_indices = [-1 for i in xrange(py3shape.structs.MAX_N_EXPOSURE)] for icutout in xrange(index_first_image,n_image+index_first_image): #find the path to a given object path = self.get_source_path(iobj, icutout) #rather than write out a whole long string just use an index into a list of the strings #because there are only e.g. 30 images and thousands of galaxies. #The max handles the fact we might be running on the coadds index = max(icutout-1, 0) source_path_indices[index] = self.i3_source_path_list.index(path) return source_path_indices 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 @staticmethod def dilate_mask(count, weight): mask = (weight!=0) ny,nx=mask.shape for iteration in xrange(count): m = np.ones_like(mask, dtype=bool) it = np.nditer(mask, flags=['multi_index']) while not it.finished: i,j = it.multi_index imin=max(i-1,0) imax=min(i+1,ny-1) jmin=max(j-1,0) jmax=min(j+1,nx-1) if (~mask[imin:imax+1,j]).any() or (~mask[i,jmin:jmax+1]).any(): m[i,j] = 0 it.iternext() mask *= m weight *= mask def select_exposures(self, iobj, psfs, first, n): #Return a list of indices of exposures to select. flags = [ self.get_source_info(iobj, i)['image_flags'] for i in xrange(first, first+n) ] #Eliminate cases where we didn't find the PSF file for i in xrange(n): if psfs[i] is None: flags[i] = 1 selected_exposures = [i for i in xrange(n) if flags[i]==0] return selected_exposures def im3shape_analyze(self, options, psf_input_method, iobj, starting_params_only=False, rescale_weight=False, first_exposure=1, max_nexposures=None, use_nearest_pix_mask=False, reject_outliers=True): """ @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 exposure images list, if there is a coadd (n_exposures>1) then skip it n_cutouts = self.get_cat()['ncutout'][iobj] ID = self.get_number(iobj) coadd_objects_id = self['id'][iobj] #Decide the first exposure. 0 in case we are just doing coadds, #1 in the normal case otherwise, and we can change it to skip some objects if n_cutouts == 1: n_image = 1 index_first_image = 0 else: n_image = n_cutouts-1 index_first_image = first_exposure #Decide the last exposure. Usually just the final one but we can cut it short if max_nexposures is None: index_last_image = None else: index_last_image = index_first_image+max_nexposures #Error check - we have a max exposures limits if n_image>py3shape.structs.MAX_N_EXPOSURE: message = "OBJECT %d (SEXTRACTOR ID %d ) HAS %d EXPOSURES IN FILE %s - MORE THAN MAX OF %d" % (iobj, ID, n_image, self._filename, py3shape.structs.MAX_N_EXPOSURE) raise ValueError(message) # Get the image data to run on images = self.get_cutout_list(iobj, type='image')[index_first_image:index_last_image] exposure_names = [] try: ccds = [] tmp_index_last_image = index_last_image if tmp_index_last_image is None: tmp_index_last_image = index_first_image+n_image for i in xrange(index_first_image, tmp_index_last_image): source_path_i = self.get_source_path(iobj, i) source_file_i = os.path.basename(source_path_i) dot_index = source_file_i.index('.') source_base_i = source_file_i[:dot_index] decam, exposure_name, ccd = source_base_i.split("_") ccd = int(ccd) exposure_names.append(decam+"_"+exposure_name) ccds.append(ccd) except: import warnings; warnings.simplefilter('once'); warnings.warn('GREAT-DES mode, using dummy ccd index') ccds.append(0) if len(images)==0: message = "OBJECT %d (SEXTRACTOR ID %d ) HAS NO EXPOSURES IN RANGE %d - %d" % (iobj, ID, index_first_image, index_last_image) raise ValueError(message) # get the exposure weights, 0-th is the coadd, so skipping try: weights = self.get_cweight_cutout_list(iobj)[index_first_image:index_last_image] #The test card images have missing weight dat aso handle this in an ugly way except IndexError: print "Warning! Temporary workaround for weights! Get rid of this!" weights = [1000*np.ones(im.shape) for im in images] if reject_outliers: nreject = meds.reject_outliers(images, weights) else: nreject=0 if use_nearest_pix_mask: #Get single epoch indixes to loop through last_se_index = index_last_image if last_se_index is None: last_se_index = len(images) SE_inds=range(index_first_image,last_se_index+1) for i,SE_ind in enumerate(SE_inds): #interpolate coadd seg to frame of each SE image SE_seg = self.interpolate_coadd_seg(iobj,SE_ind) uber_mask = uberseg(ID,weights[i],SE_seg) else: uber_mask = None #Rescale weights (i.e. check they are correctly normalized) if rescale_weight: for image,weight in zip(images,weights): utils.rescale_weight_edge(image, weight, options) # set the negative weights to zero (those are flags) for weight in weights: weight_max = weight.max() weight[weight0).astype(float) unmask_sum = (mask*abs(model)).sum() sum = abs(model).sum() fractions.append(unmask_sum/sum) return np.array(fractions) except: if fatal_errors: raise else: return np.repeat(-1.0, len(models)) def get_edge_info(models, weights, fatal_errors=False): try: means = [] stdevs = [] for (weight,model) in zip(weights, models): mask = (weight>0).astype(float) mask[1:-1, 1:-1] = 0.0 edges = model[mask>0.0] means.append(edges.mean()) stdevs.append(edges.std()) return np.array(means), np.array(stdevs) except: if fatal_errors: raise else: return np.repeat(-1.0, len(models)), np.repeat(-1.0, len(models)) def run_standard_analysis(options, meds, args, iobj, ra, dec, flag, index, output): try: result_status, i3_result, n_exposure, psfs, extra_info = meds.im3shape_analyze(options, args.psf_input_method, iobj, rescale_weight=args.rescale_weight, first_exposure=args.first_exp, max_nexposures=args.max_nexp,use_nearest_pix_mask=args.use_nearest_pix_mask, reject_outliers=(not args.no_reject_outliers)) #We try to do a fairly full output for im3shape errors, though these are rare. except Exception as error: if args.fatal_errors: raise print "An error occurred processing one object. Info below" print error.message exc_type, exc_obj, exc_tb = sys.exc_info() fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1] print exc_type, fname, exc_tb.tb_lineno return flux_fractions = get_unmasked_flux_fractions(extra_info['models'], extra_info['weights'], args.fatal_errors) img_flux_fractions = get_unmasked_flux_fractions(extra_info['images'], extra_info['weights'], args.fatal_errors) uber_flux_fractions = get_unmasked_flux_fractions(extra_info['models'], extra_info['uber_mask'], args.fatal_errors) img_uber_flux_fractions = get_unmasked_flux_fractions(extra_info['images'], extra_info['uber_mask'], args.fatal_errors) model_edge_mu, model_edge_sigma = get_edge_info(extra_info['models'], extra_info['weights'], args.fatal_errors) edge_mu, edge_sigma = get_edge_info(extra_info['images'], extra_info['weights'], args.fatal_errors) fwhm_psfs, fwhm_gals, rgpp_rps, psf_params, \ hsm_fwhm_psfs, hsm_gal_params, hsm_rgpp_rps, hsm_psf_params = \ extract_psf_outputs(psfs, n_exposure, options, meds, iobj, i3_result, fatal=args.fatal_errors) #Save output. # Append info on the PSFs to the main output file or write it out in an extra file . #This is getting a bit crazy - basically this has to match the list of headers we generated earlier. extra = [] for i in xrange(MAX_N_EXPOSURE): if i < n_exposure: for j in xrange(len(psf_params[0])): extra.append(psf_params[i][j]) extra.append(fwhm_gals[i]) extra.append(rgpp_rps[i]) for j in xrange(len(hsm_psf_params[0])): extra.append(hsm_psf_params[i][j]) for j in xrange(len(hsm_gal_params[0])): extra.append(hsm_gal_params[i][j]) extra.append(hsm_rgpp_rps[i]) extra.append(flux_fractions[i]) extra.append(img_flux_fractions[i]) extra.append(uber_flux_fractions[i]) extra.append(img_uber_flux_fractions[i]) extra.append(edge_mu[i]) extra.append(edge_sigma[i]) extra.append(model_edge_mu[i]) extra.append(model_edge_sigma[i]) extra.append(extra_info['orig_rows'][i]) extra.append(extra_info['orig_cols'][i]) extra.append(extra_info['source_path_indices'][i]) if args.psf_input_method!='great_des': extra.append(extra_info['exposure_names'][i]) else: extra.append('great-des') print 'warning: great-des: using dummy exposure_names' extra.append(extra_info['ccds'][i]) else: for j in xrange(len(psf_params[0])): extra.append(-1.0) #psf_params extra.append(-1.0) #fwhm_gals extra.append(-1.0) #rgpp_rp # --- Add stuff from hsm FWHM estimation for j in xrange(len(hsm_psf_params[0])): extra.append(-1.0) #hsm_psf for j in xrange(len(hsm_gal_params[0])): extra.append(-1.0) #hsm_gal extra.append(-1.0) #hsm_rgpp_rp extra.append(-1.0) #flux_fractions extra.append(-1.0) #img_flux_fractions extra.append(-1.0) #uber_flux_fractions extra.append(-1.0) #img_uber_flux_fractions extra.append(-1.0) #edge_mu extra.append(-1.0) #edge_sigma extra.append(-1.0) #model_edge_mu extra.append(-1.0) #model_edge_sigma extra.append(-1.0) #orig_row extra.append(-1.0) #orig_col extra.append(-1) #source_path extra.append('N/A') #exposure name extra.append(-1) #ccd extra.append(n_exposure) extra.append(extra_info['stamp_size']) extra.append(extra_info['mask_fraction']) extra.append(extra_info['nreject']) extra.append(int(flag[index])) output.write_row(i3_result, options, psfs, ra[index], dec[index], extra_output=extra) def run_psf_only_analysis(options, meds, args, iobj, ra, dec, index, output): # In this case analyse PSF image only, i.e. we do not run im3shape on galaxy images psfs = meds.get_PSFs( options, args.psf_input_method, iobj) # Try to compute PSF parameters. Note, that we do need # options variable handed in for information on upsampling settings try: psf_params = get_psf_params(psfs, options) except Exception as emsg: logger.error('PSF parameters calculation failed for object number %d with messge \n %s' % (i3_result.identifier, emsg)) psf_params = [[-1.0, -1.0, -1.0, -1.0, -1.0] for i in xrange(len(psfs))] # Append info on the PSFs to the main output file or write it out in an extra file extra = [] for i in xrange(MAX_N_EXPOSURE): if i < len(psfs): for j in xrange(len(psf_params[0])): extra.append(psf_params[i][j]) else: for j in xrange(len(psf_params[0])): extra.append(-1.0) # No write output to file output.write_row(None, options, psfs, ra[index], dec[index], extra_output=extra, id=iobj) def main(args): if args.no_clobber and os.path.exists(args.filename_output_result): logger.info('Skipping existing output file %s'%args.filename_output_result) return # load the options file options=py3shape.Options(args.file_input_ini) options.validate() options.save_output=False if args.extra_options is not None: for opt in args.extra_options: key,val = opt.split('=') setattr(options, key, val) # open the output files, as an output object output = py3shape.output.Output(args.filename_output_result, options, logger=logger) # add a few extra convenience headers, for tracking extra_lines=[ "meds: %s" % args.file_input_meds, "cat: %s" % args.file_selected_objects, "ini: %s" % args.file_input_ini, "out: %s" % args.filename_output_result, "first: %d" % args.first, "count: %d" % args.count, ] # write the header extra_cols = [] for i in xrange(MAX_N_EXPOSURE): extra_cols.append("psf_fwhm_%d"%(i+1)) extra_cols.append("psf_e1_%d"%(i+1)) extra_cols.append("psf_e2_%d"%(i+1)) extra_cols.append("psf_e1_sky_%d"%(i+1)) extra_cols.append("psf_e2_sky_%d"%(i+1)) if not args.analyse_psf_only: extra_cols.append("gal_fwhm_%d"%(i+1)) extra_cols.append("rgpp_rp_%d"%(i+1)) extra_cols.append("hsm_psf_moments_sigma_%d"%(i+1)) extra_cols.append("hsm_psf_observed_shape_e1_%d"%(i+1)) extra_cols.append("hsm_psf_observed_shape_e2_%d"%(i+1)) extra_cols.append("hsm_psf_observed_shape_e1_sky_%d"%(i+1)) extra_cols.append("hsm_psf_observed_shape_e2_sky_%d"%(i+1)) extra_cols.append("hsm_psf_moments_rho4_%d"%(i+1)) extra_cols.append("hsm_psf_moments_status_%d"%(i+1)) extra_cols.append("hsm_psf_moments_n_iter_%d"%(i+1)) if not args.analyse_psf_only: extra_cols.append("hsm_gal_moments_sigma_%d"%(i+1)) extra_cols.append("hsm_gal_moments_rho4_%d"%(i+1)) extra_cols.append("hsm_gal_moments_status_%d"%(i+1)) extra_cols.append("hsm_gal_moments_n_iter_%d"%(i+1)) extra_cols.append("hsm_rgpp_rp_%d"%(i+1)) extra_cols.append("unmasked_flux_frac_%d"%(i+1)) extra_cols.append("unmasked_img_flux_frac_%d"%(i+1)) extra_cols.append("un_uber_flux_frac_%d"%(i+1)) extra_cols.append("un_uber_img_flux_frac_%d"%(i+1)) extra_cols.append("edge_mu_%d"%(i+1)) extra_cols.append("edge_sigma_%d"%(i+1)) extra_cols.append("model_edge_mu_%d"%(i+1)) extra_cols.append("model_edge_sigma_%d"%(i+1)) extra_cols.append("orig_row_%d"%(i+1)) extra_cols.append("orig_col_%d"%(i+1)) extra_cols.append("source_path_%d"%(i+1)) extra_cols.append("exposure_name_%d"%(i+1)) extra_cols.append("ccd_%d"%(i+1)) extra_cols.append("n_exposure") extra_cols.append("stamp_size") extra_cols.append("mask_fraction") extra_cols.append("nreject") extra_cols.append("flag") output.write_header(include_radec=True, extra_cols=extra_cols, extra_lines=extra_lines, analyse_psf_only=args.analyse_psf_only) # open the meds file logger.info('opening meds file %s' % args.file_input_meds) if args.replace_psf_path: replace_psf_path = args.replace_psf_path.split(',') else: replace_psf_path = None meds = I3MEDS(args.file_input_meds, replace_psf_path=replace_psf_path, psfex_rerun_version=args.psfex_rerun_version) meds.mask_dilation = args.dilate_mask # optionally save the list of exposure paths that went into this meds file. # the objects themselves then save an index into this for their exposures if args.save_source_paths: meds.i3_save_source_paths(args.save_source_paths) n_total = meds.size # Default is to select all the objects in the MEDS file if args.file_selected_objects == None: rowids = range(0,n_total) ra = np.zeros(n_total) dec = np.zeros(n_total) flag = np.zeros(n_total) #Alternative is to read a very simple catalog else: try: selected_objects = np.atleast_2d(np.loadtxt(args.file_selected_objects)).T rowids,ximg,yimg,ra,dec,flag = selected_objects[:6] except Exception,e: logger.error('Error in opening selection file %s, format should be (rowid,ximg,yimg,ra,dec,flag)' % args.file_selected_objects) return FLAG_ERROR n_selected = len(rowids) # check if first and last object to process are within the correct range n_total = meds.size n_selected = len(rowids) #Simple error check for silly input if args.first >= n_total: logger.error('first object to process is exceeding the number of possible objects') return FLAG_ERROR last_to_process = min(args.first+args.count-1, n_total) # loop from first to last logger.info('running on objects with IDs from %d to %d' % (args.first,last_to_process)) if args.psf_input_method.lower() == 'single': logger.warning("Using PSF image %s for all galaxies" % args.filename_psf_image) if args.rerun_post: #Load in im3shape output file. Requires astropy import astropy.table logger.info("\nLoading table to re-run parameters. NO MINIMIZATION WILL BE DONE!\n") table = astropy.table.Table.read(args.rerun_post, format='fits') table.sort('identifier') options.minimizer_loops=-1 options.use_computed_start = False model_struct = getattr(py3shape.structs, options['model_name']+"_parameter_set") # set up the counter n_processed_objects = 0 #Loop through all the selected objects. #We will only run on a subset of them, though. for index in xrange(n_selected): #Determine the original MEDS ID for this selection row. #This is the value that we split up work onto, because #we do not know in advance how many galaxies will make it #through the cuts. #Downside: different numbers of gals to run for each job iobj = int(rowids[index]) #Skip if this one is not our job if iobj < args.first or iobj>last_to_process: continue if args.rerun_post: raise ValueError("Need to re-code this bit as IDs are now coadd_object_ids") ID = meds.get_number(iobj) row_index = np.searchsorted(table['identifier'], ID) if row_index>len(table) or table['identifier'][row_index]!=ID: logger.error("Failed to find re-run entry for object ID=%r, iobj=%r"%(ID,iobj)) continue row = dict(zip(table.colnames, table[row_index])) if args.unblind_rerun: search_path = os.environ['UCL_DES_SHEAR_DIR']+'/des_post' sys.path.append(search_path) import blind_catalog f = blind_catalog.get_factor(blind_catalog.code_phrase) if not (f>0.9 and f<1): raise ValueError("Blind factor broken") row['e1'] = row['e1'] / f row['e2'] = row['e2'] / f #find the row for this object in the results file for (name, _) in model_struct._fields_: val = row[name] options[name+"_start"] = row[name] #Some info logger.info('running on row=%10d from selection, MEDS row %10d, coadd number %10d' % (index,iobj,meds.get_number(iobj))) # Go to functions to do the main job. if not args.analyse_psf_only: run_standard_analysis(options, meds, args, iobj, ra, dec, flag, index, output) else: run_psf_only_analysis(options, meds, args, iobj, ra, dec, index, output) n_processed_objects += 1 #Tidy up output.close() logger.info('done analysing %d objects' % n_processed_objects) # Set up and parse the command line arguments using the nice Python argparse module description = "Run im3shape on a meds image using multiple exposures fit." parser = argparse.ArgumentParser(description=description, add_help=True) parser.add_argument('file_input_meds', type=str, help='the input meds FITS file') parser.add_argument('file_input_ini', type=str, help='the name of the existing ini file to be used with im3shape') # parser.add_argument('working_directory', type=str, help='path of working directory. If not existing, it will be created') parser.add_argument('filename_output_result', type=str, help='name of the results file with im3shape output in sky coordinates (one line per galaxy)') parser.add_argument('--save-source-paths', type=str, default="", help='Save the list that maps source path indices to source paths to this filename') parser.add_argument('-c','--file_selected_objects', default=None,type=str, help='the name of the txt file specifying the objects on which to run') parser.add_argument("--first", type=int, default=0, help='The first object BY MEDS OBJECT NUMBER to include in the output') parser.add_argument("--count", type=int, default=10000000000, help='Max index in MEDS file to be run. This is not the number actually run!') parser.add_argument("--max-nexp", type=int, default=None, help='Max number of exposures used per object. Subsequent ones are ignored.') parser.add_argument("--first-exp", type=int, default=1, help='First exposure to use, starting from 1. Ones before this are ignored.') parser.add_argument('-v', '--verbosity', type=int, action='store', default=2, choices=(0, 1, 2, 3 ), help='integer verbosity level: min=0, max=3 [default=2]') parser.add_argument('--psf_input_method', type=str, action='store', default='psfex', choices=('psfex','shapelet','single','great_des', 'single_moffat', 'none'), help='PSF input methods') parser.add_argument('--filename_psf_image', type=str, action='store', default=None, help='if supplied, then all galaxies analysed here will use PSF from this file') parser.add_argument('--rerun-post', type=str, action='store', default=None, help='if supplied, do not perform minimization; use the starting point from the specified file') parser.add_argument('--unblind-rerun', type=bool, action='store', default=False, help='Unblind the ellipticity estimates before re-running on them. Only relevant if --rerun-post is set') parser.add_argument('--no-reject-outliers', type=bool, action='store', default=False, help='Switch off outlier rejection from median image') parser.add_argument('--analyse_psf_only', type=bool, action='store', default=False, help='if set to true, then only PSF images are analysed and no galaxy images! Output file will contain info on PSF only.') parser.add_argument('--fatal-errors', action='store_true', default=False, help='errors are fatal if set') parser.add_argument('--dilate_mask', type=int, default=0, help='Enlarge completely masked out regions by this many pixels') parser.add_argument('--rescale-weight', action='store_true', default=False, help='scale the weights to match the std dev of the edge of the image') parser.add_argument('--no-clobber', action='store_true', default=False, help='If the output file already exists, do not run the analysis') parser.add_argument('--filename_psf_single_moffat', type=str, action='store', default=None, help='if supplied, then all galaxies analysed here will use PSF from this file, containing single line with beta fwhm e1 e2') parser.add_argument('--replace_psf_path', type=str, action='store', default=None, help='Comma separated path replacement for PSF') parser.add_argument('--psfex-rerun-version', type=str, action='store', default=None, help='PSF rerun version number') parser.add_argument('--use_nearest_pix_mask', action='store_true', default=False, help='Weight all pixels that are closer to the seg map of a neighbour than to the seg map of the target to zero') parser.add_argument('-p', '--option', dest='extra_options', action='append', help='Additional options to be set. Can specify more than once in form -p option=value. Overrides ini file') if __name__=="__main__": args = parser.parse_args() logging_levels = { 0: logging.CRITICAL, 1: logging.WARNING, 2: logging.INFO, 3: logging.DEBUG } logging_level = logging_levels[args.verbosity] logging.basicConfig(format="%(message)s", level=logging_level, stream=sys.stdout) logger.setLevel(logging_level) logger.info('running main') main(args) sys.exit(0)