diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..728213a --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,54 @@ +name: tests + +on: + push: + branches: + - master + pull_request: null + +jobs: + tests: + name: tests + strategy: + matrix: + pyver: [3.6, 3.7, 3.8] + + runs-on: "ubuntu-latest" + + steps: + - uses: actions/checkout@v2 + + - uses: conda-incubator/setup-miniconda@v2 + with: + python-version: ${{ matrix.pyver }} + channels: conda-forge,defaults + channel-priority: strict + show-channel-urls: true + + - name: configure conda and install code + shell: bash -l {0} + run: | + conda config --set always_yes yes + conda install -q mamba + + # requirements for meds + mamba install -q --file requirements.txt + + # installs for testing + mamba install -q \ + pip \ + setuptools \ + flake8 \ + pytest-cov + + python -m pip install -e . + + - name: lint + shell: bash -l {0} + run: | + flake8 meds + + - name: test + shell: bash -l {0} + run: | + pytest -v --cov=meds --cov-report term-missing meds diff --git a/meds/__init__.py b/meds/__init__.py index c1eb3c1..ae4814e 100644 --- a/meds/__init__.py +++ b/meds/__init__.py @@ -1,3 +1,4 @@ +# flake8: noqa from .defaults import __version__ from . import meds diff --git a/meds/bounds.py b/meds/bounds.py index 6bdd5b6..baa3e4e 100644 --- a/meds/bounds.py +++ b/meds/bounds.py @@ -5,10 +5,10 @@ class Bounds(object): Based on Bounds from deswl_shapelets """ def __init__(self, rowmin, rowmax, colmin, colmax): - self.rowmin=rowmin - self.rowmax=rowmax - self.colmin=colmin - self.colmax=colmax + self.rowmin = rowmin + self.rowmax = rowmax + self.colmin = colmin + self.colmax = colmax def contains_point(self, row, col): """ @@ -33,11 +33,10 @@ def contains_points(self, rows, cols): rows,cols: arrays The positions to check """ - return ( (rows <= self.rowmax) - & (rows >= self.rowmin) - & (cols <= self.colmax) - & (cols >= self.colmin) ) - + return ((rows <= self.rowmax) & + (rows >= self.rowmin) & + (cols <= self.colmax) & + (cols >= self.colmin)) def contains_bounds(self, bounds): """ @@ -48,7 +47,7 @@ def contains_bounds(self, bounds): bounds: Bounds The bounding box to check """ - assert isinstance(bounds,Bounds) + assert isinstance(bounds, Bounds) return (bounds.rowmin >= self.rowmin and bounds.rowmax <= self.rowmax and bounds.colmin >= self.colmin @@ -63,12 +62,12 @@ def intersects_bounds(self, bounds): bounds: Bounds The bounding box to check """ - assert isinstance(bounds,Bounds) + assert isinstance(bounds, Bounds) - return ( not (bounds.rowmin >= self.rowmax) + return (not (bounds.rowmin >= self.rowmax) and not (bounds.rowmax <= self.rowmin) and not (bounds.colmin >= self.colmax) - and not (bounds.colmax <= self.colmin) ) + and not (bounds.colmax <= self.colmin)) def expand_point(self, row, col): """ @@ -114,8 +113,6 @@ def expand_points(self, rows, cols): if colmax > self.colmax: self.colmax = colmax - - def expand_bounds(self, bounds): """ expand the bounds to include the input Bounds @@ -125,7 +122,7 @@ def expand_bounds(self, bounds): bounds: Bounds The bounding box to check """ - assert isinstance(bounds,Bounds) + assert isinstance(bounds, Bounds) if bounds.rowmin < self.rowmin: self.rowmin = bounds.rowmin @@ -138,8 +135,6 @@ def expand_bounds(self, bounds): self.colmax = bounds.colmax def __repr__(self): - mess="Bounds(rowmin=%g, rowmax=%g, colmin=%g, colmax=%g)" + mess = "Bounds(rowmin=%g, rowmax=%g, colmin=%g, colmax=%g)" mess = mess % (self.rowmin, self.rowmax, self.colmin, self.colmax) return mess - - diff --git a/meds/coadd.py b/meds/coadd.py index 427d04d..bfe1fc4 100644 --- a/meds/coadd.py +++ b/meds/coadd.py @@ -20,26 +20,19 @@ import numpy as np import fitsio from . import util -import yaml from . import maker from .meds import reject_outliers from .defaults import __version__ -try: - xrange=xrange -except: - xrange=range - - # flags for why a coadd was not made -NO_GOOD_CUTOUTS=2**0 +NO_GOOD_CUTOUTS = 2 ** 0 # flags for why a cutout wasn't included -CUTOUT_HITS_EDGE=2**1 -CUTOUT_MASK_FRAC=2**2 -CUTOUT_CENTRAL_MASKED=2**3 -CUTOUT_FULL_EDGE_MASKED=2**4 -CUTOUT_ALL_WEIGHT_ZERO=2**5 +CUTOUT_HITS_EDGE = 2 ** 1 +CUTOUT_MASK_FRAC = 2 ** 2 +CUTOUT_CENTRAL_MASKED = 2 ** 3 +CUTOUT_FULL_EDGE_MASKED = 2 ** 4 +CUTOUT_ALL_WEIGHT_ZERO = 2 ** 5 """ FLAG_MAP={ @@ -66,20 +59,21 @@ FLAG_MAP_INV[FLAG_MAP[flagname]] = flagname """ + class MEDSCoaddMaker(maker.MEDSMaker): """ Convert an existing MEDS file to a postage-stamp coadd MEDS file using the input coadder """ - def __init__(self, config, coadder, tmpdir=None): + def __init__(self, config, coadder, tmpdir=None): - self.tmpdir=tmpdir + self.tmpdir = tmpdir self._load_config(config) self._set_extra_config() - self.coadder=coadder + self.coadder = coadder # just a reference self.m = coadder.get_original_meds() @@ -88,7 +82,7 @@ def __init__(self, config, coadder, tmpdir=None): # in fits they are still arrays self._set_obj_data() - self.total_pixels = (self.m['box_size']**2).sum() + self.total_pixels = (self.m["box_size"] ** 2).sum() self._set_image_info() self._set_meta() @@ -96,11 +90,10 @@ def __init__(self, config, coadder, tmpdir=None): def _load_config(self, config): - super(MEDSCoaddMaker,self)._load_config(config) - for t in ['ormask','noise']: - if t not in self['cutout_types']: - self['cutout_types'] += [t] - + super(MEDSCoaddMaker, self)._load_config(config) + for t in ["ormask", "noise"]: + if t not in self["cutout_types"]: + self["cutout_types"] += [t] def write(self, filename, obj_range=None): """ @@ -108,8 +101,8 @@ def write(self, filename, obj_range=None): """ print("opening output MEDS file: '%s'" % filename) - with fitsio.FITS(filename,'rw',clobber=True) as fits: - self.fits=fits + with fitsio.FITS(filename, "rw", clobber=True) as fits: + self.fits = fits self._write_image_info() self._write_metadata() @@ -121,7 +114,7 @@ def write(self, filename, obj_range=None): # we fill this in as we go self._write_object_data() - print('output is in:',filename) + print("output is in:", filename) def _write_cutouts(self, obj_range=None): """ @@ -129,15 +122,15 @@ def _write_cutouts(self, obj_range=None): data """ if obj_range is None: - obj_range=[0,self.m.size] + obj_range = [0, self.m.size] - self.last_row_start=0 - self.last_size=0 - self.last_psf_row_start=0 - self.last_psf_size=0 + self.last_row_start = 0 + self.last_size = 0 + self.last_psf_row_start = 0 + self.last_psf_size = 0 - for iobj in xrange(obj_range[0], obj_range[1]): - print("%d %d/%d" % (iobj,iobj+1,obj_range[1]-1)) + for iobj in range(obj_range[0], obj_range[1]): + print("%d %d/%d" % (iobj, iobj + 1, obj_range[1] - 1)) self._write_object_cutouts(iobj) @@ -145,37 +138,37 @@ def _write_object_cutouts(self, iobj): """ make the coadd and write all the data """ - icut=0 - d=self.obj_data + icut = 0 + d = self.obj_data coadd_obs, flags = self.coadder.get_coadd_obs(iobj) - d['coadd_flags'][iobj] = flags + d["coadd_flags"][iobj] = flags if coadd_obs is not None: - assert d['box_size'][iobj]==coadd_obs.image.shape[0] + assert d["box_size"][iobj] == coadd_obs.image.shape[0] - d['ncoadd'][iobj] = coadd_obs.meta['ncoadd'] + d["ncoadd"][iobj] = coadd_obs.meta["ncoadd"] - imflags=coadd_obs.meta['imflags'] - norig=len(imflags) - d['orig_flags'][iobj,0:norig] = imflags - d['start_row'][iobj,icut]=\ - self.last_row_start+self.last_size - d['psf_start_row'][iobj,icut]=\ - self.last_psf_row_start+self.last_psf_size + imflags = coadd_obs.meta["imflags"] + norig = len(imflags) + d["orig_flags"][iobj, 0:norig] = imflags + d["start_row"][iobj, icut] = self.last_row_start + self.last_size + d["psf_start_row"][iobj, icut] = ( + self.last_psf_row_start + self.last_psf_size + ) j = coadd_obs.jacobian pj = coadd_obs.psf.jacobian - d['cutout_row'][iobj,icut] = j.row0 - d['cutout_col'][iobj,icut] = j.col0 - d['dudrow'][iobj,icut] = j.dudrow - d['dudcol'][iobj,icut] = j.dudcol - d['dvdrow'][iobj,icut] = j.dvdrow - d['dvdcol'][iobj,icut] = j.dvdcol - - d['psf_box_size'][iobj] = coadd_obs.psf.image.shape[0] - d['psf_cutout_row'][iobj,icut] = pj.row0 - d['psf_cutout_col'][iobj,icut] = pj.col0 - - for cutout_type in self['cutout_types']+['psf']: + d["cutout_row"][iobj, icut] = j.row0 + d["cutout_col"][iobj, icut] = j.col0 + d["dudrow"][iobj, icut] = j.dudrow + d["dudcol"][iobj, icut] = j.dudcol + d["dvdrow"][iobj, icut] = j.dvdrow + d["dvdcol"][iobj, icut] = j.dvdcol + + d["psf_box_size"][iobj] = coadd_obs.psf.image.shape[0] + d["psf_cutout_row"][iobj, icut] = pj.row0 + d["psf_cutout_col"][iobj, icut] = pj.col0 + + for cutout_type in self["cutout_types"] + ["psf"]: self._write_cutout( iobj, icut, @@ -183,87 +176,88 @@ def _write_object_cutouts(self, iobj): cutout_type, ) - self.last_row_start = d['start_row'][iobj,icut] - self.last_size = d['box_size'][iobj]**2 - self.last_psf_row_start = d['psf_start_row'][iobj,icut] - self.last_psf_size = d['psf_box_size'][iobj]**2 + self.last_row_start = d["start_row"][iobj, icut] + self.last_size = d["box_size"][iobj] ** 2 + self.last_psf_row_start = d["psf_start_row"][iobj, icut] + self.last_psf_size = d["psf_box_size"][iobj] ** 2 - d['ncutout'][iobj] = 1 + d["ncutout"][iobj] = 1 def _write_cutout(self, iobj, icut, coadd_obs, cutout_type): cutout_hdu = self._get_cutout_hdu(cutout_type) - if cutout_type=='psf': + if cutout_type == "psf": im_data = coadd_obs.psf.image - sname='psf_start_row' + sname = "psf_start_row" else: - sname='start_row' - if cutout_type=='image': + sname = "start_row" + if cutout_type == "image": im_data = coadd_obs.image - elif cutout_type=='weight': + elif cutout_type == "weight": im_data = coadd_obs.weight - elif cutout_type=='seg': + elif cutout_type == "seg": im_data = coadd_obs.seg - elif cutout_type=='bmask': + elif cutout_type == "bmask": im_data = coadd_obs.bmask - elif cutout_type=='ormask': + elif cutout_type == "ormask": im_data = coadd_obs.ormask - elif cutout_type=='noise': + elif cutout_type == "noise": im_data = coadd_obs.noise else: raise ValueError("bad cutout type: '%s'" % cutout_type) - start_row=self.obj_data[sname][iobj,icut] + start_row = self.obj_data[sname][iobj, icut] cutout_hdu.write(im_data.ravel(), start=start_row) def _set_obj_data(self): - nmax=2 + nmax = 2 # subtract 1 for coadd - nmax_orig=max(self.m['ncutout'].max()-1,2) - extra_fields=[ - ('number','i8'), - ('psf_box_size','i4'), - ('psf_cutout_row','f8',nmax), - ('psf_cutout_col','f8',nmax), - ('psf_start_row','i8',nmax), - ('orig_ncutout','i4'), - ('orig_flags','i4',nmax_orig), - ('ncoadd','i4'), # how many images got coadded - ('coadd_flags','i4'), + nmax_orig = max(self.m["ncutout"].max() - 1, 2) + extra_fields = [ + ("number", "i8"), + ("psf_box_size", "i4"), + ("psf_cutout_row", "f8", nmax), + ("psf_cutout_col", "f8", nmax), + ("psf_start_row", "i8", nmax), + ("orig_ncutout", "i4"), + ("orig_flags", "i4", nmax_orig), + ("ncoadd", "i4"), # how many images got coadded + ("coadd_flags", "i4"), ] - d=util.get_meds_output_struct( + d = util.get_meds_output_struct( self.m.size, nmax, extra_fields=extra_fields, ) - cat=self.m.get_cat() + cat = self.m.get_cat() - d['ncutout']=0 - d['orig_ncutout'] = self.m['ncutout']-1 # -1 for original coadd - d['ncoadd']=0 - d['file_id']=0 - d['coadd_flags']=0 + d["ncutout"] = 0 + d["orig_ncutout"] = self.m["ncutout"] - 1 # -1 for original coadd + d["ncoadd"] = 0 + d["file_id"] = 0 + d["coadd_flags"] = 0 - copy_fields=['id','box_size','ra','dec','number'] + copy_fields = ["id", "box_size", "ra", "dec", "number"] for f in copy_fields: d[f] = cat[f] # just copy coadd coadd_fields = [ - 'orig_row','orig_col', + "orig_row", + "orig_col", # needed by mof - 'orig_start_row','orig_start_col', + "orig_start_row", + "orig_start_col", ] for f in coadd_fields: - d[f][:,0] = cat[f][:,0] - + d[f][:, 0] = cat[f][:, 0] - self.obj_data=d - self.ncutout_max=nmax + self.obj_data = d + self.ncutout_max = nmax def _set_psf_layout(self): """ @@ -271,7 +265,7 @@ def _set_psf_layout(self): """ raise NotImplementedError("implement psf specifics in a sub class") - ''' + """ obj_data = self.obj_data max_psf_size = self['max_psf_size'] @@ -280,14 +274,14 @@ def _set_psf_layout(self): self.total_psf_pixels = max_npixels_per*nobj # fake the data for later self.psf_data=1 - ''' + """ def _set_image_info(self): """ fake image info """ - #self.image_info=util.get_image_info_struct(1, 10) - self.image_info=self.m.get_image_info()[0:1] + # self.image_info=util.get_image_info_struct(1, 10) + self.image_info = self.m.get_image_info()[0:1] def _set_meta(self): """ @@ -295,56 +289,53 @@ def _set_meta(self): """ import esutil as eu - ometa=self.m.get_meta() + ometa = self.m.get_meta() - dt=[] - nv=np.__version__ - ev=eu.__version__ - fv=fitsio.__version__ - mv=__version__ - mfv=maker.MEDS_FMT_VERSION + dt = [] + nv = np.__version__ + ev = eu.__version__ + fv = fitsio.__version__ + mv = __version__ + mfv = maker.MEDS_FMT_VERSION dt += [ - ('numpy_version','S%d' % len(nv)), - ('esutil_version','S%d' % len(ev)), - ('fitsio_version','S%d' % len(fv)), - ('meds_version','S%d' % len(mv)), - ('meds_fmt_version','S%d' % len(mfv)), + ("numpy_version", "S%d" % len(nv)), + ("esutil_version", "S%d" % len(ev)), + ("fitsio_version", "S%d" % len(fv)), + ("meds_version", "S%d" % len(mv)), + ("meds_fmt_version", "S%d" % len(mfv)), ] - names=[ d[0] for d in dt] + names = [d[0] for d in dt] for d in ometa.dtype.descr: - n=d[0] + n = d[0] if n not in names: - dt.append( d ) + dt.append(d) meta = np.zeros(1, dtype=dt) eu.numpy_util.copy_fields(ometa, meta) - meta['numpy_version'] = nv - meta['esutil_version'] = ev - meta['fitsio_version'] = fv - meta['meds_version'] = mv - meta['meds_fmt_version'] = mfv + meta["numpy_version"] = nv + meta["esutil_version"] = ev + meta["fitsio_version"] = fv + meta["meds_version"] = mv + meta["meds_fmt_version"] = mfv + + self.meta_data = meta + - self.meta_data=meta - class MEDSCoadder(dict): """ generate postage stamp coadds and PSF coadds from the input MEDS object """ - def __init__(self, - config, - meds_obj, - psfmap, - seed, - make_plots=False): + + def __init__(self, config, meds_obj, psfmap, seed, make_plots=False): self._set_config(config) - self.make_plots=make_plots + self.make_plots = make_plots self._set_bad_flags() self._set_edge_flags() - self.rng=np.random.RandomState(seed) + self.rng = np.random.RandomState(seed) self.m = meds_obj self.psfmap = psfmap @@ -362,24 +353,22 @@ def get_coadd_obs(self, iobj): """ import psc - flags=0 + flags = 0 obslist, imflags = self._get_obslist(iobj) - - ncoadd=len(obslist) - if ncoadd==0: + + ncoadd = len(obslist) + if ncoadd == 0: flags |= NO_GOOD_CUTOUTS print(" nothing to coadd") return None, flags - - bmask=obslist[0].bmask*0 - ormask=obslist[0].bmask*0 + + bmask = obslist[0].bmask * 0 + ormask = obslist[0].bmask * 0 for obs in obslist: ormask |= obs.bmask - coadder=psc.Coadder( - obslist, - jacobian=self.target_jacobian, - **self['coadd'] + coadder = psc.Coadder( + obslist, jacobian=self.target_jacobian, **self["coadd"], ) coadd_obs = coadder.get_coadd() coadd_obs.bmask = bmask @@ -387,15 +376,16 @@ def get_coadd_obs(self, iobj): coadd_obs.seg = self._get_interpolated_coadd_seg(iobj) if self.make_plots: - original_coadd=self.m.get_cutout(iobj, 0) - original_seg=self.m.get_cutout(iobj, 0,type='seg') - original_wt=self.m.get_cutout(iobj, 0,type='weight') - self._show_coadd(iobj, coadd_obs, - original_coadd, original_wt, original_seg) - - meta={ - 'ncoadd':ncoadd, - 'imflags':imflags, + original_coadd = self.m.get_cutout(iobj, 0) + original_seg = self.m.get_cutout(iobj, 0, type="seg") + original_wt = self.m.get_cutout(iobj, 0, type="weight") + self._show_coadd( + iobj, coadd_obs, original_coadd, original_wt, original_seg, + ) + + meta = { + "ncoadd": ncoadd, + "imflags": imflags, } coadd_obs.update_meta_data(meta) return coadd_obs, flags @@ -405,116 +395,118 @@ def _get_interpolated_coadd_seg(self, iobj): return self.m._interpolate_coadd_seg_image(iobj, jmatrix) def _get_jacobian_matrix(self): - j=self.target_jacobian + j = self.target_jacobian return np.matrix( - [ [j.dudrow, j.dudcol], - [j.dvdrow, j.dvdcol], ] + [ + [j.dudrow, j.dudcol], + [j.dvdrow, j.dvdcol], + ] ) def _set_config(self, config): self.update(config) - self['coadd'] = self.get('coadd',{}) - - self['dither_psfs'] = self['coadd'].pop('dither_psfs',True) + self["coadd"] = self.get("coadd", {}) + self["dither_psfs"] = self["coadd"].pop("dither_psfs", True) def _set_bad_flags(self): - if 'bmask_flags' not in self: - bad_flags=None + if "bmask_flags" not in self: + bad_flags = None else: - bad_flags = sum(self['bmask_flags']) + bad_flags = sum(self["bmask_flags"]) - self.bad_flags=bad_flags + self.bad_flags = bad_flags def _set_edge_flags(self): - if 'edge_flags' not in self: - edge_flags=None + if "edge_flags" not in self: + edge_flags = None else: - edge_flags = sum(self['edge_flags']) + edge_flags = sum(self["edge_flags"]) - self.edge_flags=edge_flags + self.edge_flags = edge_flags def _set_target_jacobian(self): import ngmix + # center doesn't matter - w,= np.where(self.m['ncutout'] > 1) - dudrow = np.median(self.m['dudrow'][w,1]) - dudcol = np.median(self.m['dudcol'][w,1]) - dvdrow = np.median(self.m['dvdrow'][w,1]) - dvdcol = np.median(self.m['dvdcol'][w,1]) + (w,) = np.where(self.m["ncutout"] > 1) + dudrow = np.median(self.m["dudrow"][w, 1]) + dudcol = np.median(self.m["dudcol"][w, 1]) + dvdrow = np.median(self.m["dvdrow"][w, 1]) + dvdcol = np.median(self.m["dvdcol"][w, 1]) - self.target_jacobian=ngmix.Jacobian( - row=15, col=15, + self.target_jacobian = ngmix.Jacobian( + row=15, + col=15, dudrow=dudrow, dudcol=dudcol, dvdrow=dvdrow, dvdcol=dvdcol, ) - def _get_obslist(self, iobj): import ngmix - m=self.m + m = self.m - make_plots=self.make_plots + make_plots = self.make_plots - obslist=ngmix.ObsList() + obslist = ngmix.ObsList() imlist = m.get_cutout_list(iobj) - if len(imlist)==1: + if len(imlist) == 1: # there is just the coadd; return the # empty obslist - imflags=None + imflags = None return obslist, imflags - wtlist = m.get_cutout_list(iobj, type='weight') - bmlist = m.get_cutout_list(iobj, type='bmask') - + wtlist = m.get_cutout_list(iobj, type="weight") + bmlist = m.get_cutout_list(iobj, type="bmask") - imlist=imlist[1:] - wtlist=wtlist[1:] - bmlist=bmlist[1:] + imlist = imlist[1:] + wtlist = wtlist[1:] + bmlist = bmlist[1:] - if self['reject_outliers']: + if self["reject_outliers"]: nreject = reject_outliers( imlist, wtlist, ) if nreject > 0: - print(" rejected",nreject) - obslist.meta['nreject'] = nreject + print(" rejected", nreject) + obslist.meta["nreject"] = nreject else: - obslist.meta['nreject']=0 + obslist.meta["nreject"] = 0 ncutout = len(imlist) - imflags = np.zeros(ncutout, dtype='i4') + imflags = np.zeros(ncutout, dtype="i4") - seglist = [m.interpolate_coadd_seg(iobj,i) - for i in range(1,ncutout+1)] + seglist = [ + m.interpolate_coadd_seg(iobj, i) for i in range(1, ncutout + 1) + ] # note starting at 1 assuming coadd is at 0 for i in range(ncutout): - icut=i+1 + icut = i + 1 - print(" cutout %d/%d" % (i+1,ncutout)) - im=imlist[i] - wt=wtlist[i] - seg=seglist[i] - bmask=bmlist[i] + print(" cutout %d/%d" % (i + 1, ncutout)) + im = imlist[i] + wt = wtlist[i] + seg = seglist[i] + bmask = bmlist[i] - file_id=m['file_id'][iobj, icut] - orig_row=m['orig_row'][iobj, icut] - orig_col=m['orig_col'][iobj, icut] + file_id = m["file_id"][iobj, icut] + orig_row = m["orig_row"][iobj, icut] + orig_col = m["orig_col"][iobj, icut] - if self['symmetrize_mask']: + if self["symmetrize_mask"]: self._symmetrize_weight(wt) self._symmetrize_bmask(bmask) if make_plots: - self._show_epoch(iobj,icut,im,seg,wt,bmask) + self._show_epoch(iobj, icut, im, seg, wt, bmask) if self._hits_the_edge(bmask): print(" cutout hits the edge") @@ -522,9 +514,11 @@ def _get_obslist(self, iobj): continue mask_frac = self._get_mask_frac(bmask, wt) - if mask_frac > self['max_masked_frac']: - print(" mask fraction %g exceeds " - "maximum %g" % (mask_frac, self['max_masked_frac'])) + if mask_frac > self["max_masked_frac"]: + print( + " mask fraction %g exceeds " + "maximum %g" % (mask_frac, self["max_masked_frac"]) + ) imflags[i] |= CUTOUT_MASK_FRAC continue @@ -539,7 +533,7 @@ def _get_obslist(self, iobj): continue wt_logic = wt > 0.0 - w=np.where(wt_logic) + w = np.where(wt_logic) if w[0].size == 0: print(" all weight is zero") imflags[i] |= CUTOUT_ALL_WEIGHT_ZERO @@ -547,40 +541,34 @@ def _get_obslist(self, iobj): jdict = m.get_jacobian(iobj, icut) jacobian = ngmix.Jacobian( - row=jdict['row0'], - col=jdict['col0'], - dudrow=jdict['dudrow'], - dudcol=jdict['dudcol'], - dvdrow=jdict['dvdrow'], - dvdcol=jdict['dvdcol'], + row=jdict["row0"], + col=jdict["col0"], + dudrow=jdict["dudrow"], + dudcol=jdict["dudcol"], + dvdrow=jdict["dvdrow"], + dvdcol=jdict["dvdcol"], ) - ccen = (np.array(im.shape)-1.0)/2.0 - offset_pixels=dict( - row_offset=jacobian.row0-ccen[0], - col_offset=jacobian.col0-ccen[1], + ccen = (np.array(im.shape) - 1.0) / 2.0 + offset_pixels = dict( + row_offset=jacobian.row0 - ccen[0], + col_offset=jacobian.col0 - ccen[1], ) - meta={'offset_pixels':offset_pixels} + meta = {"offset_pixels": offset_pixels} if False: import ngmix - fac=jacobian.get_scale()**2 - T = 4.0*fac - gmpsf = ngmix.GMixModel([0., 0., 0., 0.0, T, 1.0],"gauss") - gm0 = ngmix.GMixModel([0., 0., 0.025, 0.0, T, 1.0],"gauss") + + fac = jacobian.get_scale() ** 2 + T = 4.0 * fac + gmpsf = ngmix.GMixModel([0.0, 0.0, 0.0, 0.0, T, 1.0], "gauss") + gm0 = ngmix.GMixModel([0.0, 0.0, 0.025, 0.0, T, 1.0], "gauss") gm = gm0.convolve(gmpsf) im = gm.make_image(im.shape, jacobian=jacobian) noise = 0.0001 im += self.rng.normal(scale=noise, size=im.shape) - wt = wt*0 + 1.0/noise**2 - - - """ - print("image ccen:",ccen) - print("orig_row: %.3f orig_col: %.3f" % (orig_row, orig_col)) - print("stamp_row: %.3f stamp_col: %.3f" % (jdict['row0'],jdict['col0'])) - """ + wt = wt * 0 + 1.0 / noise ** 2 obs = ngmix.Observation( im, @@ -590,27 +578,32 @@ def _get_obslist(self, iobj): meta=meta, ) - obs.psf=self._get_psf_obs( - obs, file_id, meta, orig_row, orig_col, + obs.psf = self._get_psf_obs( + obs, + file_id, + meta, + orig_row, + orig_col, ) - obs.seg=seg + obs.seg = seg # interpolate im, weight, and a noise map # .noise attribute is added - obs.meta['ninterp']=self._interp_bad_pixels(obs) + obs.meta["ninterp"] = self._interp_bad_pixels(obs) - if self.make_plots and obs.meta['ninterp'] > 0: - self._show_epoch(iobj,icut,obs.image,seg, - obs.weight,obs.bmask, - extra='interp') + if self.make_plots and obs.meta["ninterp"] > 0: + self._show_epoch( + iobj, icut, obs.image, seg, + obs.weight, obs.bmask, extra="interp" + ) obslist.append(obs) return obslist, imflags def _hits_the_edge(self, bmask): if self.edge_flags is not None: - w=np.where( ((bmask & self.edge_flags) != 0) ) + w = np.where(((bmask & self.edge_flags) != 0)) hits_edge = w[0].size > 0 else: hits_edge = False @@ -618,7 +611,7 @@ def _hits_the_edge(self, bmask): return hits_edge def _get_mask_frac(self, bmask, wt): - logic=(wt <= 0.0) + logic = wt <= 0.0 if self.bad_flags is not None: logic = logic | ((bmask & self.bad_flags) != 0) @@ -633,8 +626,8 @@ def _get_mask_frac(self, bmask, wt): print(" %s %d/%d %g" % tup) """ - w=np.where(logic) - mask_frac = float(w[0].size)/bmask.size + w = np.where(logic) + mask_frac = float(w[0].size) / bmask.size return mask_frac def _symmetrize_weight(self, wt): @@ -643,7 +636,7 @@ def _symmetrize_weight(self, wt): """ assert wt.shape[0] == wt.shape[1] - wt_rot=np.rot90(wt) + wt_rot = np.rot90(wt) wzero = np.where(wt_rot == 0.0) if wzero[0].size > 0: @@ -660,26 +653,28 @@ def _symmetrize_bmask(self, bmask): bm_rot = np.rot90(bmask) - wbad = np.where( (bm_rot & self.bad_flags) != 0) + wbad = np.where((bm_rot & self.bad_flags) != 0) if wbad[0].size > 0: bmask[wbad] = bm_rot[wbad] def _central_region_is_masked(self, bmask, wt): - cconf=self['central_region'] - if cconf['check']: - rad=cconf['size']/2.0 + cconf = self["central_region"] + if cconf["check"]: + rad = cconf["size"] / 2.0 - cen = (np.array(wt.shape)-1.0)/2.0 - low=int(round(cen[0]-rad)) - high=int(round(cen[0]+rad+1)) + cen = (np.array(wt.shape) - 1.0) / 2.0 + low = int(round(cen[0] - rad)) + high = int(round(cen[0] + rad + 1)) - logic = (wt[low:high, low:high] > 0.0) + logic = wt[low:high, low:high] > 0.0 if self.bad_flags is not None: - logic = logic & \ - ((bmask[low:high, low:high] & self.bad_flags) == 0 ) + logic = ( + logic & + ((bmask[low:high, low:high] & self.bad_flags) == 0) + ) ok = np.all(logic) else: - ok=True + ok = True return not ok @@ -688,36 +683,28 @@ def _ps_edge_pixels_are_masked(self, bmask, wt): we can't interpolate near the edge, so cut extreme cases """ - conf=self['stamp_edge'] - if conf['check']: - bad_flags=self.bad_flags + conf = self["stamp_edge"] + if conf["check"]: + bad_flags = self.bad_flags sides_bad = ( - np.all(wt[0, :]==0.0) - or - np.all(wt[-1,:]==0.0) - or - np.all(wt[:, 0]==0.0) - or - np.all(wt[:,-1]==0.0) + np.all(wt[0, :] == 0.0) + or np.all(wt[-1, :] == 0.0) + or np.all(wt[:, 0] == 0.0) + or np.all(wt[:, -1] == 0.0) ) if not sides_bad: if bad_flags is not None: sides_bad = ( - np.all( (bmask[0, :] & bad_flags) != 0 ) - or - np.all( (bmask[-1,:] & bad_flags) != 0 ) - or - np.all( (bmask[:, 0] & bad_flags) != 0 ) - or - np.all( (bmask[:,-1] & bad_flags) != 0 ) + np.all((bmask[0, :] & bad_flags) != 0) + or np.all((bmask[-1, :] & bad_flags) != 0) + or np.all((bmask[:, 0] & bad_flags) != 0) + or np.all((bmask[:, -1] & bad_flags) != 0) ) - return sides_bad - def _get_psf_obs(self, obs, file_id, meta, row, col): """ psfex specific code here @@ -734,44 +721,46 @@ def _interp_bad_pixels(self, obs): add a noise image that is also interpolated, if needed """ - iconf=self['interp'] - assert iconf['type']=="cubic","only cubic interpolation for now" + iconf = self["interp"] + assert iconf["type"] == "cubic", "only cubic interpolation for now" - im=obs.image - weight=obs.weight + im = obs.image + weight = obs.weight if not obs.has_bmask(): - obs.bmask = np.zeros(im.shape, dtype='i4') + obs.bmask = np.zeros(im.shape, dtype="i4") bmask = obs.bmask bmravel = bmask.ravel() wtravel = weight.ravel() - bad_logic = (wtravel <= 0.0) + bad_logic = wtravel <= 0.0 if self.bad_flags is not None: - bad_logic = bad_logic | ( (bmravel & self.bad_flags) != 0 ) + bad_logic = bad_logic | ((bmravel & self.bad_flags) != 0) - wbad,=np.where(bad_logic) + (wbad,) = np.where(bad_logic) if wbad.size > 0: - print(" interpolating %d/%d masked or zero weight " - "pixels" % (wbad.size,im.size)) + print( + " interpolating %d/%d masked or zero weight " + "pixels" % (wbad.size, im.size) + ) - yy, xx = np.mgrid[0:im.shape[0], 0:im.shape[1]] + yy, xx = np.mgrid[0: im.shape[0], 0: im.shape[1]] x = xx.ravel() y = yy.ravel() - yx = np.zeros( (x.size, 2) ) - yx[:,0] = y - yx[:,1] = x + yx = np.zeros((x.size, 2)) + yx[:, 0] = y + yx[:, 1] = x - wgood, = np.where(bad_logic == False) + (wgood,) = np.where(~bad_logic) # make this a config option? - medwt=np.median(wtravel[wgood]) + medwt = np.median(wtravel[wgood]) wtravel[wbad] = medwt im_interp = self._do_interp(yx, im, wgood, wbad) @@ -781,12 +770,12 @@ def _interp_bad_pixels(self, obs): noise_interp = self._do_interp(yx, tmp_noise, wgood, wbad) - obs.image = im_interp - obs.noise = noise_interp + obs.image = im_interp + obs.noise = noise_interp else: obs.noise = self._make_noise_image(weight) - + return wbad.size def _make_noise_image(self, weight): @@ -794,20 +783,19 @@ def _make_noise_image(self, weight): create a noise image based on the input weight map """ - err_image = np.sqrt( 1.0/weight ) + err_image = np.sqrt(1.0 / weight) noise_image = self.rng.normal(size=weight.shape) noise_image *= err_image return noise_image - def _do_interp(self, yx, im, wgood, wbad): import scipy.interpolate im_ravel = im.ravel() ii = scipy.interpolate.CloughTocher2DInterpolator( - yx[wgood,:], + yx[wgood, :], im_ravel[wgood], fill_value=0.0, ) @@ -815,109 +803,111 @@ def _do_interp(self, yx, im, wgood, wbad): im_interp = im.copy() im_interp_ravel = im_interp.ravel() - vals = ii(yx[wbad,:]) + vals = ii(yx[wbad, :]) im_interp_ravel[wbad] = vals return im_interp - def _show_epoch(self, iobj, icut, im, seg, wt, bmask, extra=None): import biggles import images - tab=biggles.Table(2,2,aspect_ratio=1.0) + tab = biggles.Table(2, 2, aspect_ratio=1.0) + + tab[0, 0] = images.view(im, title="im", show=False) + tab[0, 1] = images.view(seg, title="seg", show=False) + tab[1, 0] = images.view(wt, title="weight", show=False) + tab[1, 1] = images.view(bmask, title="bmask", show=False) - tab[0,0] = images.view(im, title='im', show=False) - tab[0,1] = images.view(seg, title='seg', show=False) - tab[1,0] = images.view(wt, title='weight', show=False) - tab[1,1] = images.view(bmask, title='bmask', show=False) - - d='plots' + d = "plots" if not os.path.exists(d): try: os.makedirs(d) - except: + except FileExistsError: pass - fname=[ - 'object', - '%06d' % iobj, - '%02d' % icut, + fname = [ + "object", + "%06d" % iobj, + "%02d" % icut, ] if extra is not None: fname += [extra] - fname = '-'.join(fname) - fname += '.png' - fname=os.path.join(d, fname) + fname = "-".join(fname) + fname += ".png" + fname = os.path.join(d, fname) - #tab.show() - print(" writing:",fname) - tab.write_img(1500,1500,fname) + # tab.show() + print(" writing:", fname) + tab.write_img(1500, 1500, fname) - def _show_coadd(self, iobj, obs, original_coadd, original_wt, original_seg): + def _show_coadd(self, iobj, obs, original_coadd, + original_wt, original_seg): import biggles import images - tab=biggles.Table(3,3,aspect_ratio=1.0) + tab = biggles.Table(3, 3, aspect_ratio=1.0) - nl=None - tab[0,0] = images.view( - np.rot90( obs.image.transpose(),k=2), - #obs.image, + nl = None + tab[0, 0] = images.view( + np.rot90(obs.image.transpose(), k=2), + # obs.image, nonlinear=nl, - title='coadd (trans/rot)', + title="coadd (trans/rot)", show=False, ) - tab[0,1] = images.view(original_coadd, nonlinear=nl, title='orig coadd', show=False) + tab[0, 1] = images.view( + original_coadd, nonlinear=nl, title="orig coadd", show=False + ) - diff = np.rot90( obs.image.transpose(),k=2) - original_coadd - tab[0,2] = images.view(diff, title='diff', show=False) + diff = np.rot90(obs.image.transpose(), k=2) - original_coadd + tab[0, 2] = images.view(diff, title="diff", show=False) - cmin,cmax = obs.weight.min(), obs.weight.max() - tab[1,0] = images.view( - np.rot90( obs.weight.transpose(),k=2), - #obs.weight, + cmin, cmax = obs.weight.min(), obs.weight.max() + tab[1, 0] = images.view( + np.rot90(obs.weight.transpose(), k=2), + # obs.weight, nonlinear=nl, - title='weight (trans/rot) minmax: %.3g/%.3g' % (cmin,cmax), + title="weight (trans/rot) minmax: %.3g/%.3g" % (cmin, cmax), show=False, ) - cmin,cmax = original_wt.min(), original_wt.max() - tab[1,1] = images.view(original_wt, - nonlinear=nl, - title='orig weight minmax: %.3g/%.3g' % (cmin,cmax), - show=False) - - tab[1,2] = images.view(obs.noise, title='noise', show=False) - - #tab[0,2] = images.view(obs.seg, title='seg', show=False) - tab[2,0] = images.view( - np.rot90( obs.seg.transpose(),k=2), - #obs.seg, - title='seg (trans/rot)', + cmin, cmax = original_wt.min(), original_wt.max() + tab[1, 1] = images.view( + original_wt, + nonlinear=nl, + title="orig weight minmax: %.3g/%.3g" % (cmin, cmax), show=False, ) - tab[2,1] = images.view(original_seg, title='orig seg', show=False) - tab[2,2] = images.view(obs.psf.image, title='coadded psf', show=False) - - d='plots' + tab[1, 2] = images.view(obs.noise, title="noise", show=False) + + # tab[0,2] = images.view(obs.seg, title='seg', show=False) + tab[2, 0] = images.view( + np.rot90(obs.seg.transpose(), k=2), + # obs.seg, + title="seg (trans/rot)", + show=False, + ) + tab[2, 1] = images.view(original_seg, title="orig seg", show=False) + + tab[2, 2] = images.view(obs.psf.image, title="coadded psf", show=False) + + d = "plots" if not os.path.exists(d): try: os.makedirs(d) - except: + except FileExistsError: pass - fname=[ - 'object', - '%06d' % iobj, - 'coadd', + fname = [ + "object", + "%06d" % iobj, + "coadd", ] - fname = '-'.join(fname) - fname += '.png' - fname=os.path.join(d, fname) - - #tab.show() - print(" writing:",fname) - tab.write_img(1500,1500,fname) - + fname = "-".join(fname) + fname += ".png" + fname = os.path.join(d, fname) + # tab.show() + print(" writing:", fname) + tab.write_img(1500, 1500, fname) diff --git a/meds/compare.py b/meds/compare.py index 0cb7bf6..2618c1a 100644 --- a/meds/compare.py +++ b/meds/compare.py @@ -7,20 +7,23 @@ import numpy from .meds import MEDS + class Comparator(object): - def __init__(self, - file1, - file2, - id_name='id', - image_rms_tol=0.8, - same_nepoch=True, - png_prefix=False): + def __init__( + self, + file1, + file2, + id_name="id", + image_rms_tol=0.8, + same_nepoch=True, + png_prefix=False, + ): self.m1 = MEDS(file1) self.m2 = MEDS(file2) self.id_name = id_name - self.require_same_nepoch=same_nepoch - self.png_prefix=png_prefix + self.require_same_nepoch = same_nepoch + self.png_prefix = png_prefix self.image_rms_tol = image_rms_tol self.weight_tol = 1.0e-4 @@ -33,8 +36,8 @@ def go(self, nrand=None): self._match() self.compare_object_data() - for type in ['image','weight','bmask','seg']: - print("checking",type) + for type in ["image", "weight", "bmask", "seg"]: + print("checking", type) self.compare_images(type, nrand=nrand) def compare_object_data(self): @@ -57,12 +60,15 @@ def compare_object_data(self): ) if not check: import esutil as eu - #raise ValueError("field %s does not match" % n) + + # raise ValueError("field %s does not match" % n) diff = (c1[n][self.ind1] - c2[n][self.ind2]).ravel() rms = diff.std() - mn_clip, rms_clip = eu.stat.sigma_clip(diff,silent=True) - print("field %s does not match. rms %g " - "clipped: %g" % (n, rms,rms_clip)) + mn_clip, rms_clip = eu.stat.sigma_clip(diff, silent=True) + print( + "field %s does not match. rms %g " + "clipped: %g" % (n, rms, rms_clip) + ) except ValueError as err: print(err) @@ -85,85 +91,98 @@ def compare_images(self, type, nrand=None): else: ids = numpy.arange(self.ind1.size) - means=[] + means = [] ntot = ids.size for icount, index in enumerate(ids): - print(" %d/%d %d" % (icount+1,ntot,index)) + print(" %d/%d %d" % (icount + 1, ntot, index)) i1 = self.ind1[index] i2 = self.ind2[index] - ncut=self._check_ncutout(i1, i2) + ncut = self._check_ncutout(i1, i2) for icut in range(ncut): im2 = self.m2.get_cutout(i2, icut, type=type) im1 = self.m1.get_cutout(i1, icut, type=type) - #if type=='image': + # if type=='image': # wt=self.m2.get_cutout(i2, icut, type='weight') # err=numpy.sqrt(1.0/wt.max()) # print(" compare to noise:",err) - if type=='image': + if type == "image": ok, mean, err, std = self._compare_images(im1, im2) means.append(mean) - elif type == 'weight': + elif type == "weight": self._compare_weights(im1, im2) else: self._compare_images_exact(im1, im2) - if type=='image' and self.png_prefix is not None: + if type == "image" and self.png_prefix is not None: import images - fname=self.png_prefix + '-imdiff-%06d-%03d.png' % (index, icut) - print("writing:",fname) + + fname = ( + self.png_prefix + "-imdiff-%06d-%03d.png" % (index, icut) + ) + print("writing:", fname) if ok: - oks='OK' + oks = "OK" else: - oks='Not OK' + oks = "Not OK" + title = ( + "rms: %.3g mean: %.3g +/- %.3g %s" % (std, mean, err, oks) + ) images.compare_images( im1, im2, - dims=(1500,1500), - #width=1500, - #height=1500, - title='rms: %.3g mean: %.3g +/- %.3g %s' % (std, mean, err, oks), + dims=(1500, 1500), + # width=1500, + # height=1500, + title=title, file=fname, ) - #if 'q'==raw_input("hit a key: (q to quit)"): - # stop - if type=='image': + if type == "image": import esutil as eu + print(len(means)) means = numpy.array(means) mdiff = means.mean() - mdiff_err = means.std()/numpy.sqrt(means.size) + mdiff_err = means.std() / numpy.sqrt(means.size) print("average mean diff: %g +/- %g" % (mdiff, mdiff_err)) - mdiff, mdiff_std, mdiff_err = eu.stat.sigma_clip(means, get_err=True) + mdiff, mdiff_std, mdiff_err = eu.stat.sigma_clip( + means, get_err=True, + ) print("clipped average mean diff: %g +/- %g" % (mdiff, mdiff_err)) if self.png_prefix is not None: import biggles - binsize=mdiff_std*0.1 - fname=self.png_prefix + '-imdiff-means.png' + + binsize = mdiff_std * 0.1 + fname = self.png_prefix + "-imdiff-means.png" biggles.plot_hist( - means, binsize=binsize, visible=False, file=fname, + means, + binsize=binsize, + visible=False, + file=fname, ) def _check_ncutout(self, i1, i2): """ check the number of cutouts agrees """ - n1 = self.m1['ncutout'][i1] - n2 = self.m2['ncutout'][i2] + n1 = self.m1["ncutout"][i1] + n2 = self.m2["ncutout"][i2] if not self.require_same_nepoch: n = min(n1, n2) else: if n1 != n2: - raise ValueError("ncutout disagrees for objects " - "%d:%d %d:%d" % (i1,n1,i2,n2)) - n=n1 + raise ValueError( + "ncutout disagrees for objects " + "%d:%d %d:%d" % (i1, n1, i2, n2) + ) + n = n1 return n @@ -177,16 +196,18 @@ def _compare_images(self, im1, im2): raise ValueError("shapes do not " "match: %s %s" % (im1.shape, im2.shape)) - diff = im1-im2 + diff = im1 - im2 mean = diff.mean() std = diff.std() - err = std/numpy.sqrt( diff.size ) + err = std / numpy.sqrt(diff.size) if std > self.image_rms_tol: - print(" rms %g greater than " - "tolerance %g" % (std,self.image_rms_tol)) - ok=False + print( + " rms %g greater than " + "tolerance %g" % (std, self.image_rms_tol) + ) + ok = False else: - ok=True + ok = True return ok, mean, err, std @@ -198,13 +219,14 @@ def _compare_weights(self, im1, im2): raise ValueError("shapes do not " "match: %s %s" % (im1.shape, im2.shape)) - diff = im1-im2 + diff = im1 - im2 check = numpy.abs(diff) < self.weight_tol if not numpy.all(check): - import images - wbad = numpy.where(check == False) - print(" %d matched worse than %g. rms is %g" % \ - (wbad[0].size, self.weight_tol, diff.std())) + wbad = numpy.where(~check) + print( + " %d matched worse than %g. rms is %g" + % (wbad[0].size, self.weight_tol, diff.std()) + ) def _compare_images_exact(self, im1, im2): """ @@ -214,7 +236,7 @@ def _compare_images_exact(self, im1, im2): raise ValueError("shapes do not " "match: %s %s" % (im1.shape, im2.shape)) - wbad=numpy.where(im1 != im2) + wbad = numpy.where(im1 != im2) if wbad[0].size > 0: raise ValueError("%d/%d pixels did " @@ -225,6 +247,7 @@ def _get_random_subset(self, nrand): get random unique subset of range [0,n) """ import esutil as eu + return eu.random.random_indices(self.ind1.size, nrand) def _match(self): @@ -232,6 +255,7 @@ def _match(self): match by id """ import esutil as eu + self.ind1, self.ind2 = eu.numpy_util.match( self.m1[self.id_name], self.m2[self.id_name], diff --git a/meds/defaults.py b/meds/defaults.py index 135a6c4..f7ab7e6 100644 --- a/meds/defaults.py +++ b/meds/defaults.py @@ -1,15 +1,15 @@ -__version__="0.9.12" +__version__ = "0.9.14" -BMASK_EDGE=2**30 +BMASK_EDGE = 2**30 # default values in each image type. If the # cutout crosses an edge the value will get # filled in with this number default_values = { - 'image':0.0, - 'weight':0.0, - 'seg':0, - 'bmask':BMASK_EDGE, + 'image': 0.0, + 'weight': 0.0, + 'seg': 0, + 'bmask': BMASK_EDGE, } default_config = { @@ -18,7 +18,7 @@ 'bounds_buffer_uv': 16.0, # buffer for "coadd" (first entry) image in the meds in pixels - # objects within the coadd image bounds plus this buffer are + # objects within the coadd image bounds plus this buffer are # "in" the coadd image 'coadd_bounds_buffer_rowcol': 1e-3, @@ -28,9 +28,9 @@ # If True, make sure objects all are within the coadd bounds 'check_in_coadd': True, - # objects within 'coadd_bounds_buffer_rowcol' buffer around - # the bounds of the coadd are forced to be at the edge of - # the chip if 'force_into_coadd_bounds' is True + # objects within 'coadd_bounds_buffer_rowcol' buffer around + # the bounds of the coadd are forced to be at the edge of + # the chip if 'force_into_coadd_bounds' is True 'force_into_coadd_bounds': True, # allowed values in the bitmask image @@ -41,11 +41,10 @@ 'cutout_types': [], # default output data types for images - 'image_dtype':'f4', - 'weight_dtype':'f4', - 'seg_dtype':'i4', - 'bmask_dtype':'i4', - 'ormask_dtype':'i4', - 'noise_dtype':'f4', + 'image_dtype': 'f4', + 'weight_dtype': 'f4', + 'seg_dtype': 'i4', + 'bmask_dtype': 'i4', + 'ormask_dtype': 'i4', + 'noise_dtype': 'f4', } - diff --git a/meds/extractor.py b/meds/extractor.py index b7eb0a5..b466df8 100644 --- a/meds/extractor.py +++ b/meds/extractor.py @@ -12,16 +12,18 @@ import fitsio import numpy + def extract_range(meds_file, start, end, sub_file): """ Extract a subset of objects and write a new meds file. - If you want this as a temporary file, which will be cleaned + If you want this as a temporary file, which will be cleaned when you are done with it, use a MEDSExtractor object with cleanup=True """ - extractor=MEDSExtractor(meds_file, start, end, sub_file) + MEDSExtractor(meds_file, start, end, sub_file) + def extract_catalog(meds_file, sub_file): """ @@ -32,7 +34,7 @@ def extract_catalog(meds_file, sub_file): done with it, use a MEDSCatalogExtractor object with cleanup=True """ - extractor=MEDSCatalogExtractor(meds_file, sub_file) + MEDSCatalogExtractor(meds_file, sub_file) class MEDSExtractor(object): @@ -41,27 +43,28 @@ class MEDSExtractor(object): Optionally clean up the new file when the object is destroyed. """ + def __init__(self, meds_file, start, end, sub_file, cleanup=False, copy_all=False): - self.meds_file=meds_file - self.start=start - self.end=end - self.sub_file=sub_file - self.cleanup=cleanup - self.copy_all=copy_all + self.meds_file = meds_file + self.start = start + self.end = end + self.sub_file = sub_file + self.cleanup = cleanup + self.copy_all = copy_all # keep track of these in case copy_all is sent, so we # know what extra extensions to copy - self.default_ext=[ - 'object_data', - 'image_info', - 'metadata', - 'image_cutouts', - 'weight_cutouts', - 'seg_cutouts', - 'bmask_cutouts', - 'psf', - 'noise', + self.default_ext = [ + "object_data", + "image_info", + "metadata", + "image_cutouts", + "weight_cutouts", + "seg_cutouts", + "bmask_cutouts", + "psf", + "noise", ] self._check_inputs() self._extract() @@ -78,42 +81,42 @@ def __del__(self): def close(self): if self.cleanup: if os.path.exists(self.sub_file): - print('removing sub file:',self.sub_file) + print("removing sub file:", self.sub_file) os.remove(self.sub_file) def _extract(self): - + with fitsio.FITS(self.meds_file) as infits: - print('opening sub file:',self.sub_file) - with fitsio.FITS(self.sub_file,'rw',clobber=True) as outfits: + print("opening sub file:", self.sub_file) + with fitsio.FITS(self.sub_file, "rw", clobber=True) as outfits: # # subset of object data table # - obj_data = infits['object_data'][self.start:self.end+1] + obj_data = infits["object_data"][self.start: self.end + 1] cstart, cend = self._get_row_range(obj_data) if cstart != -1: # adjust to new start. If cstart==-1 will all be -9999 - obj_data['start_row'] -= cstart + obj_data["start_row"] -= cstart # psfs can have different dimensions - if 'psf' in infits: + if "psf" in infits: psf_cstart, psf_cend = get_psf_row_range(obj_data) if psf_cstart != -1: # adjust to new start. If cstart==-1 will all be -9999 - obj_data['psf_start_row'] -= psf_cstart + obj_data["psf_start_row"] -= psf_cstart - outfits.write(obj_data, extname='object_data') + outfits.write(obj_data, extname="object_data") # # copy all metadata and image info # - iinfo=infits['image_info'][:] - outfits.write(iinfo, extname='image_info') + iinfo = infits["image_info"][:] + outfits.write(iinfo, extname="image_info") - meta=infits['metadata'][:] - outfits.write(meta, extname='metadata') + meta = infits["metadata"][:] + outfits.write(meta, extname="metadata") # # extract cutouts for the requested objects @@ -121,90 +124,91 @@ def _extract(self): if cstart == -1: self._write_dummy(outfits) else: - image_cutouts=infits['image_cutouts'][cstart:cend] - outfits.write(image_cutouts, extname='image_cutouts') + image_cutouts = infits["image_cutouts"][cstart:cend] + outfits.write(image_cutouts, extname="image_cutouts") del image_cutouts - weight_cutouts=infits['weight_cutouts'][cstart:cend] - outfits.write(weight_cutouts, extname='weight_cutouts') + weight_cutouts = infits["weight_cutouts"][cstart:cend] + outfits.write(weight_cutouts, extname="weight_cutouts") del weight_cutouts - seg_cutouts=infits['seg_cutouts'][cstart:cend] - outfits.write(seg_cutouts, extname='seg_cutouts') + seg_cutouts = infits["seg_cutouts"][cstart:cend] + outfits.write(seg_cutouts, extname="seg_cutouts") del seg_cutouts - if 'bmask_cutouts' in infits: - bmask_cutouts=infits['bmask_cutouts'][cstart:cend] - outfits.write(bmask_cutouts, extname='bmask_cutouts') + if "bmask_cutouts" in infits: + bmask_cutouts = infits["bmask_cutouts"][cstart:cend] + outfits.write(bmask_cutouts, extname="bmask_cutouts") del bmask_cutouts - if 'noise_cutouts' in infits: - cutouts=infits['noise_cutouts'][cstart:cend] - outfits.write(cutouts, extname='noise_cutouts') + if "noise_cutouts" in infits: + cutouts = infits["noise_cutouts"][cstart:cend] + outfits.write(cutouts, extname="noise_cutouts") del cutouts - if 'psf' in infits: - psfs=infits['psf'][psf_cstart:psf_cend] - outfits.write(psfs, extname='psf') + if "psf" in infits: + psfs = infits["psf"][psf_cstart:psf_cend] + outfits.write(psfs, extname="psf") del psfs if self.copy_all: for hdu in infits: if hdu.has_data(): - extname=hdu.get_extname() + extname = hdu.get_extname() if extname not in self.default_ext: - h=hdu.read_header() - data=hdu.read(header=True) + h = hdu.read_header() + data = hdu.read(header=True) outfits.write(data, header=h, extname=extname) def _write_dummy(self, outfits): - print('no objects with cutouts, writing dummy data') - dummy=numpy.zeros(2, dtype='f4') + -9999 - outfits.write(dummy, extname='image_cutouts') - dummy=numpy.zeros(2, dtype='f4') - outfits.write(dummy, extname='weight_cutouts') - dummy=numpy.zeros(2, dtype='i4') + -9999 - outfits.write(dummy, extname='seg_cutouts') - dummy=numpy.zeros(2, dtype='i4') + -9999 - outfits.write(dummy, extname='bmask_cutouts') + print("no objects with cutouts, writing dummy data") + dummy = numpy.zeros(2, dtype="f4") + -9999 + outfits.write(dummy, extname="image_cutouts") + dummy = numpy.zeros(2, dtype="f4") + outfits.write(dummy, extname="weight_cutouts") + dummy = numpy.zeros(2, dtype="i4") + -9999 + outfits.write(dummy, extname="seg_cutouts") + dummy = numpy.zeros(2, dtype="i4") + -9999 + outfits.write(dummy, extname="bmask_cutouts") def _get_row_range(self, data): """ get pixel range for this subset """ - w,=numpy.where( data['ncutout'] > 0) - if w.size==0: + (w,) = numpy.where(data["ncutout"] > 0) + if w.size == 0: return -1, -1 - ifirst = w[0] - ilast = w[-1] + ilast = w[-1] - cstart = data['start_row'][ifirst,0] + cstart = data["start_row"][ifirst, 0] - ncutout_last = data['ncutout'][ilast] - npix_last = data['box_size'][ilast]**2 - cend = data['start_row'][ilast,ncutout_last-1] + npix_last + ncutout_last = data["ncutout"][ilast] + npix_last = data["box_size"][ilast] ** 2 + cend = data["start_row"][ilast, ncutout_last - 1] + npix_last return cstart, cend - def _check_inputs(self): - if self.meds_file==self.sub_file: + if self.meds_file == self.sub_file: raise ValueError("output file name equals input") if self.start > self.end: - raise ValueError("found start > end: %d %d" % (self.start,self.end) ) + raise ValueError( + "found start > end: %d %d" % (self.start, self.end) + ) class MEDSCatalogExtractor(object): """ Class to extract the catalog and meta data and write a new file """ + def __init__(self, meds_file, new_file, cleanup=False): - self.meds_file=meds_file - self.new_file=new_file - self.cleanup=cleanup + self.meds_file = meds_file + self.new_file = new_file + self.cleanup = cleanup self._check_inputs() self._extract() @@ -221,66 +225,63 @@ def __del__(self): def close(self): if self.cleanup: if os.path.exists(self.new_file): - print('removing cat only file:',self.new_file) + print("removing cat only file:", self.new_file) os.remove(self.new_file) def _extract(self): - + with fitsio.FITS(self.meds_file) as infits: - print('opening cat only file:',self.new_file) - with fitsio.FITS(self.new_file,'rw',clobber=True) as outfits: + print("opening cat only file:", self.new_file) + with fitsio.FITS(self.new_file, "rw", clobber=True) as outfits: # # object data table # - obj_data = infits['object_data'][:] + obj_data = infits["object_data"][:] - outfits.write(obj_data, extname='object_data') + outfits.write(obj_data, extname="object_data") # # copy all metadata and image info # - iinfo=infits['image_info'][:] - outfits.write(iinfo, extname='image_info') + iinfo = infits["image_info"][:] + outfits.write(iinfo, extname="image_info") - meta=infits['metadata'][:] - outfits.write(meta, extname='metadata') + meta = infits["metadata"][:] + outfits.write(meta, extname="metadata") def _check_inputs(self): - if self.meds_file==self.new_file: + if self.meds_file == self.new_file: raise ValueError("output file name equals input") + def get_psf_row_range(data): """ get pixel range for this subset """ - w,=numpy.where( data['ncutout'] > 0) - if w.size==0: + (w,) = numpy.where(data["ncutout"] > 0) + if w.size == 0: return -1, -1 - ifirst = w[0] - ilast = w[-1] + ilast = w[-1] - cstart = data['psf_start_row'][ifirst,0] + cstart = data["psf_start_row"][ifirst, 0] - ncutout_last = data['ncutout'][ilast] + ncutout_last = data["ncutout"][ilast] - if 'psf_box_size' in data.dtype.names: - if len(data['psf_box_size'].shape) > 1: + if "psf_box_size" in data.dtype.names: + if len(data["psf_box_size"].shape) > 1: # in this case the box size is allowed to be different for # different epochs - npix_last = data['psf_box_size'][ilast,ncutout_last-1]**2 + npix_last = data["psf_box_size"][ilast, ncutout_last - 1] ** 2 else: - npix_last = data['psf_box_size'][ilast]**2 + npix_last = data["psf_box_size"][ilast] ** 2 else: - rs = data['psf_row_size'][ilast,ncutout_last-1] - cs = data['psf_col_size'][ilast,ncutout_last-1] - npix_last = rs*cs + rs = data["psf_row_size"][ilast, ncutout_last - 1] + cs = data["psf_col_size"][ilast, ncutout_last - 1] + npix_last = rs * cs - cend = data['psf_start_row'][ilast,ncutout_last-1] + npix_last + cend = data["psf_start_row"][ilast, ncutout_last - 1] + npix_last return cstart, cend - - - diff --git a/meds/meds.py b/meds/meds.py index 2928451..d1c1476 100644 --- a/meds/meds.py +++ b/meds/meds.py @@ -2,11 +2,6 @@ import numpy import fitsio -try: - xrange # noqa: F821 -except Exception: - xrange = range - try: from . import _uberseg _have_c_ubserseg = True @@ -308,11 +303,10 @@ def get_psf(self, iobj, icutout): """ if not self.has_psf(): - raise RuntimeError("this MEDS file has no 'psf' extension") + raise ValueError("this MEDS file has no 'psf' extension") self._check_indices(iobj, icutout=icutout) - cat = self._cat shape = self._get_psf_shape(iobj, icutout) npix = shape[0]*shape[1] @@ -328,8 +322,8 @@ def _get_psf_shape(self, iobj, icutout): if 'psf_row_size' in cat.dtype.names: if len(cat['psf_row_size'].shape) > 1: shape = ( - cat['psf_row_size'][iobj,icutout], - cat['psf_col_size'][iobj,icutout], + cat['psf_row_size'][iobj, icutout], + cat['psf_col_size'][iobj, icutout], ) else: shape = ( @@ -360,7 +354,7 @@ def get_psf_list(self, iobj): A list of the PSFs as numpy arrays. """ ncut = self['ncutout'][iobj] - return [self.get_psf(iobj, icut) for icut in xrange(ncut)] + return [self.get_psf(iobj, icut) for icut in range(ncut)] def get_cweight_cutout(self, iobj, icutout, restrict_to_seg=False): """Composite the weight and seg maps, interpolating seg map from the @@ -532,7 +526,7 @@ def get_uberseg_list(self, iobj, fast=True): """ useg_list = [] - for i in xrange(self['ncutout'][iobj]): + for i in range(self['ncutout'][iobj]): uberseg = self.get_uberseg( iobj, i, @@ -920,7 +914,7 @@ def get_jacobian_list(self, iobj): """ self._check_indices(iobj) jlist = [] - for icutout in xrange(self['ncutout'][iobj]): + for icutout in range(self['ncutout'][iobj]): j = self.get_jacobian(iobj, icutout) jlist.append(j) @@ -1018,12 +1012,6 @@ def _make_composite_image( crow = crow.round().astype('i8') ccol = ccol.round().astype('i8') - ''' - wbad=numpy.where( (crow < 0) | (crow >= coadd_seg.shape[0]) - & (ccol < 0) | (ccol >= coadd_seg.shape[1]) ) - if wbad[0].size != 0: - cim[wbad] = 0 - ''' # clipping makes the notation easier crow = crow.clip(0, coadd_seg.shape[0]-1) ccol = ccol.clip(0, coadd_seg.shape[1]-1) @@ -1042,24 +1030,10 @@ def _make_composite_image( return cim def _get_extension_name(self, type): - if type == 'image': - return "image_cutouts" - elif type == "weight": - return "weight_cutouts" - elif type == "seg": - return "seg_cutouts" - elif type == "bmask": - return "bmask_cutouts" - elif type == "ormask": - return "ormask_cutouts" - elif type == "noise": - return "noise_cutouts" - else: - ext = "%s_cutouts" % type - if ext not in self._fits: - raise ValueError("bad cutout type '%s'" % type) - else: - return ext + ext = "%s_cutouts" % type + if ext not in self._fits: + raise ValueError("bad cutout type '%s'" % type) + return ext def _check_indices(self, iobj, icutout=None): if iobj >= self._cat.size: @@ -1112,7 +1086,7 @@ def split_mosaic(mosaic): ncutout = mosaic.shape[0]//box_size imlist = [] - for i in xrange(ncutout): + for i in range(ncutout): r1 = i*box_size r2 = (i+1)*box_size imlist.append(mosaic[r1:r2, :]) @@ -1178,7 +1152,7 @@ def reject_outliers(imlist, wtlist, nsigma=5.0, A=0.3): med = numpy.median(imstack, axis=0) - for i in xrange(nim): + for i in range(nim): im = imlist[i] wt = wtlist[i] diff --git a/meds/number_extractor.py b/meds/number_extractor.py index 742ee20..86886a1 100644 --- a/meds/number_extractor.py +++ b/meds/number_extractor.py @@ -1,4 +1,5 @@ -"""MEDSNumberExtractor +""" +MEDSNumberExtractor A class to extract a subset of the objects in a MEDS file and write to a new file using only the object numbers """ @@ -9,16 +10,18 @@ import numpy from .extractor import get_psf_row_range + def extract_numbers(meds_file, numbers, sub_file): """ Extract a subset of objects and write a new meds file. - If you want this as a temporary file, which will be cleaned + If you want this as a temporary file, which will be cleaned when you are done with it, use a MEDSNumberExtractor object with cleanup=True """ - extractor=MEDSNumbersExtractor(meds_file, numbers, sub_file) + MEDSNumberExtractor(meds_file, numbers, sub_file) + class MEDSNumberExtractor(object): """ @@ -26,19 +29,23 @@ class MEDSNumberExtractor(object): Optionally clean up the new file when the object is destroyed. """ + def __init__(self, meds_file, numbers, sub_file, cleanup=False): - self.meds_file=meds_file - nums = numpy.array(numbers,dtype=int) + self.meds_file = meds_file + nums = numpy.array(numbers, dtype=int) unums = numpy.unique(nums) - assert len(nums) == len(unums), "Input numbers must be unique! len(numbers) = %ld, len(unique(numbers)) = %ld" % (len(nums),len(unums)) + assert len(nums) == len(unums), ( + "Input numbers must be unique! len(numbers) = %ld, len(unique(numbers)) = %ld" # noqa + % (len(nums), len(unums)) + ) self.numbers = nums - self.sub_file=sub_file - self.cleanup=cleanup + self.sub_file = sub_file + self.cleanup = cleanup self._check_inputs() - + q = numpy.argsort(self.numbers) self.numbers = self.numbers[q] - + self._extract() def __enter__(self): @@ -53,110 +60,133 @@ def __del__(self): def close(self): if self.cleanup: if os.path.exists(self.sub_file): - print('removing sub file:',self.sub_file) + print("removing sub file:", self.sub_file) os.remove(self.sub_file) def _get_inds(self, data): inds = [] for number in self.numbers: - if number <= len(data['number']) and number >= 1 and number == data['number'][number-1]: - q = [number-1] + if ( + number <= len(data["number"]) + and number >= 1 + and number == data["number"][number - 1] + ): + q = [number - 1] else: - q, = numpy.where(number == data['number']) - assert len(q) == 1, "Could not find or found duplicate number: number = %ld, found %ld" % (number,len(q)) + (q,) = numpy.where(number == data["number"]) + assert ( + len(q) == 1 + ), "Could not find or found duplicate number: number = %ld, found %ld" % ( # noqa + number, + len(q), + ) inds.append(q[0]) - inds = numpy.array(inds,dtype=int) + inds = numpy.array(inds, dtype=int) return inds - + def _get_row_ranges(self, data): """ get pixel range for this subset """ - w,=numpy.where( data['ncutout'] > 0) - if w.size==0: + (w,) = numpy.where(data["ncutout"] > 0) + if w.size == 0: return [(-1, -1)] ranges = [] start = True - for w in xrange(len(data)): - if data['ncutout'][w] > 0: + for w in range(len(data)): + if data["ncutout"][w] > 0: if start: - ranges.append([data['start_row'][w,0], \ - data['start_row'][w,0]+data['box_size'][w]**2 * data['ncutout'][w]]) + ranges.append( + [ + data["start_row"][w, 0], + data["start_row"][w, 0] + + data["box_size"][w] ** 2 * data["ncutout"][w], + ] + ) start = False else: - if data['start_row'][w,0] != ranges[-1][-1]: - #add new if not back to back - ranges.append([data['start_row'][w,0], \ - data['start_row'][w,0]+data['box_size'][w]**2 * data['ncutout'][w]]) + if data["start_row"][w, 0] != ranges[-1][-1]: + # add new if not back to back + ranges.append( + [ + data["start_row"][w, 0], + data["start_row"][w, 0] + + data["box_size"][w] ** 2 * data["ncutout"][w], # noqa + ] + ) else: - #just expand range of pixels if back to back - ranges[-1][-1] = data['start_row'][w,0]+data['box_size'][w]**2 * data['ncutout'][w] + # just expand range of pixels if back to back + ranges[-1][-1] = ( + data["start_row"][w, 0] + + data["box_size"][w] ** 2 * data["ncutout"][w] + ) if len(ranges) == 0: - ranges.append([-1,-1]) - - #send back as list of tuples + ranges.append([-1, -1]) + + # send back as list of tuples tup_ranges = [tuple(rng) for rng in ranges] return tup_ranges - def _get_data_from_ranges(self, ranges, data): out_data = [] - dtype=None + dtype = None for rng in ranges: - if dtype == None: - arr = data[rng[0]:rng[1]] + if dtype is None: + arr = data[rng[0]: rng[1]] dtype = arr.dtype - out_data.extend(list(data[rng[0]:rng[1]])) - out_data = numpy.array(out_data,dtype=dtype) + out_data.extend(list(data[rng[0]: rng[1]])) + out_data = numpy.array(out_data, dtype=dtype) return out_data - + def _extract(self): - + with fitsio.FITS(self.meds_file) as infits: - print('opening sub file:',self.sub_file) - with fitsio.FITS(self.sub_file,'rw',clobber=True) as outfits: + print("opening sub file:", self.sub_file) + with fitsio.FITS(self.sub_file, "rw", clobber=True) as outfits: # # subset of object data table # - inds = self._get_inds(infits['object_data'][:]) - obj_data = infits['object_data'][inds] + inds = self._get_inds(infits["object_data"][:]) + obj_data = infits["object_data"][inds] ranges = self._get_row_ranges(obj_data) if ranges[0][0] != -1: - # adjust to new start. If ranges[0][0]==-1 will all be -9999 + # adjust to new start. + # If ranges[0][0]==-1 will all be -9999 loc = 0 - for i in xrange(len(obj_data)): - for j in xrange(obj_data['ncutout'][i]): - obj_data['start_row'][i,j] = loc - loc += obj_data['box_size'][i]**2 + for i in range(len(obj_data)): + for j in range(obj_data["ncutout"][i]): + obj_data["start_row"][i, j] = loc + loc += obj_data["box_size"][i] ** 2 - if 'psf' in infits: + if "psf" in infits: psf_ranges = get_psf_row_ranges(obj_data) if psf_ranges[0][0] != -1: - # adjust to new start. If ranges[0][0]==-1 will all be -9999 + # adjust to new start. + # If ranges[0][0]==-1 will all be -9999 loc = 0 - for iobj in xrange(len(obj_data)): - for icut in xrange(obj_data['ncutout'][iobj]): + for iobj in range(len(obj_data)): + for icut in range(obj_data["ncutout"][iobj]): + + obj_data["psf_start_row"][iobj, icut] = loc - obj_data['psf_start_row'][iobj,icut] = loc + rs = obj_data["psf_row_size"][iobj, icut] + cs = obj_data["psf_col_size"][iobj, icut] + loc += rs * cs - rs = obj_data['psf_row_size'][iobj,icut] - cs = obj_data['psf_col_size'][iobj,icut] - loc += rs*cs - - outfits.write(obj_data, extname='object_data') + outfits.write(obj_data, extname="object_data") # # copy all metadata and image info # - iinfo=infits['image_info'][:] - outfits.write(iinfo, extname='image_info') + iinfo = infits["image_info"][:] + outfits.write(iinfo, extname="image_info") - meta=infits['metadata'][:] - outfits.write(meta, extname='metadata') + meta = infits["metadata"][:] + outfits.write(meta, extname="metadata") # # extract cutouts for the requested objects @@ -164,65 +194,73 @@ def _extract(self): if ranges[0][0] == -1: self._write_dummy(outfits) else: - image_cutouts = self._get_data_from_ranges(ranges, infits['image_cutouts']) - outfits.write(image_cutouts, extname='image_cutouts') + image_cutouts = self._get_data_from_ranges( + ranges, infits["image_cutouts"] + ) + outfits.write(image_cutouts, extname="image_cutouts") del image_cutouts - - weight_cutouts = self._get_data_from_ranges(ranges, infits['weight_cutouts']) - outfits.write(weight_cutouts, extname='weight_cutouts') + + weight_cutouts = self._get_data_from_ranges( + ranges, infits["weight_cutouts"] + ) + outfits.write(weight_cutouts, extname="weight_cutouts") del weight_cutouts - seg_cutouts = self._get_data_from_ranges(ranges, infits['seg_cutouts']) - outfits.write(seg_cutouts, extname='seg_cutouts') + seg_cutouts = self._get_data_from_ranges( + ranges, infits["seg_cutouts"] + ) + outfits.write(seg_cutouts, extname="seg_cutouts") del seg_cutouts - bmask_cutouts = self._get_data_from_ranges(ranges, infits['bmask_cutouts']) - outfits.write(bmask_cutouts, extname='bmask_cutouts') + bmask_cutouts = self._get_data_from_ranges( + ranges, infits["bmask_cutouts"] + ) + outfits.write(bmask_cutouts, extname="bmask_cutouts") del bmask_cutouts - if 'psf' in infits: - psf_cutouts = self._get_data_from_ranges(psf_ranges, infits['psf']) - outfits.write(psf_cutouts, extname='psf') + if "psf" in infits: + psf_cutouts = self._get_data_from_ranges( + psf_ranges, infits["psf"] + ) + outfits.write(psf_cutouts, extname="psf") del psf_cutouts - def _write_dummy(self, outfits): - print('no objects with cutouts, writing dummy data') - dummy=numpy.zeros(2, dtype='f4') + -9999 - outfits.write(dummy, extname='image_cutouts') - dummy=numpy.zeros(2, dtype='f4') - outfits.write(dummy, extname='weight_cutouts') - dummy=numpy.zeros(2, dtype='i4') + -9999 - outfits.write(dummy, extname='seg_cutouts') - dummy=numpy.zeros(2, dtype='i4') + -9999 - outfits.write(dummy, extname='bmask_cutouts') + print("no objects with cutouts, writing dummy data") + dummy = numpy.zeros(2, dtype="f4") + -9999 + outfits.write(dummy, extname="image_cutouts") + dummy = numpy.zeros(2, dtype="f4") + outfits.write(dummy, extname="weight_cutouts") + dummy = numpy.zeros(2, dtype="i4") + -9999 + outfits.write(dummy, extname="seg_cutouts") + dummy = numpy.zeros(2, dtype="i4") + -9999 + outfits.write(dummy, extname="bmask_cutouts") def _check_inputs(self): - if self.meds_file==self.sub_file: + if self.meds_file == self.sub_file: raise ValueError("output file name equals input") if len(self.numbers) == 0: raise ValueError("one must extract at least one object") + def get_psf_row_ranges(data): """ get pixel range for this subset """ - w,=numpy.where( data['ncutout'] > 0) - if w.size==0: + (w,) = numpy.where(data["ncutout"] > 0) + if w.size == 0: return [(-1, -1)] ranges = [] - for w in xrange(len(data)): - if data['ncutout'][w] > 0: - rr = get_psf_row_range(data[w:w+1]) + for w in range(len(data)): + if data["ncutout"][w] > 0: + rr = get_psf_row_range(data[w: w + 1]) ranges.append(rr) if len(ranges) == 0: - ranges.append([-1,-1]) - - #send back as list of tuples + ranges.append([-1, -1]) + + # send back as list of tuples tup_ranges = [tuple(rng) for rng in ranges] return tup_ranges - - diff --git a/meds/tests/__init__.py b/meds/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/meds/tests/_fakemeds.py b/meds/tests/_fakemeds.py new file mode 100644 index 0000000..52e047b --- /dev/null +++ b/meds/tests/_fakemeds.py @@ -0,0 +1,367 @@ +import copy +import numpy as np +import fitsio +import meds + +DUDROW = 0.263 +DUDCOL = -0.01 +DVDROW = 0.01 +DVDCOL = 0.263 + +DET = DVDROW*DUDCOL - DVDCOL*DUDROW +PIXEL_SCALE = np.sqrt(abs(DET)) +PIXEL_AREA = PIXEL_SCALE**2 + + +CUTOUT_TYPES = [ + 'image', + 'weight', + 'seg', + 'bmask', + 'psf', +] + + +def make_fake_meds( + fname, + rng, + box_size=24, + ncutout_max=10, + nobj=10, + model='gauss', + flux=100.0, + fwhm=0.9, + noise=1.0, + with_psf=False, + psf_fwhm=0.9, + cutout_types=None, +): + + if cutout_types is None: + cutout_types = copy.deepcopy(CUTOUT_TYPES) + + with_psf = 'psf' in cutout_types + + print('writing:', fname) + with fitsio.FITS(fname, 'rw', clobber=True) as fits: + print('writing object_data') + object_data, total_pixels = make_object_data( + rng=rng, nobj=nobj, box_size=box_size, + ncutout_max=ncutout_max, + with_psf=with_psf, + ) + fits.write(object_data, extname='object_data') + + print('writing image_info') + image_info = make_image_info(nimage=ncutout_max) + fits.write(image_info, extname='image_info') + + print('writing metadata') + metadata = make_metadata() + fits.write(metadata, extname='metadata') + + print('reserving mosaic images') + reserve_mosaic_images( + fits=fits, cutout_types=cutout_types, total_pixels=total_pixels, + ) + + print('writing cutouts') + write_cutouts( + fits=fits, rng=rng, + object_data=object_data, box_size=box_size, + cutout_types=cutout_types, + fwhm=fwhm, + flux=flux, + noise=noise, + psf_fwhm=psf_fwhm, + ) + + +def write_cutouts( + fits, rng, object_data, box_size, cutout_types, + fwhm, + flux, + noise, + psf_fwhm, +): + for cutout_type in cutout_types: + + dtype = get_dtype(cutout_type) + extname = get_extname(cutout_type) + hdu = fits[extname] + + for iobj in range(object_data.size): + for icut in range(object_data['ncutout'][iobj]): + if cutout_type == 'image': + + row = object_data['cutout_row'][iobj, icut] + col = object_data['cutout_col'][iobj, icut] + + imdata = make_model_image( + row=row, col=col, + dims=[box_size]*2, + fwhm=fwhm, + flux=flux, + ) + imdata += rng.normal(scale=noise, size=imdata.shape) + + elif cutout_type == 'noise': + imdata = rng.normal(scale=noise, size=imdata.shape) + elif cutout_type == 'psf': + row, col = [(box_size - 1)/2]*2 + imdata = make_model_image( + row=row, col=col, + dims=[box_size]*2, + fwhm=psf_fwhm, + flux=1, + ) + + imdata += rng.normal(scale=1.0e-6, size=imdata.shape) + elif cutout_type == 'seg': + imdata = np.zeros((box_size, box_size), dtype=dtype) + imdata[:, :] = object_data['number'][iobj] + else: + imdata = np.zeros((box_size, box_size), dtype=dtype) + if cutout_type == 'weight': + imdata[:, :] = 1.0/noise**2 + + hdu.write(imdata, start=object_data['start_row'][iobj, icut]) + + +def make_model_image(row, col, dims, fwhm, flux): + + T = fwhm_to_T(fwhm) + irr = T/2 + icc = T/2 + irc = 0.0 + + gauss2d = np.zeros(1, dtype=_gauss2d_dtype) + gauss2d_set(gauss2d, flux, row, col, irr, irc, icc) + + return gauss2d_make_image(gauss2d, dims) + + +def get_extname(cutout_type): + if cutout_type == 'psf': + extname = cutout_type + else: + extname = '%s_cutouts' % cutout_type + return extname + + +def get_dtype(cutout_type): + if cutout_type in ['seg', 'bmask']: + dtype = 'i2' + else: + dtype = 'f4' + return dtype + + +def reserve_mosaic_images(fits, cutout_types, total_pixels): + for cutout_type in cutout_types: + + dtype = get_dtype(cutout_type) + extname = get_extname(cutout_type) + + fits.create_image_hdu( + img=None, + dtype=dtype, + dims=[total_pixels], + extname=extname, + ) + + +def make_object_data( + rng, + nobj, + box_size, + ncutout_max, + with_psf=False, +): + import meds + + extra_fields = [ + ('number', 'i4'), + ('flux_auto', 'f4'), + ('x2', 'f4'), + ('y2', 'f4'), + ] + if with_psf: + extra_fields += [ + ('psf_box_size', 'i4'), + ('psf_start_row', 'i4', ncutout_max), + ('psf_cutout_row', 'f4', ncutout_max), + ('psf_cutout_col', 'f4', ncutout_max), + ] + + data = meds.util.get_meds_output_struct( + nobj=nobj, ncutout_max=ncutout_max, + extra_fields=extra_fields, + ) + + for col in data.dtype.names: + data[col] = -9999 + + data['id'] = 1 + np.arange(nobj) + data['number'] = 1 + np.arange(nobj) + data['box_size'] = box_size + + data['ra'] = rng.uniform(low=0, high=3, size=nobj) + data['dec'] = rng.uniform(low=0, high=3, size=nobj) + data['ncutout'] = rng.randint(low=1, high=ncutout_max+1, size=nobj) + + # doesn't matter + data['orig_row'] = rng.uniform(0, 500, size=(nobj, ncutout_max)) + data['orig_col'] = rng.uniform(0, 500, size=(nobj, ncutout_max)) + data['orig_start_row'] = data['orig_row'] - box_size/2 + data['orig_start_col'] = data['orig_col'] - box_size/2 + + cen = (box_size - 1)/2 + data['cutout_row'] = cen + rng.uniform( + low=-0.5, high=0.5, size=(nobj, ncutout_max), + ) + data['cutout_col'] = cen + rng.uniform( + low=-0.5, high=0.5, size=(nobj, ncutout_max), + ) + data['dudrow'] = DUDROW + data['dudcol'] = DUDCOL + data['dvdrow'] = DVDROW + data['dvdcol'] = DVDCOL + + if with_psf: + data['psf_box_size'] = box_size + data['psf_cutout_row'] = cen + data['psf_cutout_col'] = cen + + ids = np.arange(ncutout_max) + total_pixels = 0 + current_row = 0 + npixels_per = box_size * box_size + for iobj in range(nobj): + ncutout = data['ncutout'][iobj] + data['file_id'][iobj, :ncutout] = ids[:ncutout] + + for icut in range(data['ncutout'][iobj]): + data['start_row'][iobj, icut] = current_row + + if with_psf: + data['psf_start_row'][iobj, icut] = current_row + + current_row += npixels_per + total_pixels += npixels_per + + return data, total_pixels + + +def make_image_info(nimage): + ii = meds.util.get_image_info_struct(nimage=nimage, path_len=10) + ii['scale'] = 1.0 + ii['image_path'] = 'blah.fits' + return ii + + +def make_metadata(): + metadata = np.zeros(1, dtype=[('medsconf', 'S10')]) + metadata['medsconf'] = 'meds01' + return metadata + + +def fwhm_to_T(fwhm): + """ + convert fwhm to T for a gaussian + """ + sigma = fwhm_to_sigma(fwhm) + return 2 * sigma ** 2 + + +def fwhm_to_sigma(fwhm): + """ + convert fwhm to sigma for a gaussian + """ + return fwhm / 2.3548200450309493 + + +def gauss2d_make_image(gauss2d, dims): + rows, cols = np.mgrid[ + 0:dims[0], + 0:dims[1], + ] + + rowdiff = rows - gauss2d['row'] + coldiff = cols - gauss2d['col'] + + v = DVDROW * rowdiff + DVDCOL * coldiff + u = DUDROW * rowdiff + DUDCOL * coldiff + + chi2 = ( + gauss2d["dcc"] * v * v + + gauss2d["drr"] * u * u + - 2.0 * gauss2d["drc"] * v * u + ) + + return gauss2d["pnorm"] * np.exp(-0.5 * chi2) * PIXEL_AREA + + +def gauss2d_set(gauss, p, row, col, irr, irc, icc): + """ + set the gaussian, clearing normalizations + """ + gauss["norm_set"] = 0 + + gauss["p"] = p + gauss["row"] = row + gauss["col"] = col + gauss["irr"] = irr + gauss["irc"] = irc + gauss["icc"] = icc + + gauss["det"] = irr * icc - irc * irc + + gauss2d_set_norm(gauss) + + +def gauss2d_set_norm(gauss): + """ + set the normalization, and nromalized variances + + A GMixRangeError is raised if the determinant is too small + + parameters + ---------- + gauss: a 2-d gaussian structure + See gmix.py + """ + + if gauss["det"] < 1.0e-200: + raise ValueError("det too low") + + T = gauss["irr"] + gauss["icc"] + if T <= 1.0e-200: + raise ValueError("T too low") + + idet = 1.0 / gauss["det"] + gauss["drr"] = gauss["irr"] * idet + gauss["drc"] = gauss["irc"] * idet + gauss["dcc"] = gauss["icc"] * idet + gauss["norm"] = 1.0 / (2 * np.pi * np.sqrt(gauss["det"])) + + gauss["pnorm"] = gauss["p"] * gauss["norm"] + + gauss["norm_set"] = 1 + + +_gauss2d_dtype = [ + ("p", "f8"), + ("row", "f8"), + ("col", "f8"), + ("irr", "f8"), + ("irc", "f8"), + ("icc", "f8"), + ("det", "f8"), + ("norm_set", "i8"), + ("drr", "f8"), + ("drc", "f8"), + ("dcc", "f8"), + ("norm", "f8"), + ("pnorm", "f8"), +] diff --git a/meds/tests/test_read.py b/meds/tests/test_read.py new file mode 100644 index 0000000..0f8dd8c --- /dev/null +++ b/meds/tests/test_read.py @@ -0,0 +1,190 @@ +import os +import pytest +import tempfile +import numpy as np +import meds +from ._fakemeds import make_fake_meds + + +@pytest.mark.parametrize('with_psf', [False, True]) +@pytest.mark.parametrize('with_bmask', [False, True]) +@pytest.mark.parametrize('with_ormask', [False, True]) +@pytest.mark.parametrize('with_noise', [False, True]) +def test_medsreaders_smoke(with_psf, with_bmask, with_ormask, with_noise): + + rng = np.random.RandomState(542) + + # assert meds.meds._have_c_ubserseg + + cutout_types = ['image', 'weight', 'seg'] + if with_bmask: + cutout_types += ['bmask'] + if with_ormask: + cutout_types += ['ormask'] + if with_noise: + cutout_types += ['noise'] + if with_psf: + cutout_types += ['psf'] + + with tempfile.TemporaryDirectory() as tdir: + fname = os.path.join(tdir, 'test-meds.fits') + make_fake_meds( + fname=fname, rng=rng, + cutout_types=cutout_types, + ) + + # test reading from single meds file + m = meds.MEDS(fname) + + image_info = m.get_image_info() + meta = m.get_meta() + assert 'medsconf' in meta.dtype.names + + for iobj in range(m.size): + box_size = m['box_size'][iobj] + ncutout = m['ncutout'][iobj] + + imlist = m.get_cutout_list(iobj, type='image') + mos = m.get_mosaic(iobj, type='image') + assert mos.shape == (ncutout * box_size, box_size) + + cseg_mos = m.interpolate_coadd_seg_mosaic(iobj) + assert cseg_mos.shape == (ncutout * box_size, box_size) + + cw_mos = m.get_cweight_mosaic(iobj) + assert cw_mos.shape == (ncutout * box_size, box_size) + + assert len(imlist) == ncutout, ('iobj: %s' % iobj) + + if with_psf: + assert m.has_psf() + psf_list = m.get_cutout_list(iobj, type='psf') + assert len(psf_list) == ncutout + + psf_list = m.get_psf_list(iobj) + assert len(psf_list) == ncutout + else: + with pytest.raises(ValueError): + m.get_cutout_list(iobj, type='psf') + + if with_bmask: + bmask_list = m.get_cutout_list(iobj, type='bmask') + assert len(bmask_list) == ncutout + else: + with pytest.raises(ValueError): + m.get_cutout_list(iobj, type='bmask') + + if with_ormask: + ormask_list = m.get_cutout_list(iobj, type='ormask') + assert len(ormask_list) == ncutout + else: + with pytest.raises(ValueError): + m.get_cutout_list(iobj, type='ormask') + + if with_noise: + noise_list = m.get_cutout_list(iobj, type='noise') + assert len(noise_list) == ncutout + else: + with pytest.raises(ValueError): + m.get_cutout_list(iobj, type='noise') + + for ctype in cutout_types: + for icut in range(ncutout): + cut = m.get_cutout(iobj, icut, type=ctype) + if ctype == 'psf': + psf_box_size = m['psf_box_size'][iobj] + assert cut.shape == (psf_box_size, ) * 2 + else: + assert cut.shape == (box_size, ) * 2 + + if ctype == 'seg' and icut == 0: + # we put this in, not required + assert np.any(cut == m['number'][iobj]) + + # uberseg + + for fast in [True, False]: + # hmm... fast c code is inaccessible for some reason + useglist = m.get_uberseg_list(iobj, fast=fast) + assert len(useglist) == ncutout + for icut in range(ncutout): + uberseg = m.get_uberseg(iobj, icut, fast=fast) + assert uberseg.shape == (box_size, ) * 2 + assert np.any(uberseg != 0) + + cseg_list = m.get_cseg_cutout_list(iobj) + assert len(cseg_list) == ncutout + + cwt_list = m.get_cweight_cutout_list(iobj) + assert len(cwt_list) == ncutout + + for icut in range(ncutout): + cseg_cutout = m.get_cseg_cutout(iobj, icut) + assert cseg_cutout.shape == (box_size, ) * 2 + + cs = m.interpolate_coadd_seg(iobj, icut) + assert cs.shape == (box_size, ) * 2 + + cw = m.get_cweight_cutout(iobj, icut) + assert cw.shape == (box_size, ) * 2 + + for use_canonical_cen in [True, False]: + cseg_weight = m.get_cseg_weight( + iobj, icut, use_canonical_cen=use_canonical_cen, + ) + assert cseg_weight.shape == (box_size, ) * 2 + + jlist = m.get_jacobian_list(iobj) + assert len(jlist) == ncutout + + for icut in range(ncutout): + source_info = m.get_source_info(iobj, icut) + ifile = m['file_id'][iobj, icut] + for f in source_info.dtype.names: + assert source_info[f] == image_info[f][ifile] + source_path = m.get_source_path(iobj, icut) + assert source_path == image_info['image_path'][ifile] + + j = m.get_jacobian(iobj, icut) + row0, col0 = m.get_cutout_rowcol(iobj, icut) + assert j['row0'] == row0 + assert j['col0'] == col0 + assert j['dudrow'] == m['dudrow'][iobj, icut] + assert j['dudcol'] == m['dudcol'][iobj, icut] + assert j['dvdrow'] == m['dvdrow'][iobj, icut] + assert j['dvdcol'] == m['dvdcol'][iobj, icut] + + jm = m.get_jacobian_matrix(iobj, icut) + assert jm[0, 0] == j['dudrow'] + assert jm[0, 1] == j['dudcol'] + assert jm[1, 0] == j['dvdrow'] + assert jm[1, 1] == j['dvdcol'] + + # test an error that can occur + with pytest.raises(ValueError): + m.get_cutout_list(0, type='blah') + with pytest.raises(ValueError): + m.get_cutout_list(1000, 0) + with pytest.raises(ValueError): + m.get_cutout_list(0, 1000) + + +def test_outlier_rejection(): + rng = np.random.RandomState(11) + with tempfile.TemporaryDirectory() as tdir: + fname = os.path.join(tdir, 'test-meds.fits') + make_fake_meds(fname=fname, rng=rng) + + m = meds.MEDS(fname) + imax = m['ncutout'].argmax() + + imlist = m.get_cutout_list(imax) + wtlist = m.get_cutout_list(imax, type='weight') + + dims = imlist[0].shape + rowbad, colbad = ((np.array(dims)-1)/2).astype('i4') + + imlist[0][rowbad, colbad] = 1.e9 + + meds.meds.reject_outliers(imlist, wtlist) + assert wtlist[0][rowbad, colbad] == 0.0 diff --git a/meds/util.py b/meds/util.py index 00a689c..5f31fc7 100644 --- a/meds/util.py +++ b/meds/util.py @@ -1,18 +1,18 @@ from __future__ import print_function -import os import numpy -import tempfile -import subprocess DEFVAL = -9999 -IMAGE_INFO_TYPES = ['image','weight','seg','bmask','bkg'] +IMAGE_INFO_TYPES = ["image", "weight", "seg", "bmask", "bkg"] + class MEDSCreationError(Exception): def __init__(self, value): self.value = value + def __str__(self): return repr(self.value) + def validate_meds(filename): """ validate the input MEDS file @@ -29,65 +29,65 @@ def validate_meds(filename): """ from .meds import MEDS - m=MEDS(filename) + m = MEDS(filename) - fits=m._fits + fits = m._fits print("checking for required extensions") - required_ext=['object_data','image_info','image_cutouts'] + required_ext = ["object_data", "image_info", "image_cutouts"] - nbad=0 + nbad = 0 for re in required_ext: - mess=" required extension named '%s' was not found" % re + mess = " required extension named '%s' was not found" % re if re not in fits: print(mess) nbad += 1 if nbad != 0: - print(" %d/%d were missing" % (nbad,len(required_ext))) + print(" %d/%d were missing" % (nbad, len(required_ext))) else: print(" OK") # require all fields print() print("checking for required object_data columns") - dt = numpy.dtype( get_meds_output_dtype(10) ) + dt = numpy.dtype(get_meds_output_dtype(10)) names = dt.names - onames = [c.lower() for c in fits['object_data'].get_colnames()] + onames = [c.lower() for c in fits["object_data"].get_colnames()] - nbad=0 + nbad = 0 for n in names: - mess=" required object_data field named '%s' was not found" % n + mess = " required object_data field named '%s' was not found" % n if n.lower() not in onames: print(mess) nbad += 1 if nbad != 0: - print(" %d/%d were missing" % (nbad,len(names))) + print(" %d/%d were missing" % (nbad, len(names))) else: print(" OK") # require only a subset print() print("checking for required image_info columns") - names=[ - 'image_path', - 'image_ext', - 'image_id', - 'image_flags', - 'magzp', - 'scale', + names = [ + "image_path", + "image_ext", + "image_id", + "image_flags", + "magzp", + "scale", ] - inames = [c.lower() for c in fits['image_info'].get_colnames()] + inames = [c.lower() for c in fits["image_info"].get_colnames()] - nbad=0 + nbad = 0 for n in names: - mess=" required image_info field named '%s' was not found" % n + mess = " required image_info field named '%s' was not found" % n if n not in inames: print(mess) nbad += 1 if nbad != 0: - print(" %d/%d were missing" % (nbad,len(names))) + print(" %d/%d were missing" % (nbad, len(names))) else: print(" OK") @@ -106,17 +106,18 @@ def get_meds_output_struct(nobj, ncutout_max, extra_fields=None): extra_fields: numpy descriptor A numpy type descriptor to add to the default set """ - dtype=get_meds_output_dtype(ncutout_max, extra_fields=extra_fields) + dtype = get_meds_output_dtype(ncutout_max, extra_fields=extra_fields) data = numpy.zeros(nobj, dtype=dtype) - noset=['ncutout'] + noset = ["ncutout"] for d in dtype: - name=d[0] + name = d[0] if name not in noset: data[name] = DEFVAL return data + def get_meds_input_struct(nobj, extra_fields=None): """ get the minimal object_data structure for input to the @@ -127,14 +128,14 @@ def get_meds_input_struct(nobj, extra_fields=None): nobj: int number of objects (length of array0 extra_fields: numpy descriptor, optional - optional extra fields to add + optional extra fields to add """ - dtype=get_meds_input_dtype(extra_fields=extra_fields) + dtype = get_meds_input_dtype(extra_fields=extra_fields) data = numpy.zeros(nobj, dtype=dtype) - noset=['ncutout'] + noset = ["ncutout"] for d in dtype: - name=d[0] + name = d[0] if name not in noset: data[name] = DEFVAL @@ -151,10 +152,10 @@ def get_meds_input_dtype(extra_fields=None): A numpy type descriptor to add to the default set """ dtype = [ - ('id', 'i8'), - ('box_size', 'i8'), - ('ra','f8'), - ('dec','f8'), + ("id", "i8"), + ("box_size", "i8"), + ("ra", "f8"), + ("dec", "f8"), ] if extra_fields is not None: @@ -163,7 +164,6 @@ def get_meds_input_dtype(extra_fields=None): return dtype - def get_meds_output_dtype(ncutout_max, extra_fields=None): """ Full dtype for the object_data structure and file extension @@ -177,23 +177,23 @@ def get_meds_output_dtype(ncutout_max, extra_fields=None): """ dtype = [ - ('id', 'i8'), - ('box_size', 'i8'), - ('ra','f8'), - ('dec','f8'), - ('ncutout', 'i8'), - ('file_id', 'i8', (ncutout_max,)), - ('start_row', 'i8', (ncutout_max,)), - ('orig_row', 'f8', (ncutout_max,)), - ('orig_col', 'f8', (ncutout_max,)), - ('orig_start_row', 'i8', (ncutout_max,)), - ('orig_start_col', 'i8', (ncutout_max,)), - ('cutout_row', 'f8', (ncutout_max,)), - ('cutout_col', 'f8', (ncutout_max,)), - ('dudrow', 'f8', (ncutout_max,)), - ('dudcol', 'f8', (ncutout_max,)), - ('dvdrow', 'f8', (ncutout_max,)), - ('dvdcol', 'f8', (ncutout_max,)), + ("id", "i8"), + ("box_size", "i8"), + ("ra", "f8"), + ("dec", "f8"), + ("ncutout", "i8"), + ("file_id", "i8", (ncutout_max,)), + ("start_row", "i8", (ncutout_max,)), + ("orig_row", "f8", (ncutout_max,)), + ("orig_col", "f8", (ncutout_max,)), + ("orig_start_row", "i8", (ncutout_max,)), + ("orig_start_col", "i8", (ncutout_max,)), + ("cutout_row", "f8", (ncutout_max,)), + ("cutout_col", "f8", (ncutout_max,)), + ("dudrow", "f8", (ncutout_max,)), + ("dudcol", "f8", (ncutout_max,)), + ("dvdrow", "f8", (ncutout_max,)), + ("dvdcol", "f8", (ncutout_max,)), ] if extra_fields is not None: @@ -202,11 +202,10 @@ def get_meds_output_dtype(ncutout_max, extra_fields=None): return dtype -def get_image_info_struct(nimage, path_len, - image_id_len=None, - wcs_len=None, - ext_len=None, - extra_dtype=None): +def get_image_info_struct( + nimage, path_len, image_id_len=None, + wcs_len=None, ext_len=None, extra_dtype=None +): """ get the image info structure @@ -236,15 +235,14 @@ def get_image_info_struct(nimage, path_len, data = numpy.zeros(nimage, dtype=dt) - data['scale'] = 1.0 + data["scale"] = 1.0 return data -def get_image_info_dtype(path_len, - image_id_len=None, - wcs_len=None, - ext_len=None, - extra_dtype=None): + +def get_image_info_dtype( + path_len, image_id_len=None, wcs_len=None, ext_len=None, extra_dtype=None +): """ get the image_info dtype for the specified path string length and wcs string length @@ -261,38 +259,38 @@ def get_image_info_dtype(path_len, instead of an integer, and this is the length """ - path_fmt = 'S%d' % path_len + path_fmt = "S%d" % path_len if image_id_len is None: - image_id_descr = 'i8' + image_id_descr = "i8" else: - image_id_descr = 'S%d' % image_id_len + image_id_descr = "S%d" % image_id_len if ext_len is not None: - ext_descr = 'S%d' % ext_len + ext_descr = "S%d" % ext_len else: - ext_descr = 'i2' - dt=[] + ext_descr = "i2" + dt = [] for ctype in IMAGE_INFO_TYPES: - path_name = '%s_path' % ctype - ext_name = '%s_ext' % ctype + path_name = "%s_path" % ctype + ext_name = "%s_ext" % ctype dt += [ (path_name, path_fmt), - (ext_name,ext_descr), + (ext_name, ext_descr), ] dt += [ - ('image_id', image_id_descr), - ('image_flags', 'i8'), - ('magzp', 'f4'), - ('scale', 'f4'), - ('position_offset','f8'), + ("image_id", image_id_descr), + ("image_flags", "i8"), + ("magzp", "f4"), + ("scale", "f4"), + ("position_offset", "f8"), ] if wcs_len is not None: - wcs_fmt = 'S%d' % wcs_len + wcs_fmt = "S%d" % wcs_len dt += [ - ('wcs',wcs_fmt), + ("wcs", wcs_fmt), ] if extra_dtype is not None: @@ -301,7 +299,6 @@ def get_image_info_dtype(path_len, return dt - def make_wcs_positions(row, col, offset, inverse=False): """ make a structure holding both the original wcs and zero-offset positions. @@ -321,29 +318,26 @@ def make_wcs_positions(row, col, offset, inverse=False): inverse: bool, optional Set to True if the input are zero based """ - n=row.size - dt=[('wcs_row','f8'), - ('wcs_col','f8'), - ('zrow','f8'), - ('zcol','f8')] + dt = [("wcs_row", "f8"), ("wcs_col", "f8"), ("zrow", "f8"), ("zcol", "f8")] - data=numpy.zeros(row.size, dtype=dt) + data = numpy.zeros(row.size, dtype=dt) if inverse: - data['zrow'] = row - data['zcol'] = col + data["zrow"] = row + data["zcol"] = col - data['wcs_row'] = row + offset - data['wcs_col'] = col + offset + data["wcs_row"] = row + offset + data["wcs_col"] = col + offset else: - data['wcs_row'] = row - data['wcs_col'] = col + data["wcs_row"] = row + data["wcs_col"] = col - data['zrow'] = row - offset - data['zcol'] = col - offset + data["zrow"] = row - offset + data["zcol"] = col - offset return data + # coordinates # ra = -u # ra = -phi @@ -354,30 +348,35 @@ def make_wcs_positions(row, col, offset, inverse=False): # u = -ra on sphere = +phi on sphere # v = dec on sphere = -theta on sphere -def radec_to_uv(ra,dec,ra_cen,dec_cen): - rhat_cen,uhat_cen,vhat_cen = radec_to_unitvecs_ruv(ra_cen,dec_cen) - rhat,uhat,vhat = radec_to_unitvecs_ruv(ra,dec) - cosang = numpy.dot(rhat,rhat_cen) - u = numpy.dot(rhat,uhat_cen)/cosang/numpy.pi*180.0*60.0*60.0 # arcsec - v = numpy.dot(rhat,vhat_cen)/cosang/numpy.pi*180.0*60.0*60.0 # arcsec - return u,v -def radec_to_unitvecs_ruv(ra,dec): - theta,phi = radec_to_thetaphi(ra,dec) - return thetaphi_to_unitvecs_ruv(theta,phi) +def radec_to_uv(ra, dec, ra_cen, dec_cen): + rhat_cen, uhat_cen, vhat_cen = radec_to_unitvecs_ruv(ra_cen, dec_cen) + rhat, uhat, vhat = radec_to_unitvecs_ruv(ra, dec) + cosang = numpy.dot(rhat, rhat_cen) + + fac = 1 / numpy.pi * 180.0 * 60.0 * 60.0 + u = numpy.dot(rhat, uhat_cen) / cosang * fac # arcsec + v = numpy.dot(rhat, vhat_cen) / cosang * fac # arcsec + return u, v + + +def radec_to_unitvecs_ruv(ra, dec): + theta, phi = radec_to_thetaphi(ra, dec) + return thetaphi_to_unitvecs_ruv(theta, phi) -def radec_to_thetaphi(ra,dec): - return (90.0-dec)/180.0*numpy.pi,-1.0*ra/180.0*numpy.pi -def thetaphi_to_unitvecs_ruv(theta,phi): +def radec_to_thetaphi(ra, dec): + return (90.0 - dec) / 180.0 * numpy.pi, -1.0 * ra / 180.0 * numpy.pi + + +def thetaphi_to_unitvecs_ruv(theta, phi): sint = numpy.sin(theta) cost = numpy.cos(theta) sinp = numpy.sin(phi) cosp = numpy.cos(phi) - rhat = numpy.array([sint*cosp,sint*sinp,cost]).T - that = numpy.array([cost*cosp,cost*sinp,-1.0*sint]).T - phat = numpy.array([-1.0*sinp,cosp,numpy.zeros_like(sinp)]).T - - return rhat,phat,-1.0*that + rhat = numpy.array([sint * cosp, sint * sinp, cost]).T + that = numpy.array([cost * cosp, cost * sinp, -1.0 * sint]).T + phat = numpy.array([-1.0 * sinp, cosp, numpy.zeros_like(sinp)]).T + return rhat, phat, -1.0 * that diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..33423df --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +numpy +fitsio +esutil +joblib diff --git a/setup.py b/setup.py index 232fbdc..0d681e1 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ setup( name="meds", - version="0.9.12", + version="0.9.14", description="Python and C libraries for reading MEDS files", license="GNU GPLv3", author="Erin Scott Sheldon",