Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
plt.clf()
img = tim.getImage()
mn,mx = [np.percentile(img,p) for p in [25,99]]
dimshow(img, vmin=mn, vmax=mx)
xx,yy = [],[]
for src in srcs:
x,y = tim.getWcs().positionToPixel(src.getPosition())
xx.append(x)
yy.append(y)
ax = plt.axis()
plt.plot(xx, yy, 'r+')
plt.axis(ax)
plt.savefig('tim-%s%i.png' % (band, ifield))
# Populate the FITS header we will put in the output file.
hdr = fitsio.FITSHDR()
hdr.add_record(dict(name='RUN', value=run, comment='SDSS run'))
hdr.add_record(dict(name='CAMCOL', value=camcol, comment='SDSS camcol'))
hdr.add_record(dict(name='FIELD', value=field, comment='SDSS field'))
hdr.add_record(dict(name='BAND', value=band, comment='SDSS band'))
# Copy from input "frame" header
orighdr = tinfo['hdr']
for key in ['NMGY']:
hdr.add_record(dict(name=key, value=orighdr[key],
comment=orighdr.get_comment(key)))
hdr.add_record(dict(name='CALIB', value=tim.sdss_calib,
comment='Mean "calibvec" value for this image'))
hdr.add_record(dict(name='SKY', value=tim.sdss_sky,
comment='SDSS sky estimate at image center'))
hdr.add_record(dict(name='GAIN', value=tim.sdss_gain,
for (key, value, comment) in addToHeader:
hdr.add_record(dict(name=key, value=value, comment=comment))
# Convert tim's WCS to a local Tan.
cd = tim.wcs.cdAtPixel(x, y)
tanwcs = Tan(ra, dec, fx, fy,
cd[0,0], cd[0,1], cd[1,0], cd[1,1],
x1-x0, y1-y0)
tim.wcs = ConstantFitsWcs(tanwcs)
imghdr = fitsio.FITSHDR()
imghdr.add_record(dict(name='DATATYPE', value='IMAGE',
comment='SDSS image pixels in nanomaggies'))
tanwcs.add_to_header(imghdr)
ivhdr = fitsio.FITSHDR()
ivhdr.add_record(dict(name='DATATYPE', value='INVVAR',
comment='SDSS image pixel inverse variance in 1/nmgy^2'))
tanwcs.add_to_header(ivhdr)
fn = stamp_pattern % band
print('writing', fn)
F = fitsio.FITS(fn, mode='rw', clobber=True)
tim.toFits(F, primheader=hdr,
imageheader=imghdr, invvarheader=ivhdr)
continue
# Resample to common grid
th,tw = tim.shape
for iband,band in enumerate(bands):
i = allbands.index(band)
TT.decam_apflux[:,i,:] = AP.get('apflux_img_%s' % band)
TT.decam_apflux_ivar[:,i,:] = AP.get('apflux_img_ivar_%s' % band)
TT.decam_apflux_resid[:,i,:] = AP.get('apflux_resid_%s' % band)
cat.thawAllRecursive()
hdr = None
T2,hdr = prepare_fits_catalog(cat, invvars, TT, hdr, bands, fs,
allbands=allbands)
# mod
T2.ra += (T2.ra < 0) * 360.
T2.ra -= (T2.ra > 360) * 360.
primhdr = fitsio.FITSHDR()
for r in version_header.records():
primhdr.add_record(r)
primhdr.add_record(dict(name='ALLBANDS', value=allbands,
comment='Band order in array values'))
for i,ap in enumerate(apertures_arcsec):
primhdr.add_record(dict(name='APRAD%i' % i, value=ap,
comment='Aperture radius, in arcsec'))
bits = CP_DQ_BITS.values()
bits.sort()
bitmap = dict((v,k) for k,v in CP_DQ_BITS.items())
for i in range(16):
bit = 1<
hdrs = []
#s = ''
hdr = None
for line in f.readlines():
#print('Read %i chars' % len(line), ': ' + line[:25]+'...')
line = line.strip()
#print('Stripped to %i' % len(line))
#line = line + ' ' * (80 - len(line))
#assert(len(line) == 80)
#s += line
if line.startswith('END'):
#print('Found END')
hdr = None
continue
if hdr is None:
hdr = fitsio.FITSHDR()
hdrs.append(hdr)
#print('Started header number', len(hdrs))
hdr.add_record(line)
return hdrs
def write_fits(self, filename, hdr=None):
import fitsio
tt = type(self)
sky_type = '%s.%s' % (tt.__module__, tt.__name__)
if hdr is None:
hdr = fitsio.FITSHDR()
hdr.add_record(dict(name='SKY', value=sky_type,
comment='Sky class'))
self.toFitsHeader(hdr, prefix='SKY_')
fitsio.write(filename, None, header=hdr, clobber=True)
#version_hdr.add_record(dict(name='CPFILE', value=im.imgfn, comment='DECam comm.pipeline file'))
version_hdr.add_record(dict(name='CPFILE', value=fname, comment='DECam comm.pipeline file'))
version_hdr.add_record(dict(name='CPHDU', value=im.hdu, comment='DECam comm.pipeline ext'))
version_hdr.add_record(dict(name='CAMERA', value='DECam', comment='Dark Energy Camera'))
version_hdr.add_record(dict(name='EXPNUM', value=im.expnum, comment='DECam exposure num'))
version_hdr.add_record(dict(name='CCDNAME', value=im.extname, comment='DECam CCD name'))
version_hdr.add_record(dict(name='FILTER', value=tim.band, comment='Bandpass of this image'))
version_hdr.add_record(dict(name='EXPOSURE', value='decam-%s-%s' % (im.expnum, im.extname), comment='Name of this image'))
keys = ['TELESCOP','OBSERVAT','OBS-LAT','OBS-LONG','OBS-ELEV',
'INSTRUME']
for key in keys:
if key in tim.primhdr:
version_hdr.add_record(dict(name=key, value=tim.primhdr[key]))
hdr = fitsio.FITSHDR()
units = {'mjd':'sec', 'exptime':'sec', 'flux':'nanomaggy',
'flux_ivar':'1/nanomaggy^2'}
columns = F.get_columns()
for i,col in enumerate(columns):
if col in units:
hdr.add_record(dict(name='TUNIT%i' % (i+1), value=units[col]))
outdir = os.path.dirname(outfn)
if len(outdir):
trymakedirs(outdir)
fitsio.write(outfn, None, header=version_hdr, clobber=True)
F.writeto(outfn, header=hdr, append=True)
print 'Wrote', outfn
print 'Finished forced phot:', Time()-t0
plt.ylabel('Tractor W1 mag')
plt.title('WISE catalog vs Tractor forced photometry')
plt.axis([cathi,lo,cathi,lo])
ps.savefig()
print 'Tile', tile.coadd_id, 'band', wband, 'took', Time()-tb0
T.cut(inbounds)
T.delete_column('psfflux')
T.delete_column('cmodelflux')
T.delete_column('devflux')
T.delete_column('expflux')
T.treated_as_pointsource = T.treated_as_pointsource.astype(np.uint8)
hdr = fitsio.FITSHDR()
hdr.add_record(dict(name='SEQ_VER', value=version['Revision'],
comment='SVN revision'))
hdr.add_record(dict(name='SEQ_URL', value=version['URL'], comment='SVN URL'))
hdr.add_record(dict(name='SEQ_DATE', value=datetime.datetime.now().isoformat(),
comment='forced phot run time'))
hdr.add_record(dict(name='SEQ_NNEG', value=opt.nonneg, comment='non-negative?'))
hdr.add_record(dict(name='SEQ_SKY', value=opt.sky, comment='fit sky?'))
for band in bands:
minsig = getattr(opt, 'minsig%i' % band)
hdr.add_record(dict(name='SEQ_MNS%i' % band, value=minsig,
comment='min surf brightness in sig, band %i' % band))
hdr.add_record(dict(name='SEQ_BL', value=opt.blocks, comment='image blocks'))
hdr.add_record(dict(name='SEQ_CERE', value=opt.ceres, comment='use Ceres?'))
hdr.add_record(dict(name='SEQ_ERRF', value=opt.errfrac, comment='error flux fraction'))
if opt.ceres:
hdr.add_record(dict(name='SEQ_CEBL', value=opt.ceresblock,
def main(opt, cs82field):
t0 = Time()
bands = opt.bands
if opt.plots:
ps = PlotSequence(opt.prefix)
else:
ps = None
version = get_git_version()
print('git version info:', version)
hdr = fitsio.FITSHDR()
hdr.add_record(dict(name='PHO_VER', value=version['commit'],
comment='cs82.py code version (git commit)'))
if 'describe' in version:
hdr.add_record(dict(name='PHO_GIT', value=version['describe'],
comment='cs82.py code version (git describe)'))
hdr.add_record(dict(name='PHO_DATE', value=datetime.datetime.now().isoformat(),
comment='cs82.py photometry run time'))
hdr.add_record(dict(name='PHO_MACH', value=socket.getfqdn(),
comment='machine where phot was run'))
T,F = getTables(cs82field, enclosed=False)
masks = get_cs82_masks(cs82field)
sdss = DR9(basedir='data/unzip')
if opt.local:
# Running out of memory can cause failure to converge
# and term status = 2.
# Fail completely in this case.
if term != 0:
raise RuntimeError(
'Ceres terminated with status %i' % term)
if wantims:
ims0 = R.ims0
ims1 = R.ims1
IV, fs = R.IV, R.fitstats
if save_fits:
import fitsio
(dat, mod, ie, chi, roi) = ims1[0]
wcshdr = fitsio.FITSHDR()
tim.wcs.wcs.add_to_header(wcshdr)
tag = 'fit-%s-w%i' % (tile.coadd_id, band)
fitsio.write('%s-data.fits' %
tag, dat, clobber=True, header=wcshdr)
fitsio.write('%s-mod.fits' % tag, mod,
clobber=True, header=wcshdr)
fitsio.write('%s-chi.fits' % tag, chi,
clobber=True, header=wcshdr)
if get_models:
(dat, mod, ie, chi, roi) = ims1[0]
models[(tile.coadd_id, band)] = (mod, tim.roi)
if ps:
tag = '%s W%i' % (tile.coadd_id, band)
def getStandardFitsHeader(self, header=None):
if header is None:
import fitsio
hdr = fitsio.FITSHDR()
else:
hdr = header
psf = self.getPsf()
wcs = self.getWcs()
sky = self.getSky()
pcal = self.getPhotoCal()
psf.toStandardFitsHeader(hdr)
wcs.toStandardFitsHeader(hdr)
sky.toStandardFitsHeader(hdr)
pcal.toStandardFitsHeader(hdr)
return hdr