Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
itune1 = 5
itune2 = 5
ntune = 0
IRLS_scale = 25.
bands=['r','g','u','i','z']
bandname = 'r'
flipBands = ['r']
rerun = 0
TI = []
sources = []
for run,field,camcol,roi in zip(runs,fields,camcols,rois):
TI.extend([st.get_tractor_image(run, camcol, field, bandname,roi=roi,useMags=True) for bandname in bands])
sources.append(st.get_tractor_sources(run, camcol, field,bandname,roi=roi,bands=bands))
timg,info = TI[0]
photocal = timg.getPhotoCal()
wcs = timg.getWcs()
lvl = logging.DEBUG
logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)
tims = [timg for timg,tinf in TI]
tractor = st.SDSSTractor(tims)
for source in sources:
tractor.addSources(source)
zr = np.array([-5.,+5.]) * info['skysig']
print bands
TI = []
for i,(run,camcol,field) in enumerate(RCFS):
for bandname in bands:
im,inf = st.get_tractor_image(run, camcol, field, bandname,
useMags=True,
roiradecsize=(ra,dec,S))
#im.dxdy = (fxc - xc, fyc - yc)
TI.append((im,inf))
tims = [im for im,inf in TI]
skyvals = [(info['skysig'], info['sky']) for timg,info in TI]
zrs = [np.array([-1.,+6.]) * std + sky for std,sky in skyvals]
# Grab sources from the LAST RCF
roi = inf['roi']
sources = st.get_tractor_sources(run, camcol, field, roi=roi,
bands=bands)
# DevGalaxy
s0 = sources[0]
s0.re = 0.01
s0.brightness.r = 22.
tractor = Tractor(tims)
tractor.addSources(sources)
lnp0 = tractor.getLogProb()
print 'Catalog:'
print tractor.catalog
print tractor.catalog.getParams()
ntune = 0
IRLS_scale = 25.
bands=['r','g','u','i','z']
bandname = 'r'
flipBands = ['r']
rerun = 0
TI = []
sources = []
for rcf in rcfs:
print rcf
try:
TI.extend([st.get_tractor_image(rcf[0], rcf[1], rcf[2], bandname,useMags=True) for bandname in bands])
sources.append(st.get_tractor_sources(rcf[0], rcf[1], rcf[2],bandname,bands=bands))
except:
print "Sources not in data! %s" % (rcf)
timg,info = TI[0]
photocal = timg.getPhotoCal()
wcs = timg.getWcs()
lvl = logging.DEBUG
logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)
tims = [timg for timg,tinf in TI]
tractor = st.SDSSTractor(tims)
for source in sources:
tractor.addSources(source)
zr = np.array([-5.,+5.]) * info['skysig']
if not os.path.exists(fn):
sdss.retrieve('tsField', run, camcol, field)
tsf = sdss.readTsField(run, camcol, field, rerun)
astrans = tsf.getAsTrans(canonband)
x,y = astrans.radec_to_pixel(ra, dec)
print 'x,y', x,y
roi = [np.clip(x-pixr, 0, W), np.clip(x+pixr, 0, W),
np.clip(y-pixr, 0, H), np.clip(y+pixr, 0, H)]
if roi[0] == roi[1] or roi[2] == roi[3]:
print 'Skipping', run, camcol, field
continue
TI.extend([st.get_tractor_image(run, camcol, field, band, useMags=True,
roiradecsize=(ra,dec,pixr), sdssobj=sdss)
for band in bands])
sources.extend(st.get_tractor_sources(run, camcol, field, bandname=canonband,
bands=bands, roi=roi))
# (could also get 'roi' from the st.get_tractor_image info dict)
tractor = Tractor([tim for tim,tinf in TI], sources, mp)
for im in tractor.getImages():
im.freezeParams('wcs', 'photocal', 'psf', 'sky')
#tractor.freezeParam('catalog')
for i in range(5):
if True or (i % 5 == 0):
print 'Thawing sky...'
tractor.images.thawParamsRecursive('sky')
roi = [x0,x1,y0,y1]
ra = 133.175
dec = 33.4167
itune = 10
ntune = 5
bands=['r','g','i','u','z']
bandname = 'r'
flipBands = ['r']
rerun = 0
TI = []
TI.extend([st.get_tractor_image(run, camcol, field, bandname,roi=roi,useMags=True) for bandname in bands])
sources = st.get_tractor_sources(run, camcol, field,bandname,roi=roi, bands=bands)
timg,info = TI[0]
photocal = timg.getPhotoCal()
wcs = timg.getWcs()
lvl = logging.DEBUG
logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)
tims = [timg for timg,tinf in TI]
tractor = st.SDSSTractor(tims)
tractor.addSources(sources)
zr = np.array([-5.,+5.]) * info['skysig']
print bands
wcs = im.getWcs()
# im,info = TI[bands.index('r')]
#wcs = im.getWcs()
# this is a shifted WCS; S,S is the center.
# rd = wcs.pixelToPosition(None, (S,S))
# ra,dec = rd.ra,rd.dec
# print 'RA,Dec', ra,dec
cd = wcs.cdAtPixel(xc,yc)
pixscale = np.sqrt(np.abs(np.linalg.det(cd)))
print('pixscale', pixscale * 3600.)
extent = pixscale * S
#for timg,info in TI:
# print timg.hashkey()
sources = st.get_tractor_sources(run, camcol, field, roi=roi,
bands=bands)
print('Sources:')
for s in sources:
print(s)
print()
print('Tractor images:', TI)
cftimg,cfsky,cfstd = get_cfht_img(ra,dec, extent)
# Create fake tractor.Image
# psf = NCircularGaussianPSF([0.1], [1.])
# sky = timg.sky
# print 'SDSS Sky:', timg.sky
# photocal = timg.photocal
x0,y0 = roi[0],roi[2]
# rd1 is near the position of the SDSS-catalog object
# rd2 is four pixels above.
rd1 = RaDecPos(225.67954, 11.265948)
rd2 = RaDecPos(225.67991, 11.265852)
lvl = logging.INFO
logging.basicConfig(level=lvl, format='%(message)s', stream=sys.stdout)
bandname = band
sdssprefix = '%06i-%s%i-%04i' % (run, bandname, camcol, field)
timg,info = st.get_tractor_image(run, camcol, field, bandname,
roi=roi)
sources = st.get_tractor_sources(run, camcol, field, bandname,
roi=roi)
wcs = timg.getWcs()
for source in sources:
x,y = wcs.positionToPixel(source, source.getPosition())
print(' (%.2f, %.2f):' % (x+x0,y+y0), source)
tractor = st.SDSSTractor([timg])
tractor.addSources(sources)
lnlorig = tractor.getLogLikelihood()
zr = np.array([-5.,+20.]) * info['skysig']
save(sdssprefix+'-orig', tractor, zr)
# find closest source
mini = -1