In [1]:
import numpy as np
import os
import pandas as pd
import itertools
from PIL import Image
import seaborn as sns
import fitsio
import skimage.io
import galsim
from astrometry.util.fits import fits_table, merge_tables
# to make this notebook's output stable across runs
np.random.seed(7)
# To plot pretty figures
%matplotlib inline
#%matplotlib notebook
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
%load_ext autoreload
%autoreload 2
In [2]:
from obiwan.qa.visual import readImage,sliceImage,plotImage, flux2mag
In [3]:
REPO_DIR= os.path.join(os.environ['HOME'],
'myrepo/obiwan')
DATA_DIR= os.path.join(os.environ['HOME'],
'mydata')
In [9]:
gal_flux = 1.e5 # total counts on the image
rhalf = 2. # arcsec
psf_sigma = 1. # arcsec
pixel_scale = 0.2 # arcsec / pixel
noise = 30. # standard
gal = galsim.Sersic(1,flux=gal_flux, half_light_radius=rhalf)
gal = gal.shear(e1=0.3, e2=-0.8)
rot_gal= gal.rotate(np.pi/2.* galsim.radians)
psf = galsim.Gaussian(flux=1., sigma=psf_sigma) # PSF flux should always = 1
final = galsim.Convolve([gal, psf])
image = final.drawImage(scale=pixel_scale)
image.addNoise(galsim.GaussianNoise(sigma=noise))
In [11]:
fig,ax= plt.subplots(2,3,figsize=(6,3))
kw={'scale':0.26}
plotImage().imshow(gal.drawImage(**kw).array,ax[0,0])
plotImage().imshow(rot_gal.drawImage(**kw).array,ax[0,1])
plotImage().imshow(psf.drawImage(**kw).array,ax[0,2])
plotImage().imshow(final.drawImage(**kw).array,ax[1,0])
plotImage().imshow(image.array,ax[1,1])
In [12]:
ba= np.random.rand(1000)
ba= ba[ba > 0.1]
beta= np.random.rand(len(ba))*180
In [13]:
def ellip(q):
"""q: b/a ratio"""
return (1-q**2)/(1+q**2)
def e1_e2(q,beta):
"""
q: b/a ratio
beta: rotation
"""
e= ellip(q)
return e*np.cos(2*beta), e*np.sin(2*beta)
fig,ax= plt.subplots(1,2,figsize=(6,3))
ax[0].scatter(ba,beta,alpha=0.5)
e1,e2= e1_e2(ba,beta)
ax[1].scatter(e1,e2,alpha=0.5)
Out[13]:
In [4]:
from legacypipe.decam import DecamImage
from legacypipe.survey import LegacySurveyData
In [5]:
ccds= fits_table(os.path.join(DATA_DIR,
'1741p242/dr5/legacysurvey-1741p242-ccds.fits'))
t= ccds[((pd.Series(ccds.image_filename).str.contains('c4d_160116_084245_oki_z_v1.fits')) &
(ccds.ccdname == 'S17'))]
#for col in ['image_filename','ccdname','filter','camera']:
# t.set(col,t.get(col)[0])
#t
In [6]:
print(type(ccds))
print(type(t))
for ccd in t:
print(type(ccd),ccd)
In [7]:
# Generator class for LegacySurveyImage() objects
survey = LegacySurveyData(ccds=ccds)
#survey= DecamImage(ccds=)
for ccd in t:
# a LegacySurveyImage
im = survey.get_image_object(ccd)
#X = im.get_good_image_subregion()
kwargs = dict(pixPsf=True, splinesky=True, subsky=False, hybridPsf=True,
pixels=True, dq=True, invvar=True)
tim = im.get_tractor_image(**kwargs)
In [12]:
tim.getImage().max(),tim.invvar.shape,tim.data.shape,tim.dq.shape
Out[12]:
In [36]:
def get_wcs(tim):
from astropy.io import fits
temp_hdr = fitsio.FITSHDR()
subwcs = tim.wcs.wcs.get_subimage(tim.wcs.x0, tim.wcs.y0,
int(tim.wcs.wcs.get_width())-tim.wcs.x0,
int(tim.wcs.wcs.get_height())-tim.wcs.y0)
subwcs.add_to_header(temp_hdr)
# Galsim uses astropy header, not fitsio
hdr = fits.Header()
for key in temp_hdr.keys():
hdr[key]=temp_hdr[key]
return galsim.GSFitsWCS(header=hdr)
psf = tim.psf.getPointSourcePatch(500,500)
const_psf= tim.psf.constantPsfAt(500,500)
psf= galsim.Image(psf.getImage(),wcs=get_wcs(tim))
pos = galsim.PositionD(500,500)
psf_interp= galsim.InterpolatedImage(psf, wcs=get_wcs(tim))
star = psf_interp.drawImage(scale=0.262) #wcs=get_wcs(tim).local(image_pos=pos), method='no_pixel')
In [71]:
gal.drawImage?
In [75]:
galsim.InterpolatedImage?
In [58]:
gal = galsim.Sersic(4, half_light_radius=0.37,\
flux=1., gsparams= galsim.GSParams(maximum_fft_size=65536))
gal = gal.shear(e1=0.33,e2=0.27)
gal_nopix = gal.drawImage(scale=tim.wcs.pixel_scale(),method='no_pixel')
gal_auto= gal.drawImage(scale=tim.wcs.pixel_scale(),method='auto')
In [68]:
fig,ax= plt.subplots(1,2,figsize=(4,3))
xslc,yslc= slice(35,55),slice(35,55)
plotImage().imshow(sliceImage(gal_nopix.array,xslice=xslc,yslice=yslc),ax[0])
plotImage().imshow(sliceImage(gal_auto.array,xslice=xslc,yslice=yslc),ax[1])
print(gal_nopix.array.sum(),gal_auto.array.sum())
In [70]:
tim.wcs.pixscale_at(500,500)
Out[70]:
In [66]:
galsim.InterpolatedImage?
In [34]:
fig,ax=plt.subplots(5,5,figsize=(10,10))
for i,x in enumerate(np.linspace(100,4000,num=5).astype(int)):
for j,y in enumerate(np.linspace(100,2046,num=5).astype(int)):
psf=tim.psf.getPointSourcePatch(x,y).patch
plotImage().imshow(sliceImage(psf,xslice=slice(20,44),yslice=slice(20,44)),
ax[i,j],qs=None)
print(psf.sum())
In [44]:
fig,ax=plt.subplots(5,5,figsize=(10,10))
for i,x in enumerate(np.linspace(100,4000,num=5).astype(int)):
for j,y in enumerate(np.linspace(100,2046,num=5).astype(int)):
psf=tim.psf.getPointSourcePatch(x,y).patch
psf= galsim.Image(psf,wcs=get_wcs(tim))
plotImage().imshow(sliceImage(psf.array,xslice=slice(20,44),yslice=slice(20,44)),
ax[i,j],qs=None)
print(psf.array.sum())
In [29]:
gal = galsim.Sersic(1, half_light_radius=2,flux=1.)
gal = galsim.Convolve([gal, psf_interp])
In [13]:
psf = tim.psf
wcs = tim.wcs.wcs
sky = tim.sky
In [ ]:
In [76]:
galsim.InterpolatedImage?
In [77]:
def magToNanomaggies(mag):
nmgy = 10. ** ((mag - 22.5) / -2.5)
return nmgy
def global_wcs(tim):
from astropy.io import fits
#self.wcs = tim.getWcs()
#self.psf = tim.getPsf()
# Tractor wcs object -> galsim wcs object
temp_hdr = fitsio.FITSHDR()
subwcs = tim.wcs.wcs.get_subimage(tim.wcs.x0, tim.wcs.y0,
int(tim.wcs.wcs.get_width())-tim.wcs.x0,
int(tim.wcs.wcs.get_height())-tim.wcs.y0)
subwcs.add_to_header(temp_hdr)
# Galsim uses astropy header, not fitsio
hdr = fits.Header()
for key in temp_hdr.keys():
hdr[key]=temp_hdr[key]
return galsim.GSFitsWCS(header=hdr)
def local_wcs(global_wcs,x,y):
return global_wcs.local(image_pos= galsim.PositionD(x,y))
def galsim_psf(x,y,tim,global_wcs,gsparams=None):
psf_tim= tim.psf.getPointSourcePatch(x,y)
psf_tim.patch /= psf_tim.patch.sum()
return galsim.InterpolatedImage(galsim.Image(psf_tim.getImage()),
wcs= global_wcs, gsparams=gsparams)
def galsim_galaxy(x,y,tim,global_wcs,
n=1,rhalf=2.,mag=22.,
e1=0.3,e2=-0.8,angle=np.pi/2,
gsparams=None):
psf= galsim_psf(x,y,tim,global_wcs,gsparams=gsparams)
gal = galsim.Sersic(n,flux=magToNanomaggies(mag), half_light_radius=rhalf,
gsparams=gsparams)
gal= gal.shear(e1=e1, e2=e2)
gal= gal.rotate(angle* galsim.radians)
return galsim.Convolve([gal, psf],gsparams=gsparams)
def Draw(obj,x,y,global_wcs):
"""obj: galsim object having method drawImage()"""
image= obj.drawImage(wcs=local_wcs(global_wcs,x,y), method='no_pixel')
# assign to a location
image.setCenter(x,y)
return image
wcs= global_wcs(tim)
xpos,ypos,hw= 1000,1000,30
xslc=slice(xpos-hw,xpos+hw)
yslc=slice(ypos-hw,ypos+hw)
psf= galsim_psf(xpos,ypos,tim,wcs)
psf= Draw(psf,xpos,ypos,wcs)
gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=1,rhalf=1.,mag=20.,
e1=0.3,e2=-0.8,angle=np.pi/2)
gal= Draw(gal,xpos,ypos,wcs)
tim_galsim = galsim.Image(tim.getImage())
overlap= {'gal':tim_galsim.bounds & gal.bounds,
'psf':tim_gablsim.bounds & psf.bounds}
psf.bounds,overlap['psf'],gal.bounds,overlap['gal']
Out[77]:
In [82]:
galsim.Image?
In [15]:
hdu=fitsio.FITS('images/decam/DECam_CP/CP20160107/c4d_160116_084245_oki_z_v1.fits.fz')
cp=hdu['S17'].read()
fig,ax= plt.subplots(2,3,figsize=(7,4))
##
olap= overlap['gal']
plotImage().imshow(sliceImage(cp,xslice=slice(olap.xmin,olap.xmax),
yslice=slice(olap.ymin,olap.ymax)),
ax[0,0])
plotImage().imshow(tim_galsim[olap].array,ax[0,1])
plotImage().imshow(gal[olap].array,ax[0,2])
##
olap= overlap['psf']
plotImage().imshow(sliceImage(cp,xslice=slice(olap.xmin,olap.xmax),
yslice=slice(olap.ymin,olap.ymax)),
ax[1,0])
plotImage().imshow(tim_galsim[olap].array,ax[1,1])
plotImage().imshow(psf[olap].array,ax[1,2])
In [8]:
def addNoise(image,zpscale=None,gain=4.,seed=7):
"""image: source after it has undergone drawImage()"""
nano2e= zpscale * gain
src_e= image.copy() * nano2e
var_e= (image.copy() * nano2e)**2
rng = galsim.BaseDeviate(7)
src_e.addNoise(galsim.VariableGaussianNoise(rng,var_e))
src_cnts= src_e / nano2e
var_cnts= var_e / nano2e**2
return src_cnts,var_cnts
In [17]:
olap= overlap['gal']
tim_wgal= tim_galsim.copy()
tim_wgal[olap] += gal[olap]
gal_wnoise, gal_wnoise_var= addNoise(gal,zpscale=tim.zpscale, gain=1.)
tim_wgalnoise= tim_galsim.copy()
tim_wgalnoise[olap] += gal_wnoise[olap]
fig,ax= plt.subplots(2,4,figsize=(8,4))
plotImage().imshow(tim_galsim[olap].array,ax[0,0])
plotImage().imshow(tim_wgal[olap].array,ax[0,1])
plotImage().imshow(tim_galsim[olap].array - tim_wgal[olap].array,ax[0,2])
plotImage().imshow(gal[olap].array,ax[0,3])
##
plotImage().imshow(tim_galsim[olap].array,ax[1,0])
plotImage().imshow(tim_wgalnoise[olap].array,ax[1,1])
plotImage().imshow(tim_galsim[olap].array - tim_wgalnoise[olap].array,ax[1,2])
plotImage().imshow(gal_wnoise[olap].array,ax[1,3])
In [18]:
fig,ax= plt.subplots(figsize=(2,2))
noise_only= gal_wnoise[olap].array - gal[olap].array
plotImage().imshow(noise_only,ax)
In [19]:
fig,ax= plt.subplots(1,2,figsize=(6,3))
ax[0].scatter(noise_only,gal[olap].array.flatten())
ax[1].scatter(noise_only,np.sqrt(gal[olap].array.flatten()))
Out[19]:
In [14]:
dr5tractor= fits_table(os.path.join(DATA_DIR,
'1741p242','dr5',
'tractor-1741p242.fits'))
img_jpg= readImage(os.path.join(DATA_DIR,
'1741p242','dr5',
'legacysurvey-1741p242-image.jpg'),
jpeg=True)
ccds= fits_table(os.path.join(DATA_DIR,
'1741p242','dr5',
'ccds-annotated-dr5-near-1741p242.fits'))
# 'legacysurvey-1741p242-ccds.fits'))
In [15]:
oneBand= (dr5tractor.nobs_g + dr5tractor.nobs_r + dr5tractor.nobs_z) == 1
isGal= dr5tractor.get('type') != 'PSF '
dr5tractor.cut(((oneBand) &
(isGal)))
len(dr5tractor)
Out[15]:
In [16]:
hw=30
for trac in dr5tractor:
fig,ax=plt.subplots(figsize=(2,2))
x,y= int(trac.bx),int(trac.by)
print(x,y,trac.flux_g,trac.flux_r,trac.flux_z)
plotImage().imshow(sliceImage(img_jpg,xslice=slice(x-hw,x+hw),yslice=slice(y-hw,y+hw)),
ax)
In [17]:
trac= dr5tractor[-2]
fig,ax=plt.subplots(figsize=(3,3))
x,y= int(trac.bx),int(trac.by)
plotImage().imshow(sliceImage(img_jpg,xslice=slice(x-hw,x+hw),yslice=slice(y-hw,y+hw)),
ax)
In [18]:
def flux2mag(flux):
return -2.5*np.log10(1e-9 * flux)
def mag2flux(mag):
return 10**(9-mag/2.5)
print(flux2mag(trac.flux_z), trac.mw_transmission_z, trac.ra,trac.dec)
print(trac.type, trac.shapedev_r, trac.shapedev_e1, trac.shapedev_e2)
In [19]:
# Which ccd is this source in?
ra,dec= trac.ra,trac.dec
ramax= np.max([ccds.ra1,ccds.ra2],axis=0)
ramin= np.min([ccds.ra0,ccds.ra3],axis=0)
decmax= np.max([ccds.dec1,ccds.dec0],axis=0)
decmin= np.min([ccds.dec3,ccds.dec3],axis=0)
ccds.cut((ra > ramin) &
(ra < ramax) &
(dec > decmin) &
(dec < decmax))
len(ccds),ccds.image_filename[0],ccds.ccdname[0]
Out[19]:
In [22]:
ccds= fits_table(os.path.join(DATA_DIR,
'1741p242/dr5/legacysurvey-1741p242-ccds.fits'))
t= ccds[((pd.Series(ccds.image_filename).str.contains('c4d_150329_031648_ooi_z_v1.fits.fz')) &
(ccds.ccdname == 'S16'))]
survey = LegacySurveyData(ccds=ccds)
#survey= DecamImage(ccds=)
for ccd in t:
# a LegacySurveyImage
im = survey.get_image_object(ccd)
#X = im.get_good_image_subregion()
kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
pixels=True, dq=True, invvar=True)
tim = im.get_tractor_image(**kwargs)
In [23]:
tim.wcs.wcs.radec2pixelxy(trac.ra,trac.dec), tim.getImage().shape
Out[23]:
In [24]:
dpix=25
x,y=1040.23,3722.59
newra,newdec=[],[]
for dx,dy in zip([dpix,-dpix,0,0],[0,0,dpix,-dpix]):
newx,newy=x+dx,y+dy
ra,dec= tim.wcs.wcs.pixelxy2radec(newx,newy)
print(newx,newy,ra,dec)
newra.append(ra)
newdec.append(dec)
print(newra)
print(newdec)
In [31]:
dpix=50
x,y=1040.23,3722.59
newra,newdec=[],[]
for dx,dy in zip([dpix,-dpix,0,0],[0,0,dpix,-dpix]):
newx,newy=x+dx,y+dy
ra,dec= tim.wcs.wcs.pixelxy2radec(newx,newy)
print(newx,newy,ra,dec)
newra.append(ra)
newdec.append(dec)
print(newra)
print(newdec)
In [32]:
print(trac.objid,flux2mag(trac.flux_z), trac.mw_transmission_z, trac.ra,trac.dec)
print(trac.type, trac.shapedev_r, trac.shapedev_e1, trac.shapedev_e2)
In [25]:
wcs= global_wcs(tim)
xpos,ypos,hw= 1040,3723,30
psf= galsim_psf(xpos,ypos,tim,wcs)
psf= Draw(psf,xpos,ypos,wcs)
gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=4,rhalf=0.3696,mag=18.543,
e1=0.334,e2=0.273,angle=0)
gal= Draw(gal,xpos,ypos,wcs)
gal_wnoise,_= addNoise(gal,zpscale=tim.zpscale, gain=4.)
# Move fake galaxy over by an offset so not on top of real one
off= 0
gal.setCenter(xpos+off,ypos)
gal_wnoise.setCenter(xpos+off,ypos)
tim_galsim = galsim.Image(tim.getImage())
olap= tim_galsim.bounds & gal.bounds
psf.bounds,gal.bounds,olap
In [ ]:
fig,ax= plt.subplots(1,2,figsize=(6,3))
hw=5
zoom= galsim.BoundsI(xmin=xpos-hw,xmax=xpos+hw,ymin=ypos-hw,ymax=ypos+hw)
plotImage().imshow(tim_galsim[zoom].array,ax[0],qs=None)
plotImage().imshow(gal[zoom].array,ax[1],qs=None)
In [18]:
# Correct noise model
one_std_per_pix= gal.array.copy()
one_std_per_pix= np.sqrt(one_std_per_pix)
num_stds= np.random.randn(one_std_per_pix.shape[0],one_std_per_pix.shape[1])
#one_std_per_pix.shape, num_stds.shape
noise= one_std_per_pix * num_stds
gal_noise= galsim.Image(noise)
gal_noise.setCenter(xpos,ypos)
In [19]:
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(gal.array,ax[0],qs=None)
plotImage().imshow(gal_noise.array,ax[1],qs=None)
In [20]:
# Fake galaxy can have negative pixel vals (~ few percent of the pixels)
gal.array[gal.array < 0].shape, gal.array.flatten().shape
Out[20]:
In [21]:
# Noise model + no negative image vals when compute noise
one_std_per_pix= gal.array.copy()
one_std_per_pix[one_std_per_pix < 0]=0
one_std_per_pix= np.sqrt(one_std_per_pix)
num_stds= np.random.randn(one_std_per_pix.shape[0],one_std_per_pix.shape[1])
#one_std_per_pix.shape, num_stds.shape
noise= one_std_per_pix * num_stds
gal_noise= galsim.Image(noise)
In [22]:
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(gal.array,ax[0],qs=None)
plotImage().imshow(gal_noise.array,ax[1],qs=None)
In [35]:
plot3d(gal_noise.array)
In [36]:
plot3d((gal+gal_noise).array)
In [23]:
gal_wnoise= gal + gal_noise
gal_wnoise.setCenter(xpos+15,ypos)
tim_wgal_wnoise= tim_galsim.copy()
olap= tim_wgal.bounds & gal_wnoise.bounds
tim_wgal_wnoise[olap] += gal_wnoise[olap]
tim_wgal= tim_galsim.copy()
gal2= gal.copy()
gal2.setCenter(xpos+15,ypos)
tim_wgal[olap] += gal2[olap]
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(tim_wgal[olap].array,ax[0],qs=None)
plotImage().imshow(tim_wgal_wnoise[olap].array,ax[1],qs=None)
for i in range(2):
ax[i].set_xlim(25,60)
ax[i].set_ylim(40,60)
#plotImage().imshow(gal_noise.array,ax[1],qs=None)
plt.savefig('2d.png',dpi=150)
In [38]:
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)
import plotly.graph_objs as go
def plot3d(image,zlim=None):
data = [
go.Surface(
z=image
)
]
layout=go.Layout()
if zlim:
layout = go.Layout(
scene= dict(
zaxis = dict(
nticks=10, range = [zlim[0],zlim[1]]),
)
)
fig = go.Figure(data=data, layout=layout)
#fig = go.Figure(data=data)
iplot(fig, filename='image')
In [40]:
plot3d(tim_wgal[olap].array)
In [41]:
X=np.random.randn(100,100)
Y= 4*X
np.mean(Y)/np.mean(X),np.std(Y)/np.std(X), tim.zpscale
Out[41]:
In [59]:
Out[59]:
In [27]:
def nano2e(zpscale,gain=1.):
return zpscale * gain
def noise_for_galaxy(gal,zpscale,gain=1.):
"""gal: galsim stamp for galaxy"""
f= nano2e(zpscale,gain=gain)
# Noise model + no negative image vals when compute noise
one_std_per_pix= gal.array.copy() # nanomaggies
one_std_per_pix[one_std_per_pix < 0]=0
# rescale
one_std_per_pix *= f # e-
one_std_per_pix= np.sqrt(one_std_per_pix)
num_stds= np.random.randn(one_std_per_pix.shape[0],one_std_per_pix.shape[1])
#one_std_per_pix.shape, num_stds.shape
noise= one_std_per_pix * num_stds
# rescale
noise /= f #nanomaggies
return galsim.Image(noise)
def var_for_galaxy(gal,zpscale,gain=1.):
assert(type(gal) == galsim.Image)
f= nano2e(zpscale,gain=gain)
# Noise model + no negative image vals when compute noise
var_= gal.array.copy() # nanomaggies
var_per_pix[var_per_pix < 0]=0
# rescale
var_per_pix *= f # variance in nanomaggies
return galsim.Image(var_per_pix)
def remap_invvar_var_galaxy(gal,invvar,dq,
zpscale,gain=1.):
"""add variance from simulated galaxy to invvar map
Args:
gal,invvar,dq: galsim image objects overlapping region of gal
invvar,dq from CP oow and ood
zpscale,gain
"""
assert(type(gal) == galsim.Image)
assert(type(invvar) == galsim.Image)
assert(type(dq) == galsim.Image)
f= nano2e(zpscale,gain=gain)
var_gal= gal.array.copy() * f**2
var_gal[var_gal < 0]=0.
with np.errstate(divide='ignore'):
var = 1./invvar.array.copy() + var_gal
#var += var_for_galaxy(gal,zpscale,gain=gain).array
wt = 1./var # s/e
# Zero out NaNs and masked pixels
wt[np.isfinite(wt) == False] = 0.
wt[dq != 0] = 0.
return galsim.Image(wt)
In [58]:
galsim.VariableGaussianNoise?
In [56]:
a=galsim.Image(np.random.randint(-1,1,(100,100)).astype(np.float32))
sns.distplot(a.array.flatten(),color='b')
print(len(np.where(a.array < 0)[0]))
a.applyNonlinearity(np.abs)
sns.distplot(a.array.flatten(),color='r')
print(len(np.where(a.array < 0)[0]))
a += np.ones(a.array.shape)
sns.distplot(a.array.flatten(),color='g')
Out[56]:
In [60]:
np.random.randn?
In [59]:
sns.distplot(np.random.randn(100,100).flatten())
Out[59]:
In [41]:
#tim_wgal.array= np.abs(tim_wgal.array)
print(len(tim_wgal.array[tim_wgal < 0]))
tim_wgal.applyNonlinearity(np.abs)
In [38]:
# Gain = 1
gal_noise= noise_for_galaxy(gal,tim.zpscale,gain=1)
gal_wnoise= gal + gal_noise
gal_wnoise.setCenter(xpos+15,ypos)
tim_wgal_wnoise_1= tim_galsim.copy()
olap= tim_galsim.bounds & gal_wnoise.bounds
tim_wgal_wnoise_1[olap] += gal_wnoise[olap]
tim_wgal= tim_galsim.copy()
gal2= gal.copy()
gal2.setCenter(xpos+15,ypos)
tim_wgal[olap] += gal2[olap]
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(tim_wgal[olap].array,ax[0],qs=None)
plotImage().imshow(tim_wgal_wnoise_1[olap].array,ax[1],qs=None)
for i in range(2):
ax[i].set_xlim(25,60)
ax[i].set_ylim(40,60)
#plotImage().imshow(gal_noise.array,ax[1],qs=None)
In [34]:
# Gain = 4
gal_noise_4= noise_for_galaxy(gal,tim.zpscale,gain=4)
gal_wnoise_4= gal + gal_noise
gal_wnoise_4.setCenter(xpos+15,ypos)
tim_wgal_wnoise_4= tim_galsim.copy()
olap= tim_galsim.bounds & gal_wnoise_4.bounds
tim_wgal_wnoise_4[olap] += gal_wnoise_4[olap]
tim_wgal= tim_galsim.copy()
gal2= gal.copy()
gal2.setCenter(xpos+15,ypos)
tim_wgal[olap] += gal2[olap]
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(tim_wgal[olap].array,ax[0],qs=None)
plotImage().imshow(tim_wgal_wnoise_4[olap].array,ax[1],qs=None)
for i in range(2b
ax[i].set_xlim(25,60)
ax[i].set_ylim(40,60)
#plotImage().imshow(gal_noise.array,ax[1],qs=None)
In [24]:
wcs= global_wcs(tim)
xpos,ypos,hw= 1040,3723,30
In [149]:
import photutils
In [152]:
psf.trueCenter()
Out[152]:
In [157]:
Out[157]:
In [179]:
psf= tim.psf.getPointSourcePatch(xpos,ypos)
psf.patch.shape
Out[179]:
In [52]:
photutils.CircularAperture?
In [46]:
%matplotlib notebook
import photutils
def normed_psf(xpos,ypos,wcs):
"""
Returns:
galsim Image() of PSF normalized at 7''
"""
psf= tim.psf.getPointSourcePatch(xpos,ypos)
#psf= galsim.Image(psf.getImage(),wcs=wcs)
#psfim /= psfim.sum()
# Now make sum to 1 at 7''
psf= galsim.Image(psf.getImage(),wcs=wcs)
apers= photutils.CircularAperture((psf.trueCenter().x,psf.trueCenter().y),
r=3.5/0.262) #KEY is to harcode this # pix
apy_table = photutils.aperture_photometry(psf.array, apers)
flux_in_7= np.array(apy_table['aperture_sum'])[0]
frac_in_7= flux_in_7 / psf.array.sum()
#print('flux_in_7=',flux_in_7)
psf /= flux_in_7
#psfim /= frac_in_7
return psf
def normed_galaxy(xpos,ypos,wcs,local_wcs):
"""
Returns:
galsim Image object of normalized-PSF-convolved Galaxy profile, which iteslf
is normalized at 7''
"""
#print(psf.array.sum())
#flux=magToNanomaggies(18.543)
gal = galsim.Sersic(4.,flux=1., half_light_radius=0.37)
gal= gal.shear(e1=0.334,e2=0.273)
psf= normed_psf(xpos,ypos,wcs)
psf= galsim.InterpolatedImage(psf,wcs=wcs)
gal= galsim.Convolve(gal,psf)
gal= gal.drawImage(wcs=local_wcs,method='no_pixel')
# Normalize
apers= photutils.CircularAperture((gal.trueCenter().x,gal.trueCenter().y),
r=3.5/0.262)
apy_table = photutils.aperture_photometry(gal.array, apers)
flux_in_7= np.array(apy_table['aperture_sum'])[0]
#print('flux_in_7=',flux_in_7)
gal /= flux_in_7
#gal.setCenter(xpos,ypos)
return gal
wcs= global_wcs(tim)
psf= normed_psf(xpos,ypos,global_wcs(tim))
gal= normed_galaxy(xpos,ypos,wcs,local_wcs(wcs,xpos,ypos))
# Test that they are normalized
#psf_img= psf.drawImage(wcs=local_wcs(wcs,xpos,ypos),method='no_pixel')
#psf.trueCenter(),gal.trueCenter()
apers= photutils.CircularAperture((psf.trueCenter().x,psf.trueCenter().y),
r=3.5/0.262)
apy_table = photutils.aperture_photometry(psf.array, apers)
print(np.array(apy_table['aperture_sum'])[0])
apers= photutils.CircularAperture((gal.trueCenter().x,gal.trueCenter().y),
r=3.5/0.262)
apy_table = photutils.aperture_photometry(gal.array, apers)
print(np.array(apy_table['aperture_sum'])[0])
#gal_norm= galsim_psf_galaxy(xpos,ypos,wcs,norm_at_7=True)
fig,ax=plt.subplots(1,3,figsize=(9,3))
plotImage().imshow(psf.array,ax[0],qs=None)
plotImage().imshow(gal.array,ax[1],qs=None)
gal *= mag2flux(18.543)
apy_table = photutils.aperture_photometry(gal.array, apers)
print(np.array(apy_table['aperture_sum'])[0])
plotImage().imshow(gal.array,ax[2],qs=None)
# print(gal_norm.array.sum() - gal.array.sum())
# fig,ax=plt.subplots(1,2,figsize=(6,3))
# plotImage().imshow(gal.array,ax[0],qs=None)
# plotImage().imshow(gal_norm.array,ax[1],qs=None)
In [178]:
print(gal.array.sum())
gal[gal.bounds]= gal[gal.bounds]/2
print(gal.array.sum())
In [142]:
gal = galsim.Exponential(flux=1.e5, scale_radius=2.7)
#psf = galsim.Moffat(beta=5, flux=1., half_light_radius=1.)
conv = galsim.Convolve([gal, psf])
pxscale=0.262
gal= gal.drawImage(scale=pxscale)
psf= psf.drawImage(scale=pxscale)
conv= conv.drawImage(scale=pxscale)
fig,ax=plt.subplots(1,3,figsize=(12,3))
plotImage().imshow(gal.array,ax[0],qs=None)
plotImage().imshow(psf.array,ax[1],qs=None)
plotImage().imshow(conv.array,ax[2],qs=None)
In [118]:
no_noise,w_noise=[],[]
for i in range(100):
base_gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=4,rhalf=0.3696,mag=18.543,
e1=0.334,e2=0.273,angle=0)
gal= Draw(base_gal,xpos,ypos,wcs)
no_noise.append(gal.array.sum())
gal += noise_for_galaxy(gal,tim.zpscale,gain=4)
w_noise.append(gal.array.sum())
no_noise,w_noise= np.array(no_noise),np.array(w_noise)
fig,ax=plt.subplots(1,3,figsize=(12,3))
plotImage().imshow(gal.array,ax[0],qs=None)
sns.distplot(w_noise - mag2flux(18.543),ax=ax[1])
ax[1].set_xlabel('Flux (wnoise - Truth)')
sns.distplot(flux2mag(w_noise) - 18.543,ax=ax[2])
ax[2].set_xlabel('Mag (wnoise - Truth)')
Out[118]:
In [119]:
gsparams = galsim.GSParams(maximum_fft_size=1048576L,\
folding_threshold=1.e-5)
no_noise,w_noise=[],[]
for i in range(100):
base_gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=4,rhalf=0.3696,mag=18.543,
e1=0.334,e2=0.273,angle=0)
gal= Draw(base_gal,xpos,ypos,wcs)
no_noise.append(gal.array.sum())
gal += noise_for_galaxy(gal,tim.zpscale,gain=4)
w_noise.append(gal.array.sum())
no_noise,w_noise= np.array(no_noise),np.array(w_noise)
fig,ax=plt.subplots(1,3,figsize=(12,3))
plotImage().imshow(gal.array,ax[0],qs=None)
sns.distplot(w_noise - mag2flux(18.543),ax=ax[1])
ax[1].set_xlabel('Flux (wnoise - Truth)')
sns.distplot(flux2mag(w_noise) - 18.543,ax=ax[2])
ax[2].set_xlabel('Mag (wnoise - Truth)')
Out[119]:
In [120]:
gsparams = galsim.GSParams(folding_threshold=1.e-7)
no_noise,w_noise=[],[]
for i in range(100):
base_gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=4,rhalf=0.3696,mag=18.543,
e1=0.334,e2=0.273,angle=0)
gal= Draw(base_gal,xpos,ypos,wcs)
no_noise.append(gal.array.sum())
gal += noise_for_galaxy(gal,tim.zpscale,gain=4)
w_noise.append(gal.array.sum())
no_noise,w_noise= np.array(no_noise),np.array(w_noise)
fig,ax=plt.subplots(1,3,figsize=(12,3))
plotImage().imshow(gal.array,ax[0],qs=None)
sns.distplot(w_noise - mag2flux(18.543),ax=ax[1])
ax[1].set_xlabel('Flux (wnoise - Truth)')
sns.distplot(flux2mag(w_noise) - 18.543,ax=ax[2])
ax[2].set_xlabel('Mag (wnoise - Truth)')
Out[120]:
In [73]:
gal_wnoise_4.bounds,invvar_galsim.bounds,dq_galsim.bounds,olap
Out[73]:
In [146]:
tim_galsim[box]
Out[146]:
In [ ]:
In [40]:
tim_invvar_nosrc= invvar_galsim[olap].copy()
hw=10
box=galsim.BoundsI(xmin=xpos-hw,xmax=xpos+hw,ymin=ypos-hw,ymax=ypos+hw)
skybox=galsim.BoundsI(xmin=xpos-hw+25,xmax=xpos+hw+25,ymin=ypos-hw+25,ymax=ypos+hw+25)
tim_invvar_nosrc[box]= tim_invvar_nosrc[skybox]
var_nosrc= tim_invvar_nosrc.copy()
var_nosrc.invertSelf()
fig,ax= plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(var_nosrc.array,ax[0,0],qs=None)
var_nosrc[box] += galsim.Image( np.abs(tim_galsim[box].array) )
plotImage().imshow(var_nosrc.array,ax[0,1],qs=None)
ivar_wsrc= var_nosrc.copy()
ivar_wsrc.invertSelf()
plotImage().imshow(ivar_wsrc.array,ax[1,0],qs=None)
plotImage().imshow(invvar_galsim[olap].array,ax[1,1],qs=None)
In [41]:
sns.distplot(np.log10(var_nosrc.array.flatten()))
In [42]:
#f=tim.zpscale * 4.
fig,ax= plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(tim_galsim[olap].array,ax[0,0],qs=None)
plotImage().imshow(invvar_galsim[olap].array,ax[0,1],qs=None)
plotImage().imshow(tim_invvar_nosrc.array,ax[1,0],qs=None)
plotImage().imshow(ivar_wsrc.array,ax[1,1],qs=None)
In [128]:
plt.clf()
plt.scatter(tim_galsim[olap].array.flatten(),invvar_galsim[olap].array.flatten(),alpha=0.5)
Out[128]:
In [43]:
def remap_invvar_var_galaxy(gal,invvar,dq,
zpscale,gain=1.):
"""add variance from simulated galaxy to invvar map
Args:
gal,invvar,dq: galsim image objects overlapping region of gal
invvar,dq from CP oow and ood
zpscale,gain
"""
assert(type(gal) == galsim.Image)
assert(type(invvar) == galsim.Image)
assert(type(dq) == galsim.Image)
f= nano2e(zpscale,gain=gain)
var_gal= gal.array.copy() # nanomaggies
var_gal[var_gal < 0]=0.
with np.errstate(divide='ignore'):
var_real = 1./invvar.array.copy() # nano
var_nano= (var_real + var_gal) #* f**2 # e-
#var += var_for_galaxy(gal,zpscale,gain=gain).array
wt_nano = 1./var_nano
# Zero out NaNs and masked pixels
#wt_nano[np.isfinite(wt_nano) == False] = 0.
#wt_nano[dq != 0] = 0.
return galsim.Image(wt_nano) # nano
def invvar_galaxy(gal,zpscale,gain=1.):
assert(type(gal) == galsim.Image)
f= nano2e(zpscale,gain=gain)
var_e= np.abs(gal.array.copy()) * f # nano
var_nano= galsim.Image(var_e / f**2)
var_nano.invertSelf()
return var_nano
invvar_galsim = galsim.Image(tim.invvar)
dq_galsim = galsim.Image(tim.dq)
ivar_gal2= invvar_galaxy(gal2,tim.zpscale,gain=4.)
fig,ax= plt.subplots(1,3,figsize=(7,3))
plotImage().imshow(dq_galsim[olap].array,ax[0],qs=None)
plotImage().imshow(invvar_galsim[olap].array,ax[1],qs=None)
plotImage().imshow(ivar_gal2.array,ax[2],qs=None)
#invvar_galsim[olap] = remap_invvar_var_galaxy(gal2[olap],
# invvar_galsim[olap],dq_galsim[olap],
# tim.zpscale,gain=4.)
#plotImage().imshow(invvar_galsim[olap].array,ax[2],qs=None)
In [61]:
tim.invvar
Out[61]:
In [52]:
plot3d(tim_wgal_wnoise_1[olap].array)
In [54]:
plot3d(tim_wgal_wnoise_4[olap].array)
In [44]:
noise_1= noise_for_galaxy(gal,tim.zpscale,gain=1)
noise_4= noise_for_galaxy(gal,tim.zpscale,gain=4)
fig,ax= plt.subplots(1,2,figsize=(6,3))
plotImage().imshow(noise_1.array,ax[0],qs=None)
plotImage().imshow(noise_4.array,ax[1],qs=None)
In [ ]:
with np.errstate(divide='ignore'):
var_SR = 1./invvar # e/s
var_Astro = np.abs(img - const_sky) / expt # e/s
wt = 1./(var_SR + var_Astro) # s/e
# Zero out NaNs and masked pixels
wt[np.isfinite(wt) == False] = 0.
wt[dq != 0] = 0.
In [119]:
# How does the factor of 4 (which should be there), affect faint sources?
wcs= global_wcs(tim)
xpos,ypos,hw= 1040,3723,30
psf= galsim_psf(xpos,ypos,tim,wcs)
psf= Draw(psf,xpos,ypos,wcs)
gal= galsim_galaxy(xpos,ypos,tim,wcs,
n=4,rhalf=0.3696,mag=24,
e1=0.334,e2=0.273,angle=0)
gal= Draw(gal,xpos,ypos,wcs)
noise_1= noise_for_galaxy(gal,tim.zpscale,gain=1)
noise_4= noise_for_galaxy(gal,tim.zpscale,gain=4)
fig,ax= plt.subplots(1,3,figsize=(7,3))
plotImage().imshow(gal.array,ax[0],qs=None)
plotImage().imshow((gal+noise_4).array,ax[1],qs=None)
plotImage().imshow((gal+noise_1).array,ax[2],qs=None)
In [125]:
plot3d((gal+noise_4).array)
In [126]:
plot3d((gal+noise_1).array)
In [101]:
tim_galsim[olap].array.sum(),tim_wgal[olap].array.sum()-tim_galsim[olap].array.sum(),tim_wgal_wnoise[olap].array.sum()-tim_galsim[olap].array.sum()
Out[101]:
In [98]:
# g
In [ ]:
In [116]:
tim_galsim[zoom].array.sum(),gal[zoom].array.sum(),gal_wnoise[zoom].array.sum()
Out[116]:
In [ ]:
In [128]:
tim_galsim[olap].array.sum(),gal.array.sum(),gal_wnoise.array.sum(),(gal + gal_wnoise).array.sum()
Out[128]:
In [53]:
tim_wgal= tim_galsim.copy()
tim_wgal[olap] += gal[olap]
tim_wgal_wnoise= tim_galsim.copy()
tim_wgal_wnoise[olap] += gal_wnoise[olap]
fig,ax= plt.subplots(1,3,figsize=(6,3))
plotImage().imshow(tim_galsim[olap].array,ax[0])
plotImage().imshow(tim_wgal[olap].array,ax[1])
plotImage().imshow(tim_wgal_wnoise[olap].array,ax[2])
In [66]:
In [67]:
plot3d(tim_wgal[olap].array)
In [101]:
plot3d(tim_wgal_wnoise[olap].array,zlim=[0,3.5])
In [102]:
fig,ax=plt.subplots(1,3,figsize=(6,3))
plotImage().imshow(tim_galsim[olap].array,ax[0],qs=None)
plotImage().imshow(gal[olap].array,ax[1],qs=None)
wgal= tim_galsim.copy()
wgal[olap] += gal[olap]
plotImage().imshow(wgal[olap].array,ax[2],qs=None)
In [85]:
xpos
Out[85]:
In [94]:
fig,ax=plt.subplots(1,2,figsize=(4,3))
#tim_galsim.setCenter(xpos,ypos)
hw=30
box= galsim.BoundsI(xmin=xpos-hw,xmax=xpos+hw,ymin=ypos-hw,ymax=ypos+hw)
plotImage().imshow(tim_galsim[box].array,ax[0],qs=None)
In [96]:
gal.array.sum()
Out[96]:
In [61]:
fig,ax=plt.subplots(figsize=(3,3))
hw=15
box= galsim.BoundsI(xmin=xpos-hw,xmax=xpos+hw,ymin=ypos-hw,ymax=ypos+hw)
plotImage().imshow(tim_galsim[box].array,ax)
In [63]:
plot3d(tim_galsim[box].array)
In [64]:
tim_galsim[box].array.shape
Out[64]:
In [103]:
tim_galsim[olap].array.sum(), trac.flux_z, gal.added_flux, t.exptime,tim.zpscale
Out[103]:
In [67]:
tim_galsim[box].array.sum()/t.exptime, trac.flux_z, t.exptime
Out[67]:
In [135]:
hdu=fitsio.FITS('images/decam/DECam_CP/CP20160107/c4d_160116_084245_oki_z_v1.fits.fz')
cp=hdu['S17'].read()
cp.max(), tim.getImage().max()
Out[135]:
In [108]:
hdu=fitsio.FITS('images/decam/DECam_CP/CP20160107/c4d_160116_084245_oow_z_v1.fits.fz')
cp_ivar=hdu['S17'].read()
In [109]:
fig,ax= plt.subplots(1,2,figsize=(10,10))
plotImage().imshow(tim.getImage(),ax[0])
plotImage().imshow(cp,ax[1])
#psfim = self.psf.getPointSourcePatch(self.xpos, self.ypos).getImage()
In [110]:
slc= slice(1000,1500)
img_ratio= (cp/tim.getImage())[slc,slc].flatten()
ivar_ratio= (tim.getInvvar()/cp_ivar)[slc,slc].flatten()
set(img_ratio),set(ivar_ratio)
Out[110]:
In [111]:
tim.zpscale
Out[111]:
In [118]:
zp= t.ccdzpt + 2.5 * np.log10(t.exptime)
zpscale= 10**(zp - 22.5)/2.5
zpscale
Out[118]:
In [115]:
t.zpt
Out[115]:
In [88]:
img_ratio= f
ivar_ratio= f**2
print(img_ratio,ivar_ratio)
In [78]:
t.ccdzpt, t.ccdzpt + 2.5 * np.log10(t.exptime)
Out[78]:
In [79]:
im.ccdzpt
Out[79]: