Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
S = klpsf.shape[0]
# number of Gaussian components
K = 3
w, mu, sig = em_init_params(K, None, None, None)
II = klpsf.copy()
II /= II.sum()
# HIDEOUS HACK
II = np.maximum(II, 0)
#print('Multi-Gaussian PSF fit...')
xm, ym = -(S // 2), -(S // 2)
if savepsfimg is not None:
plt.clf()
plt.imshow(II, interpolation='nearest', origin='lower')
plt.title('PSF image to fit with EM')
plt.savefig(savepsfimg)
res = em_fit_2d(II, xm, ym, w, mu, sig)
#print('em_fit_2d result:', res)
if res == 0:
# print('w,mu,sig', w,mu,sig)
mypsf = GaussianMixturePSF(w, mu, sig)
mypsf.computeRadius()
else:
# Failed! Return 'dg' model instead?
print('PSF model fit', psf, 'failed! Returning DG model instead')
psf = 'dg'
if psf == 'dg':
print('Creating double-Gaussian PSF approximation')
(a, s1, b, s2) = dgpsf
mypsf = NCircularGaussianPSF([s1, s2], [a, b])
if brightpsf:
print('Wrapping PSF in SdssBrightPSF')
del I
del var
psfimg = pyfits.open(psffn)[0].data
print('PSF image shape', psfimg.shape)
# number of Gaussian components
K = 3
PS = psfimg.shape[0]
w,mu,sig = em_init_params(K, None, None, None)
II = psfimg.copy()
II /= II.sum()
# HACK
II = np.maximum(II, 0)
print('Multi-Gaussian PSF fit...')
xm,ym = -(PS/2), -(PS/2)
em_fit_2d(II, xm, ym, w, mu, sig)
print('w,mu,sig', w,mu,sig)
psf = GaussianMixturePSF(w, mu, sig)
if bandname is None:
# try looking up in filtermap.
filt = phdr['FILTER']
if filt in filtermap:
print('Mapping filter', filt, 'to', filtermap[filt])
bandname = filtermap[filt]
else:
print('No mapping found for filter', filt)
bandname = flit
photocal = cf.CfhtPhotoCal(hdr=phdr, bandname=bandname)
filename = phdr['FILENAME'].strip()
from tractor.emfit import em_fit_2d
from tractor.fitpsf import em_init_params
# Create Gaussian mixture model PSF approximation.
S = psf.shape[0]
# number of Gaussian components
K = 3
w,mu,sig = em_init_params(K, None, None, None)
II = psf.copy()
II /= II.sum()
# HIDEOUS HACK
II = np.maximum(II, 0)
print('Multi-Gaussian PSF fit...')
xm,ym = -(S/2), -(S/2)
em_fit_2d(II, xm, ym, w, mu, sig)
print('w,mu,sig', w,mu,sig)
mypsf = tractor.GaussianMixturePSF(w, mu, sig)
P = mypsf.getPointSourcePatch(S/2, S/2)
mn,mx = psf.min(), psf.max()
ima = dict(interpolation='nearest', origin='lower',
vmin=mn, vmax=mx)
plt.clf()
plt.subplot(1,2,1)
plt.imshow(psf, **ima)
plt.subplot(1,2,2)
pimg = np.zeros_like(psf)
P.addTo(pimg)
plt.imshow(pimg, **ima)
plt.savefig('psf.png')
tim = Image(data=image, invvar=invvar, psf=tpsf, wcs=twcs, photocal=photocal,
sky=tsky, name='CFHT coadd %s %s' % (obj, bandname))
# set "zr" for plots
sig = 1./np.median(tim.inverr)
tim.zr = np.array([-1., +20.]) * sig
if not doplots:
return tim
psfimpatch = Patch(-(PW/2), -(PH/2), psfim)
# number of Gaussian components
for K in range(1, 4):
w,mu,sig = em_init_params(K, None, None, None)
xm,ym = -(PW/2), -(PH/2)
em_fit_2d(psfim, xm, ym, w, mu, sig)
#print 'w,mu,sig', w,mu,sig
psf = GaussianMixturePSF(w, mu, sig)
patch = psf.getPointSourcePatch(0, 0)
plt.clf()
plt.subplot(1,2,1)
plt.imshow(patch.getImage(), interpolation='nearest', origin='lower')
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow((patch - psfimpatch).getImage(), interpolation='nearest', origin='lower')
plt.colorbar()
plt.savefig('copsf-%i.png' % K)
plt.clf()
plt.imshow(psfim, interpolation='nearest', origin='lower')
plt.colorbar()
# print 'Read PSF image:', psf.shape, 'range', psf.min(), psf.max()
if positive:
psf = np.maximum(psf, 0.)
psf = PixelizedPSF(psf)
cache[key] = psf
return psf
S = psf.shape[0]
# number of Gaussian components
K = 3
w, mu, sig = em_init_params(K, None, None, None)
II = psf.copy()
II = np.maximum(II, 0)
II /= II.sum()
xm, ym = -(S / 2), -(S / 2)
res = em_fit_2d(II, xm, ym, w, mu, sig)
if res != 0:
raise RuntimeError('Failed to fit PSF')
print('W1 PSF:')
print(' w', w)
print(' mu', mu)
print(' sigma', sig)
psf = GaussianMixturePSF(w, mu, sig)
psf.computeRadius()
cache[key] = psf
return psf
invvar=table[1].data
skyobj = ba.ConstantSky(header['skyval'])
psffn = 'J082742.02+212844.7-r-bpsf.fits.gz'
psfimg = pyfits.open(psffn)[0].data
print('PSF image shape', psfimg.shape)
# number of Gaussian components
PS = psfimg.shape[0]
K = 3
w,mu,sig = em_init_params(K, None, None, None)
II = psfimg.copy()
II /= II.sum()
# HACK
II = np.maximum(II, 0)
print('Multi-Gaussian PSF fit...')
xm,ym = -(PS/2), -(PS/2)
em_fit_2d(II, xm, ym, w, mu, sig)
print('w,mu,sig', w,mu,sig)
psf = GaussianMixturePSF(w, mu, sig)
sources = []
# for run,camcol,field in zip(run,camcol,field):
# sources.append(st.get_tractor_sources(run,camcol,field,bandname,bands=bands))
wcs = Tan("J082742.02+212844.7-r.fits",0)
wcs = FitsWcs(wcs)
wcs.setX0Y0(1.,1.)
photocal = ba.NasaSloanPhotoCal(bandname) #Also probably not right
TI.append(en.Image(data=data,invvar=invvar,sky=skyobj,psf=psf,wcs=wcs,photocal=photocal,name = "NASA-Sloan Test"))
lvl = logging.DEBUG
logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)
tims = [TI[0]]
tractor = st.SDSSTractor(tims)
psfimg = pyfits.open(psffn)[0].data
print('PSF image shape', psfimg.shape)
from tractor.emfit import em_fit_2d
from tractor.fitpsf import em_init_params
# number of Gaussian components
S = psfimg.shape[0]
K = 3
w,mu,sig = em_init_params(K, None, None, None)
II = psfimg.copy()
II /= II.sum()
# HACK
II = np.maximum(II, 0)
print('Multi-Gaussian PSF fit...')
xm,ym = -(S/2), -(S/2)
em_fit_2d(II, xm, ym, w, mu, sig)
print('w,mu,sig', w,mu,sig)
psf = GaussianMixturePSF(w, mu, sig)
photocal = cf.CfhtPhotoCal(hdr=phdr,
bandname='r')
cftimg = Image(data=image, invvar=invvar, psf=psf, wcs=wcs,
sky=skyobj, photocal=photocal,
name='CFHT')
return cftimg, cfsky, cfstd
plt.subplot(3, 3, 3 * i + j + 1)
plt.imshow(np.log10(np.maximum(P, 1e-8)),
interpolation='nearest', origin='lower', vmax=0.01)
# plt.colorbar()
plt.savefig('psf-w1-xy.png')
psf = pyfits.open('psf%i.fits' % 1)[0].data
S = psf.shape[0]
# number of Gaussian components
for K in range(1, 6):
w, mu, sig = em_init_params(K, None, None, None)
II = psf.copy()
II /= II.sum()
II = np.maximum(II, 0)
xm, ym = -(S / 2), -(S / 2)
res = em_fit_2d(II, xm, ym, w, mu, sig)
print('em_fit_2d result:', res)
if res != 0:
raise RuntimeError('Failed to fit PSF')
print('w,mu,sig', w, mu, sig)
mypsf = GaussianMixturePSF(w, mu, sig)
mypsf.computeRadius()
#
mypsf.radius = S / 2
mod = mypsf.getPointSourcePatch(0., 0.)
mod = mod.patch
mod /= mod.sum()
plt.clf()
plt.subplot(1, 2, 1)
ima = dict(interpolation='nearest', origin='lower',