Travis continuation integration Test Cases



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]:

    
# Testcase functions not in py/obiwan/ directory
a=os.path.join(os.getcwd(),
              "../../tests/end_to_end")
!ls $a

import sys
sys.path.append(a)









    



__init__.py				out_testcase_decals_z_dr5_elg_onedge
legacypipedir_1238p245_dataset_DR3	out_testcase_decals_z_dr5_star
legacypipedir_1238p245_dataset_DR5	outdir_1238p245_DR3
out_testcase_bassmzls_grz_dr6_elg	outdir_1238p245_DR5
out_testcase_decals_grz_dr5_elg		test_datasets.py
out_testcase_decals_z_dr5_elg		testcase_bassmzls_grz
out_testcase_decals_z_dr5_elg_allblobs	testcase_decals_grz
out_testcase_decals_z_dr5_elg_coadds	testcase_decals_z



In [5]:

    
print('hello from docker image')









    



hello from docker image

testcase_DR_z (z-band, 200x200 pix region, real galaxy at center)

all_blobs == True (fit model to all detected sources)



In [17]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg_allblobs/elg



In [13]:

    
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[ireal].bx,t.obitractor[ireal].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')



In [18]:

    
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('simcat id',t.simcat.id)
print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('img_apflux:',img_apflux)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])

print(sky_apflux + sims_apflux - t.obitractor[itrac].apflux_z[:,5])
print(t.obitractor[itrac].flux_z - t.simcat.zflux)









    



('simcat id', array([1, 2, 3, 4]))
('sky_apflux:', array([-0.8440134 , -0.56180657,  0.44062445, -1.00573623]))
('sims_apflux:', array([ 36.60767405,  36.59797994,  36.61960024,  36.56948628]))
('simcat_flux:', array([ 38.28181781,  38.28181781,  38.28181781,  38.28181781]))
('obitrac_flux:', array([ 33.39308548,  32.36734772,  34.74517059,  32.90026855], dtype=float32))
('img_apflux:', array([ 35.76366065,  36.03617336,  37.06022469,  35.56375005]))
('obitrac_apflux:', array([ 35.9729805 ,  36.11907196,  36.9563179 ,  35.5991478 ], dtype=float32))
[-0.20931985 -0.0828986   0.10390679 -0.03539775]
[-4.88873233 -5.91447009 -3.53664721 -5.38154925]



In [24]:

    
np.max([t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapeexp_r],axis=1)









    Out[24]:





array([ 0.44706032,  0.44706032], dtype=float32)



In [19]:

    
t.obitractor[ireal].apflux_z[:,5], t.obitractor[ireal].flux_z









    Out[19]:





(array([ 38.94617462], dtype=float32), array([ 38.28445053], dtype=float32))



In [20]:

    
t.simcat.rhalf,t.simcat.e1,t.simcat.e2,t.simcat.n









    Out[20]:





(array([ 0.37005,  0.37005,  0.37005,  0.37005]),
 array([ 0.33340144,  0.33340144,  0.33340144,  0.33340144]),
 array([ 0.27338931,  0.27338931,  0.27338931,  0.27338931]),
 array([4, 4, 4, 4]))



In [21]:

    
t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r









    Out[21]:





(array(['REX ', 'REX ', 'REX ', 'REX '],
       dtype='|S4'),
 array([ 0.44169348,  0.38656846,  0.44706032,  0.40121055], dtype=float32))



In [22]:

    
t.obitractor[itrac].shapeexp_r-t.simcat.rhalf









    Out[22]:





array([ 0.07164348,  0.01651846,  0.07701032,  0.03116055])



In [11]:

    
t.obitractor[itrac].shapeexp_r-t.simcat.rhalf









    Out[11]:





array([ 0.0719691 ,  0.01574509,  0.07740407,  0.03068192])



In [6]:

    
def flux2mag(flux):
    return -2.5*np.log10(1e-9 * flux)

def mag2flux(mag):
    return 10**(9-mag/2.5)

t.obitractor[ireal].type,t.obitractor[ireal].flux_z, flux2mag(t.obitractor[ireal].flux_z)









    Out[6]:





(array(['DEV '],
       dtype='|S4'),
 array([ 38.28407288], dtype=float32),
 array([ 18.54245377], dtype=float32))



In [50]:

    
t.obitractor[ireal].shapedev_r,t.obitractor[ireal].shapedev_e1,t.obitractor[ireal].shapedev_e2









    Out[50]:





(array([ 0.37025318], dtype=float32),
 array([ 0.33367831], dtype=float32),
 array([ 0.27334914], dtype=float32))



In [51]:

    
print(t.obitractor[itrac].type,t.obitractor[itrac].flux_z, flux2mag(t.obitractor[itrac].flux_z))
t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapeexp_e2









    



(array(['REX ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 33.38793564,  32.36148453,  34.74031448,  32.89451981], dtype=float32), array([ 18.69102669,  18.72492981,  18.64791489,  18.70719147], dtype=float32))






    Out[51]:





(array([ 0.44197041,  0.38681737,  0.44734058,  0.40145969], dtype=float32),
 array([ 0.,  0.,  0.,  0.], dtype=float32),
 array([ 0.,  0.,  0.,  0.], dtype=float32))



In [31]:

    
x=np.array([198,2,100,150]) + t.zoom[0]
y=np.array([2,198,2,198]) + t.zoom[2]
t.brickwcs.pixelxy2radec(x,y)









    Out[31]:





(array([ 174.23933751,  174.25500639,  174.24716429,  174.24318502]),
 array([ 24.32087412,  24.33512531,  24.32086773,  24.33513542]))



In [47]:

    
t.simcat.ra









    Out[47]:





array([ 174.24714385,  174.24715649,  174.25116781,  174.24313254])

all_blobs (coadds, top to bottom: image, simulated image, and residual)



In [49]:

    
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])
    
    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')
    
    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])



In [12]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)



In [15]:

    
plot3d(t.ivar_fits['z'])



In [50]:

    
len(isim),len(itrac),len(ireal)









    Out[50]:





(4, 4, 1)



In [51]:

    
t.simcat.zflux, t.obitractor[itrac].flux_z









    Out[51]:





(array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]),
 array([ 35.25556183,  34.26049423,  36.58589172,  34.72808075], dtype=float32))



In [52]:

    
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])


# apers= photutils.CircularAperture(([100],[100]), 
#                                    r=3.5/brickwcs.pixel_scale()) 
# apy_table = photutils.aperture_photometry(img_fits['g'], apers)
# real_apflux= np.array(apy_table['aperture_sum'])

# print(simcat.gflux,
#       sim_apflux,
#       real_apflux)
# print(flux2mag(simcat.gflux) - flux2mag(sim_apflux),flux2mag(simcat.gflux[0]),flux2mag(real_apflux))









    



('sky_apflux:', array([-0.84267303, -0.56008409,  0.44084312, -1.00828833]))
('sims_apflux:', array([ 38.55891982,  38.5594253 ,  38.56627197,  38.53363932]))
('simcat_flux:', array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]))
('obitrac_flux:', array([ 35.25556183,  34.26049423,  36.58589172,  34.72808075], dtype=float32))
('obitrac_apflux:', array([ 37.90312195,  38.06328583,  38.90242767,  37.55826569], dtype=float32))



In [53]:

    
print(sky_apflux + sims_apflux - t.obitractor[itrac].apflux_z[:,5])









    



[-0.18687516 -0.06394461  0.10468742 -0.03291469]



In [54]:

    
print(t.simcat.zflux + sky_apflux - t.obitractor[itrac].apflux_z[:,5])









    



[-0.48094632 -0.35852125 -0.19673589 -0.30170536]



In [55]:

    
apers= photutils.CircularAperture((t.obitractor[ireal].bx,t.obitractor[ireal].by), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])

print('real galaxy apflux:',img_apflux)
print('real galaxy obitrac flux:',t.obitractor[ireal].flux_z)
print('real galaxy obitrac apflux:',t.obitractor[ireal].apflux_z[0][5])
print(flux2mag(img_apflux) - flux2mag(t.obitractor[ireal].flux_z))









    



('real galaxy apflux:', array([ 38.91102173]))
('real galaxy obitrac flux:', array([ 38.26608276], dtype=float32))
('real galaxy obitrac apflux:', 38.91172)
[-0.01814652]



In [9]:

    
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 [ ]:

    
plot3d(img_fits)



In [ ]:

    
plot3d(sims_fits)



In [24]:

    
plot3d(img_fits - sims_fits)

Early coadds (just make the coadds, don't fit models)



In [4]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             all_blobs=True,early_coadds=True)
            
t= AnalyzeTestcase(**kwargs)
t.load_outputs()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg_allblobs_coadds/elg/028/0285m165/rs0/



In [5]:

    
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])
    
    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')
    
    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])

stars (inject stars into the images)



In [60]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             obj='star',
             onedge=False,add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_star_allblobs/star/174/1741p242/rs0/



In [61]:

    
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[ireal].bx,t.obitractor[ireal].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')



In [62]:

    
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print(t.simcat.id)
print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])
added_flux= np.array([39.135,39.102,39.122,39.115])
print('added_flux/1.02',added_flux/1.02)









    



[1 2 3 4]
('sky_apflux:', array([-0.84315645, -0.56062535,  0.4417466 , -1.00764566]))
('sims_apflux:', array([ 37.83162486,  37.85237236,  37.85502317,  37.80914618]))
('simcat_flux:', array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]))
('obitrac_flux:', array([ 36.89432526,  37.22283554,  37.80667496,  37.30305481], dtype=float32))
('obitrac_apflux:', array([ 37.17383194,  37.24513245,  38.22380829,  36.844944  ], dtype=float32))
('added_flux/1.02', array([ 38.36764706,  38.33529412,  38.35490196,  38.34803922]))



In [8]:

    
t.obitractor[itrac].type









    Out[8]:





array(['PSF ', 'PSF ', 'PSF ', 'PSF '],
      dtype='|S4')



In [12]:

    
dflux= t.simcat.zflux + \
        sky_apflux - t.obitractor[itrac].apflux_z[:,5]
print('diff flux: ',dflux)
dmag= flux2mag(t.simcat.zflux + sky_apflux) - flux2mag(t.obitractor[itrac].apflux_z[:,5])
print('diff mag: ',dmag)









    



('diff flux: ', array([-0.61297537, -0.36756682, -0.36557727, -0.42651694]))
('diff mag: ', array([ 0.01764039,  0.01053395,  0.01020564,  0.01235859]))

onedge (inject galaxies on the edges of the images, how good is best-fit model?)



In [19]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             onedge=True,add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()

print(t.simcat.id[isim])
print([(x,y) for x,y in zip(t.simcat.x[isim],t.simcat.y[isim])])









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg_onedge/elg
rhalf [-0.01998445  0.13738838  0.02882487  0.11826912]
apflux [ 0.1496385   0.09753577 -0.00312189  0.01799367]
skyflux [-0.26486025 -0.06903909  1.73439526 -0.05648272]
modelflux [ 0.33289523 -0.195715    5.08362384  2.63030047]
passed
[1 2 3 4]
[(197.99994482662373, 2.0000151657400238), (1.9999776837685204, 198.00002486870926), (100.00003559541528, 2.0000627804956821), (150.00001340770677, 197.9999763324513)]



In [20]:

    
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[itrac].bx,t.obitractor[itrac].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')



In [18]:

    
rhalf= np.max((t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r),
              axis=0)
e1= np.max((t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapedev_e1),
              axis=0)
e2= np.max((t.obitractor[itrac].shapeexp_e2,t.obitractor[itrac].shapedev_e2),
              axis=0)
print(t.simcat[isim].id)
print(t.obitractor[itrac].type,rhalf,e1,e2)
print(t.simcat.n[isim],t.simcat.rhalf[isim],t.simcat.e1[isim],t.simcat.e2[isim])
print(t.simcat.zflux[isim] - t.obitractor.flux_z[itrac])









    



[1 2 3 4]
(array(['DEV ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 0.34959555,  0.50696838,  0.39840487,  0.48784912], dtype=float32), array([ 0.11206593,  0.        ,  0.        ,  0.        ], dtype=float32), array([ 0.0281052,  0.       ,  0.       ,  0.       ], dtype=float32))
(array([4, 4, 4, 4]), array([ 0.36958,  0.36958,  0.36958,  0.36958]), array([ 0.33375,  0.33375,  0.33375,  0.33375]), array([ 0.27345,  0.27345,  0.27345,  0.27345]))
[ 0.33289523 -0.195715    5.08362384  2.63030047]



In [24]:

    
import galsim
a=galsim.PositionD(10,10)
a.x









    Out[24]:





10.0



In [39]:

    
a=galsim.Image(np.zeros((10,10)))
a.shift(dx=0.5,dy=0.5)
a.setCenter(10,10)



In [40]:

    
a.trueCenter(),a.center()









    Out[40]:





(galsim.PositionD(x=9.5, y=9.5), galsim.PositionI(x=10, y=10))

all_blobs == False (only fit models inside blobs that have fake sources)



In [22]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg/elg

The central galaxy is real and was not modeled



In [25]:

    
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=7./0.25,color='y')
plotImage().circles(t.obitractor[itrac].bx,t.obitractor[itrac].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.25,color='m')



In [23]:

    
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])









    



('sky_apflux:', array([-0.8424493 , -0.56035691,  0.44068548, -1.00731097]))
('sims_apflux:', array([ 35.84701954,  35.87222211,  35.87172834,  35.8296502 ]))
('simcat_flux:', array([ 38.28181781,  38.28181781,  38.28181781,  38.28181781]))
('obitrac_flux:', array([ 32.69493866,  31.70311165,  34.05732346,  32.2202034 ], dtype=float32))
('obitrac_apflux:', array([ 35.21274948,  35.39452744,  36.20796204,  34.85829544], dtype=float32))



In [24]:

    
print(t.simcat.rhalf[isim],t.simcat.e1[isim],t.simcat.e2[isim])
print(t.obitractor.type[itrac],t.obitractor.shapeexp_r[itrac],t.obitractor.shapeexp_e1[itrac],t.obitractor.shapeexp_e2[itrac])









    



(array([ 0.37005,  0.37005,  0.37005,  0.37005]), array([ 0.33340144,  0.33340144,  0.33340144,  0.33340144]), array([ 0.27338931,  0.27338931,  0.27338931,  0.27338931]))
(array(['REX ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 0.4420191 ,  0.38579509,  0.44745407,  0.40073192], dtype=float32), array([ 0.,  0.,  0.,  0.], dtype=float32), array([ 0.,  0.,  0.,  0.], dtype=float32))



In [26]:

    
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])
    
    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')
    
    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])



In [11]:

    
plot3d(t.img_fits['z'])



In [ ]:

    
all_blobs == False



In [23]:

    
### No noise
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=False)
no= AnalyzeTestcase(**kwargs)
no.load_outputs()
isim,itrac,ireal= no.match_simcat_tractor()

print(no.obitractor[itrac].apflux_z)[:,5]
plot3d(no.sims_fits['z'])









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z/elg/174/1741p242/rs0/
[ 37.90312195  38.06328583  38.90242767  37.55826569]



In [24]:

    
### w/noise original formula
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=True,all_blobs=False)
noise= AnalyzeTestcase(**kwargs)
noise.load_outputs()
isim,itrac,ireal= noise.match_simcat_tractor()

print(noise.obitractor[itrac].apflux_z)[:,5]
plot3d(noise.sims_fits['z'])









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_addnoise/elg/174/1741p242/rs0/
[ 37.78057098  38.04665375  38.87749863  37.56063843]



In [40]:

    
### w/noise mine
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=True,all_blobs=False)
my= AnalyzeTestcase(**kwargs)
my.load_outputs()
isim,itrac,ireal= my.match_simcat_tractor()

print(my.obitractor[itrac].apflux_z)[:,5]
plot3d(my.img_fits['z'])









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_addnoise/elg/174/1741p242/rs0/
[ 37.59292221  37.92685318  38.91709137  37.45846176]



In [31]:

    
sns.distplot(np.random.randn(100,100).flatten())









    Out[31]:





<matplotlib.axes._subplots.AxesSubplot at 0x7fce62556150>



In [8]:

    
print(no.obitractor[itrac].apflux_z[:,5] - my.obitractor[itrac].apflux_z[:,5])
print(no.obitractor[itrac].apflux_z[:,5] - noise.obitractor[itrac].apflux_z[:,5])









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-566c0d0c5800> in <module>()
----> 1 print(no.obitractor[itrac].apflux_z[:,5] - my.obitractor[itrac].apflux_z[:,5])
      2 print(no.obitractor[itrac].apflux_z[:,5] - noise.obitractor[itrac].apflux_z[:,5])

NameError: name 'no' is not defined

5 galaxies but only 4 blobs, meaning only the 4 fake ones were modeled



In [62]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)

testcase_DR_grz (Same as previous testcase, but three bands (grz) and a different region)

allblobs == True



In [83]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             obj='elg',
             add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg_allblobs/elg/028/0285m165/rs0/

Many sources are modeled, not just the 4 fake galaxies near center



In [84]:

    
fig,ax=plt.subplots(4,3,figsize=(7,8))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])
    
    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')
    
    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
    plotImage().imshow(t.ivar_fits[b],ax[3,i])



In [ ]:



In [23]:

    
plot3d(t.ivar_fits['r'])



In [85]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)



In [71]:

    
len(isim),len(itrac),len(ireal)









    Out[71]:





(4, 4, 1)



In [72]:

    
t.simcat.zflux, t.obitractor[itrac].flux_z









    Out[72]:





(array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]),
 array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))



In [78]:

    
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
for band in t.bands:
    print('band=',band)
    apy_table = photutils.aperture_photometry(t.img_fits[b], apers)
    img_apflux= np.array(apy_table['aperture_sum'])
    apy_table = photutils.aperture_photometry(t.sims_fits[b], apers)
    sims_apflux= np.array(apy_table['aperture_sum'])

    sky_apflux= img_apflux - sims_apflux
    #sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

    simcat_flux= t.simcat.get(b+'flux')
    trac_flux= t.obitractor[itrac].get('flux_'+b)
    trac_apflux=t.obitractor[itrac].get('apflux_'+b)[:,5]
    
    print('sky_apflux:',sky_apflux)
    #print('sky_meas:',sky_meas)
    print('sims_apflux:',sims_apflux)
    print('simcat_flux:',simcat_flux)
    print('obitrac_flux:',trac_flux)
    print('obitrac_apflux:',trac_apflux)
    print(sky_apflux + sims_apflux - trac_apflux)
    print(simcat_flux + sky_apflux - trac_apflux)









    



('band=', 'g')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]
('band=', 'r')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]
('band=', 'z')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]



In [74]:

    
t.img_fits['g'].shape









    Out[74]:





(200, 200)



In [79]:

    
t.obitractor[ireal].type,t.obitractor[ireal].shapeexp_r,t.obitractor[ireal].shapedev_r,









    Out[79]:





(array(['COMP'],
       dtype='|S4'),
 array([ 2.76125026], dtype=float32),
 array([ 0.24888119], dtype=float32))



In [80]:

    
for mod in ['shapedev','shapeexp']: 
    print(t.obitractor[ireal].get(mod+'_e1'),t.obitractor[ireal].get(mod+'_e2'))









    



(array([ 0.22466888], dtype=float32), array([-0.1372889], dtype=float32))
(array([ 0.12009902], dtype=float32), array([ 0.2517896], dtype=float32))



In [81]:

    
for b in 'grz':
    print(b,flux2mag(t.obitractor[ireal].get('flux_'+b)))









    



('g', array([ 19.44617081], dtype=float32))
('r', array([ 18.69445229], dtype=float32))
('z', array([ 18.13280869], dtype=float32))



In [86]:

    
t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r,









    Out[86]:





(array(['EXP ', 'EXP ', 'DEV ', 'DEV '],
       dtype='|S4'),
 array([ 2.30459976,  2.32127619,  0.        ,  0.        ], dtype=float32),
 array([ 0.        ,  0.        ,  2.61173296,  2.21803832], dtype=float32))



In [87]:

    
t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapeexp_e2









    Out[87]:





(array([ 0.06568176,  0.06891665,  0.        ,  0.        ], dtype=float32),
 array([-0.13783246, -0.14058098,  0.        ,  0.        ], dtype=float32))



In [88]:

    
t.obitractor[itrac].shapedev_e1,t.obitractor[itrac].shapedev_e2









    Out[88]:





(array([ 0.        ,  0.        ,  0.11277404,  0.05932122], dtype=float32),
 array([ 0.        ,  0.        , -0.12263217, -0.13615023], dtype=float32))



In [90]:

    
for b in 'grz':
    print(b,flux2mag(t.obitractor[itrac].get('flux_'+b)))









    



('g', array([ 19.4939785 ,  19.47326469,  19.43455124,  19.55050659], dtype=float32))
('r', array([ 18.72727394,  18.7271843 ,  18.68118286,  18.78072548], dtype=float32))
('z', array([ 18.16558075,  18.18332863,  18.09988785,  18.22354126], dtype=float32))



In [66]:

    
apers= photutils.CircularAperture((t.obitractor[ireal].bx,t.obitractor[ireal].by), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])

print('real galaxy apflux:',img_apflux)
print('real galaxy obitrac flux:',t.obitractor[ireal].flux_z)
print('real galaxy obitrac apflux:',t.obitractor[ireal].apflux_z[0][5])









    



('real galaxy apflux:', array([ 40.51106904]))
('real galaxy obitrac flux:', array([ 56.64278793], dtype=float32))
('real galaxy obitrac apflux:', 40.513737)

allblobs == False



In [41]:

    
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             obj='elg',
             add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()









    



Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg/elg

blobs are more extended, so less sources modeled, but still more than the 4 fake galaxies



In [42]:

    
fig,ax=plt.subplots(4,3,figsize=(7,8))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])
    
    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')
    
    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
    plotImage().imshow(t.ivar_fits[b],ax[3,i])



In [44]:

    
print(t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r)
print(t.simcat[isim].rhalf)









    



(array(['EXP ', 'EXP ', 'DEV ', 'DEV '],
      dtype='|S4'), array([ 2.30231023,  2.31898761,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        ,  2.60479164,  2.21619177], dtype=float32))
[ 2.  2.  2.  2.]



In [45]:

    
print(t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapedev_e1)
print(t.simcat[isim].e1)









    



(array([ 0.06559182,  0.06887048,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        ,  0.11189053,  0.05938281], dtype=float32))
[ 0.12408947  0.12408947  0.12408947  0.12408947]



In [46]:

    
print(t.obitractor[itrac].shapeexp_e2,t.obitractor[itrac].shapedev_e2)
print(t.simcat[isim].e2)









    



(array([-0.13772453, -0.14047767,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        , -0.12285928, -0.13606296], dtype=float32))
[ 0.26107663  0.26107663  0.26107663  0.26107663]



In [104]:

    
for b in 'grz':
    print(b,flux2mag(t.obitractor[itrac].get('flux_'+b)))









    



('g', array([ 19.49403381,  19.47325134,  19.43470192,  19.55048943], dtype=float32))
('r', array([ 18.72727203,  18.72712708,  18.68135262,  18.78062248], dtype=float32))
('z', array([ 18.16565132,  18.183321  ,  18.10007477,  18.22348976], dtype=float32))



In [28]:

    
plot3d(t.ivar_fits['g'])

only blobs containing fake sources are modeled, but this is still many sources



In [115]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)



In [90]:

    
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y), 
                                   r=3.5/t.brickwcs.pixel_scale()) 
for band in t.bands:
    print('band=',band)
    apy_table = photutils.aperture_photometry(t.img_fits[b], apers)
    img_apflux= np.array(apy_table['aperture_sum'])
    apy_table = photutils.aperture_photometry(t.sims_fits[b], apers)
    sims_apflux= np.array(apy_table['aperture_sum'])

    sky_apflux= img_apflux - sims_apflux
    #sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

    simcat_flux= t.simcat.get(b+'flux')
    trac_flux= t.obitractor[itrac].get('flux_'+b)
    trac_apflux=t.obitractor[itrac].get('apflux_'+b)[:,5]
    
    print('sky_apflux:',sky_apflux)
    #print('sky_meas:',sky_meas)
    print('sims_apflux:',sims_apflux)
    print('simcat_flux:',simcat_flux)
    print('obitrac_flux:',trac_flux)
    print('obitrac_apflux:',trac_apflux)
    print(sky_apflux + sims_apflux - trac_apflux)
    print(simcat_flux + sky_apflux - trac_apflux)









    



('band=', 'g')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]
('band=', 'r')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]
('band=', 'z')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]

Measured mags are systematcially 0.1-0.2 mags brighter than truth mags.

Why? Is it due to add_flux being systematically larger than requested flux? something else?



In [137]:

    
a=np.ones((10,10))+34
a /= a.sum()
a.sum()









    Out[137]:





0.99999999999999989



In [148]:

    
assert(len(itrac) == 4)
print(simcat[isim].zflux-obitractor[itrac].flux_z)
print(np.sqrt(1/obitractor[itrac].flux_ivar_z))
print(flux2mag(simcat[isim].zflux)-flux2mag(obitractor[itrac].flux_z))









    



[ 2.80877871  3.7907352   1.49074359  3.23818774]
[ 0.36161846  0.34976208  0.36395597  0.35320479]
[-0.08277248 -0.11326717 -0.04314349 -0.09600375]



In [102]:

    
# sizes?
print(simcat[isim].rhalf-obitractor[itrac].shapeexp_r)









    



[-0.04129598  0.0056425  -0.05493272 -0.01868245]



In [103]:

    
# Does extinction vary by this much? I doubt it...
mags= flux2mag(obitractor[itrac].flux_z)
print('measured mag variation=',mags - np.median(mags))
mags= flux2mag(obitractor[itrac].flux_z/obitractor[itrac].mw_transmission_z)
print('extinction variation=',mags - np.median(mags))









    



('measured mag variation=', array([-0.00476646,  0.02602196, -0.05817604,  0.00476646], dtype=float32))
('extinction variation=', array([-0.00471115,  0.02607346, -0.05801773,  0.00471115], dtype=float32))



In [104]:

    
np.std(obitractor.mw_transmission_z)









    Out[104]:





6.3498112e-05

sim sources should have similar flux to real galaxy at center



In [105]:

    
# Diff b/w fake truth and measure real galaxy identical
print(simcat[isim].zflux-obitractor[ireal].flux_z)
print(flux2mag(simcat[isim].zflux)-flux2mag(obitractor[ireal].flux_z))









    



[-0.00122647 -0.00122647 -0.00122647 -0.00122647]
[  3.50646973e-05   3.50646973e-05   3.50646973e-05   3.50646973e-05]



In [106]:

    
# But diff b/w measure fake and measure real galaxy are large
print(obitractor[itrac].flux_z-obitractor[ireal].flux_z)
print(flux2mag(obitractor[itrac].flux_z)-flux2mag(obitractor[ireal].flux_z))









    



[-5.3830986  -6.30246735 -3.72509003 -5.67053986]
[ 0.164608    0.19539642  0.11119843  0.17414093]



In [110]:

    
obitractor.flux_z









    Out[110]:





array([ 32.88297653,  38.26607513,  32.59553528,  34.54098511,
        31.96360779,   4.96016264,   2.12047195], dtype=float32)



In [107]:



In [113]:

    
plot3d(img_jpg[:,:,0])



In [114]:

    
plot3d(model_jpg[:,:,0])

What "zoom" is needed for this testcase?



In [22]:

    
fig,ax=plt.subplots()
plotImage().imshow(blobs,ax,qs=None)



In [27]:

    
ra=[174.24714700910047, 174.24715333044438, 174.24915898792247, 174.2451413528916]
dec=[24.326224120319477, 24.32988191327332, 24.328050866811363, 24.32805515572893]
np.mean(ra),np.mean(dec)









    Out[27]:





(174.24715017008975, 24.328053014033273)



In [41]:

    
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

from legacypipe.survey import LegacySurveyData, wcs_for_brick
survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('1741p242')
brickwcs = wcs_for_brick(brickinfo)









    



---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-41-a294eb3f18c1> in <module>()
      4 survey = LegacySurveyData()
      5 brickinfo = survey.get_brick_by_name('1741p242')
----> 6 brickwcs = wcs_for_brick(brickinfo)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.pyc in wcs_for_brick(b, W, H, pixscale)
    533     from astrometry.util.util import Tan
    534     pixscale = pixscale / 3600.
--> 535     return Tan(b.ra, b.dec, W/2.+0.5, H/2.+0.5,
    536                -pixscale, 0., 0., pixscale,
    537                float(W), float(H))

AttributeError: 'NoneType' object has no attribute 'ra'



In [28]:

    
brickwcs.radec2pixelxy(174.24715,24.32805)









    Out[28]:





(True, 190.26994464884524, 2873.6879080454523)



In [38]:

    
_,x,y= brickwcs.radec2pixelxy(174.2471,24.3262)
x,y= int(x),int(y)+25
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(blobs,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)









    














    











    



(slice(90, 290, None), slice(2773, 2973, None))

Find som isolated sources with 3 bands, make 4 testcases 1 per band and 1 all-band

1) use viewer to find good source, write down brickname, ra,dec of source and 3 good ccds touching that ra,dec by clicking on the source in the viewer

2) grab the coadd/ccds file for that brick and cut to just those three ccds

3) create a initial dir structure by running create_testcase on NERSC using the 3-ccd ccd file. This is NOT the final product yet, but is need for #4.

4) You need to read each of the ccds as tim objects and print each of their x,y pixel positions for the ra,dec above. This requires the files output by #3 and the original images (on project at NERSC). scp -r the output dir from #3 to a dir on laptop. Keep the subdirs in images/ but remove the images you scp'd over. Copy the actual oo[idw] images from project to each band's dir. Then run the "Correct Region?" and "TIM object for each CCD" and "x,y pixels for each ccd" cells below. This uses the tim.wcs.wcs.radec2pixelxy func

6) replace the ccds table with these new values for ccd_x0,_x1,_y0,_y1 (round to nearest integer) then copy the ccds table back to NERSC and rerun create_testcase but using the new ccds table.

7) Done!

DECam



In [27]:

    
# DR3: Correct Region?
from legacypipe.survey import LegacySurveyData, wcs_for_brick
dr= os.path.join(os.environ['HOME'],
                'myrepo/legacypipe/dr3_testcase')
os.environ["LEGACY_SURVEY_DIR"]= dr

survey = LegacySurveyData()
brick='2357p087'
brickinfo = survey.get_brick_by_name(brick)
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 235.6799,8.8597
print(brickwcs.radec2pixelxy(ra,dec))

# Plot 
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
img_jpg= readImage(os.path.join(dr,'legacysurvey-%s-image.jpg' % brick),
                   jpeg=True)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(x,y)
print(x-hw,x+hw,y-hw,y+hw)

for x,y in [(x-hw,y),(x+hw,y),(x,y+hw),(x,y-hw)]:
    print(brickwcs.pixelxy2radec(x,y))









    



(True, 2675.942971854415, 3307.9059576884274)






    














    











    



(2675, 3307)
(2575, 2775, 3207, 3407)
(235.68733510594618, 8.859635266642558)
(235.67260383081032, 8.85963274643333)
(235.67996820762602, 8.866911823308769)
(235.67997072695223, 8.852356330322968)



In [15]:

    
# TIM objects for each CCD
from legacypipe.decam import DecamImage
from legacypipe.survey import LegacySurveyData

ccds= fits_table(os.path.join(dr,'dr3_testcase_ccds.fits'))

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)









    



Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_ooi_g_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_oow_g_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_ood_g_v2.fits.fz hdu 35
Warning: un-setting SATUR for 3 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00349/decam-00349605.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00349/00349605/decam-00349605-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00349/decam-00349605.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00349/00349605/decam-00349605-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.09, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.114103
Galaxy norm 0.0902848326709
Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_ooi_r_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_oow_r_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_ood_r_v2.fits.fz hdu 35
Warning: un-setting SATUR for 2 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00347/decam-00347323.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00347/00347323/decam-00347323-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00347/decam-00347323.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00347/00347323/decam-00347323-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.93, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.136238
Galaxy norm 0.102938388518
Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_ooi_z_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_oow_z_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_ood_z_v2.fits.fz hdu 35
Warning: un-setting SATUR for 47 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00346/decam-00346606.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00346/00346606/decam-00346606-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00346/decam-00346606.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00346/00346606/decam-00346606-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.42, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.106902
Galaxy norm 0.0885711598239



In [24]:

    
# x,y pixels for each ccd
hw=100
for b in ccds.filter:
    _,x,y= tims[b].wcs.wcs.radec2pixelxy(ra,dec)
    print('%s: ' % b,x-hw,x+hw,y-hw,y+hw)









    



('g: ', 318.70780228177273, 518.7078022817727, 1921.286605631049, 2121.286605631049)
('r: ', 367.16452653880515, 567.1645265388051, 1968.49883705385, 2168.49883705385)
('z: ', 357.70981507257636, 557.7098150725764, 1941.7682361872326, 2141.7682361872326)



In [43]:

    
# brick 0285m165, RA,Dec = (28.4194, -16.4362) in viewer looks like good 3 colors source
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz/'))
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('0285m165')
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 28.4194, -16.4362
brickwcs.radec2pixelxy(ra,dec)









    Out[43]:





(True, 3177.4193025452664, 2676.786370753425)



In [67]:

    
brick='0285m165'
DATA_DIR= os.path.join(os.environ['HOME'],
                       'mydata',brick,'dr5')
dr5tractor= fits_table(os.path.join(DATA_DIR,
                      'tractor-%s.fits' % brick))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'dr5_3ccds.fits'))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'ccds_0285m165_dr5coadddir.fits'))
ccds= fits_table(os.path.join(IN_DIR,'survey-ccds-1.fits.gz'))

img_jpg= readImage(os.path.join(DATA_DIR,'legacysurvey-0285m165-image.jpg'),
                   jpeg=True)



In [46]:

    
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)









    














    











    



(slice(3077, 3277, None), slice(2576, 2776, None))

Mosaic, 90Prime



In [4]:

    
# DR6: Correct Region?
from legacypipe.survey import LegacySurveyData, wcs_for_brick
dr= os.path.join(os.environ['HOME'],
                'myrepo/obiwan/tests/end_to_end/testcase_bassmzls_grz')
os.environ["LEGACY_SURVEY_DIR"]= dr

survey = LegacySurveyData()
brick='2176p330'
brickinfo = survey.get_brick_by_name(brick)
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 217.5429, 33.0873
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
print('brick box pixels: x=[%d,%d],y=[%d,%d]' % \
     (x-hw,x+hw,y-hw,y+hw))

for x,y in [(x-hw,y),(x+hw,y),(x,y+hw),(x,y-hw)]:
    print(brickwcs.pixelxy2radec(x,y))

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
img_jpg= readImage(os.path.join(dr,'legacysurvey-%s-image.jpg' % brick),
                   jpeg=True)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)









    



---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-4-4d057a932720> in <module>()
      7 survey = LegacySurveyData()
      8 brick='2176p330'
----> 9 brickinfo = survey.get_brick_by_name(brick)
     10 brickwcs = wcs_for_brick(brickinfo)
     11 

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_brick_by_name(self, brickname)
   1161         Returns a brick (as one row in a table) by name (string).
   1162         '''
-> 1163         B = self.get_bricks_readonly()
   1164         I, = np.nonzero(np.array([n == brickname for n in B.brickname]))
   1165         if len(I) == 0:

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_bricks_readonly(self)
   1140         '''
   1141         if self.bricks is None:
-> 1142             self.bricks = self.get_bricks()
   1143             # Assert that bricks are the sizes we think they are.
   1144             # ... except for the two poles, which are half-sized

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_bricks(self)
   1133         uses a cached version.
   1134         '''
-> 1135         return fits_table(self.find_file('bricks'))
   1136 
   1137     def get_bricks_readonly(self):

~/myinstall/astrometry/lib/python/astrometry/util/fits.py in fits_table(dataorfn, rows, hdunum, hdu, ext, header, columns, column_map, lower, mmap, normalize, use_fitsio, tabledata_class)
    665 
    666         if fitsio:
--> 667             F = fitsio.FITS(dataorfn)
    668             data = F[hdunum]
    669             hdr = data.read_header()

/srv/py3_venv/lib/python3.5/site-packages/fitsio/fitslib.py in __init__(self, filename, mode, **keys)
    361                     create=1
    362 
--> 363         self._FITS =  _fitsio_wrap.FITS(filename, self.intmode, create)
    364 
    365     def close(self):

OSError: FITSIO status = 104: could not open the named file
failed to find or open the following file: (ffopen)
/root/myrepo/obiwan/tests/end_to_end/testcase_bassmzls_grz/survey-bricks.fits.gz



In [10]:

    
survey = LegacySurveyData(ccds=ccds)
im = survey.get_image_object(ccds[0])

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)
im.get_tractor_image(**kwargs)









    



camera=-mosaic- camera.strip=-mosaic-
Reading image slice: None
Reading image from /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ooi_zd_v1-CCD3.fits hdu 1
Reading weight map image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_oow_zd_v1-CCD3.fits ext 1
Reading data quality image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ood_zd_v1-CCD3.fits ext 1
Remapping weight map for 00110480-CCD3
Merged spline sky model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky-merged/00110/mosaic-00110480.fits
Reading sky model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky/00110/00110480/mosaic-00110480-CCD3.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex-merged/00110/mosaic-00110480.fits
Reading PsfEx model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex/00110/00110480/mosaic-00110480-CCD3.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.81, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.145167






    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-b9337c9d9c8d> in <module>()
      4 kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
      5               pixels=True, dq=True, invvar=True)
----> 6 im.get_tractor_image(**kwargs)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in get_tractor_image(self, slc, radecpoly, gaussPsf, pixPsf, hybridPsf, splinesky, nanomaggies, subsky, tiny, dq, invvar, pixels, constant_invvar)
    331         print('PSF norm:', tim.psfnorm)
    332         # Galaxy-detection norm
--> 333         tim.galnorm = self.galaxy_norm(tim)
    334         print('Galaxy norm', tim.galnorm)
    335         assert(tim.galnorm < tim.psfnorm)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in galaxy_norm(self, tim, x, y)
    437             y = h//2
    438         pos = tim.wcs.pixelToPosition(x, y)
--> 439         gal = SimpleGalaxy(pos, NanoMaggies(**{band:1.}))
    440         S = 32
    441         mm = ModelMask(int(x-S), int(y-S), 2*S+1, 2*S+1)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.py in __init__(self, *args)
     91 
     92     def __init__(self, *args):
---> 93         super(SimpleGalaxy, self).__init__(*args)
     94         self.shape = SimpleGalaxy.shape
     95 

TypeError: super(type, obj): obj must be an instance or subtype of type



In [11]:

    
# TIM objects for each CCD
from legacypipe.mosaic import MosaicImage
from legacypipe.bok import BokImage
from legacypipe.survey import LegacySurveyData

ccds= fits_table(os.path.join(dr,'survey-ccds-1.fits.gz'))

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)









    



camera=-mosaic- camera.strip=-mosaic-
Reading image slice: None
Reading image from /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ooi_zd_v1-CCD3.fits hdu 1
Reading weight map image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_oow_zd_v1-CCD3.fits ext 1
Reading data quality image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ood_zd_v1-CCD3.fits ext 1
Remapping weight map for 00110480-CCD3
Merged spline sky model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky-merged/00110/mosaic-00110480.fits
Reading sky model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky/00110/00110480/mosaic-00110480-CCD3.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex-merged/00110/mosaic-00110480.fits
Reading PsfEx model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex/00110/00110480/mosaic-00110480-CCD3.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.81, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.145167






    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-712d6a45f37e> in <module>()
     16 for ccd in ccds:
     17     im = survey.get_image_object(ccd)
---> 18     tims[ccd.filter] = im.get_tractor_image(**kwargs)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in get_tractor_image(self, slc, radecpoly, gaussPsf, pixPsf, hybridPsf, splinesky, nanomaggies, subsky, tiny, dq, invvar, pixels, constant_invvar)
    331         print('PSF norm:', tim.psfnorm)
    332         # Galaxy-detection norm
--> 333         tim.galnorm = self.galaxy_norm(tim)
    334         print('Galaxy norm', tim.galnorm)
    335         assert(tim.galnorm < tim.psfnorm)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in galaxy_norm(self, tim, x, y)
    437             y = h//2
    438         pos = tim.wcs.pixelToPosition(x, y)
--> 439         gal = SimpleGalaxy(pos, NanoMaggies(**{band:1.}))
    440         S = 32
    441         mm = ModelMask(int(x-S), int(y-S), 2*S+1, 2*S+1)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.py in __init__(self, *args)
     91 
     92     def __init__(self, *args):
---> 93         super(SimpleGalaxy, self).__init__(*args)
     94         self.shape = SimpleGalaxy.shape
     95 

TypeError: super(type, obj): obj must be an instance or subtype of type



In [24]:

    
# x,y pixels for each ccd
hw=100
for b in ccds.filter:
    _,x,y= tims[b].wcs.wcs.radec2pixelxy(ra,dec)
    print('%s: ' % b,x-hw,x+hw,y-hw,y+hw)









    



('g: ', 318.70780228177273, 518.7078022817727, 1921.286605631049, 2121.286605631049)
('r: ', 367.16452653880515, 567.1645265388051, 1968.49883705385, 2168.49883705385)
('z: ', 357.70981507257636, 557.7098150725764, 1941.7682361872326, 2141.7682361872326)



In [43]:

    
# brick 0285m165, RA,Dec = (28.4194, -16.4362) in viewer looks like good 3 colors source
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz/'))
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('0285m165')
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 28.4194, -16.4362
brickwcs.radec2pixelxy(ra,dec)









    Out[43]:





(True, 3177.4193025452664, 2676.786370753425)



In [67]:

    
brick='0285m165'
DATA_DIR= os.path.join(os.environ['HOME'],
                       'mydata',brick,'dr5')
dr5tractor= fits_table(os.path.join(DATA_DIR,
                      'tractor-%s.fits' % brick))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'dr5_3ccds.fits'))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'ccds_0285m165_dr5coadddir.fits'))
ccds= fits_table(os.path.join(IN_DIR,'survey-ccds-1.fits.gz'))

img_jpg= readImage(os.path.join(DATA_DIR,'legacysurvey-0285m165-image.jpg'),
                   jpeg=True)



In [46]:

    
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)









    














    











    



(slice(3077, 3277, None), slice(2576, 2776, None))

DECam



In [16]:

    
from legacypipe.decam import DecamImage
from legacypipe.survey import LegacySurveyData



In [80]:

    
kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)









    



Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ooi_z_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_oow_z_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ood_z_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00384/decam-00384128.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00384/00384128/decam-00384128-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00384/decam-00384128.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00384/00384128/decam-00384128-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.64, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.158392
Galaxy norm 0.113801322034
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ooi_r_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_oow_r_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ood_r_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390924.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390924/decam-00390924-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390924.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390924/decam-00390924-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.49, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.178092
Galaxy norm 0.121385678047
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ooi_g_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_oow_g_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ood_g_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390498.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390498/decam-00390498-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390498.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390498/decam-00390498-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.05, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.123636
Galaxy norm 0.0966921909813



In [81]:

    
tim=tims['g']
tim.wcs.wcs.radec2pixelxy(ra,dec)









    Out[81]:





(True, 918.0465384893796, 420.42404627704695)



In [82]:

    
_,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(1,3,figsize=(7,4))
for i,band in enumerate(['g','r','z']):
    plotImage().imshow(sliceImage(tims[band].getImage(),xslice=xslc,yslice=yslc).T,
                       ax[i])



In [87]:

    
def flux2mag(flux):
    return -2.5*np.log10(1e-9 * flux)

hw=30
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(1,3,figsize=(7,4))
for i,band in enumerate(['g','r','z']):
    img= sliceImage(tims[band].getImage(),xslice=xslc,yslice=yslc).T
    plotImage().imshow(img,ax[i])
    # Extract the nanomaggie flux for each band that we will add to the images
    print(band,'nanomags=',img.sum(),'mags=',flux2mag(img.sum()))









    














    











    



('g', 'nanomags=', 17.956371, 'mags=', 19.364453556756388)
('r', 'nanomags=', 32.756645, 'mags=', 18.711751458105404)
('z', 'nanomags=', 46.451378, 'mags=', 18.332503497939175)



In [72]:

    
# flux of this galaxy
from astrometry.libkd.spherematch import match_radec
_,itrac,d= match_radec(np.array([ra]), np.array([dec]), dr5tractor.ra, dr5tractor.dec,          
                        5./3600.0,nearest=True)
itrac









    Out[72]:





array([4962], dtype=int32)



In [74]:

    
trac=dr5tractor.copy()
trac.cut(itrac)
len(trac)









    Out[74]:





1



In [76]:

    
trac.ra,trac.dec, ra,dec









    Out[76]:





(array([ 28.4192355]), array([-16.43595489]), 28.4194, -16.4362)



In [88]:

    
# Tractor measurements on many imags should be close to above values for the three images
for band in ['g','r','z']:
    flux=trac.get('flux_'+band)
    print(band,'nanomags=',flux,'mags=',flux2mag(flux))









    



('g', 'nanomags=', array([ 16.06177902], dtype=float32), 'mags=', array([ 19.48551559], dtype=float32))
('r', 'nanomags=', array([ 31.8071003], dtype=float32), 'mags=', array([ 18.74369049], dtype=float32))
('z', 'nanomags=', array([ 52.83349609], dtype=float32), 'mags=', array([ 18.19272614], dtype=float32))



In [89]:

    
trac.ra,trac.dec,trac.type,flux2mag(trac.flux_g),flux2mag(trac.flux_r),flux2mag(trac.flux_z)









    Out[89]:





(array([ 28.4192355]), array([-16.43595489]), array(['EXP '],
       dtype='|S4'), array([ 19.48551559], dtype=float32), array([ 18.74369049], dtype=float32), array([ 18.19272614], dtype=float32))



In [90]:

    
trac.shapeexp_r,trac.shapeexp_e1,trac.shapeexp_e2









    Out[90]:





(array([ 2.4597683], dtype=float32),
 array([ 0.12245165], dtype=float32),
 array([ 0.23880503], dtype=float32))



In [104]:

    
# Lets circle ra,dec where we'll add galaxies to the r-band image
print(ra,dec)
tim= tims['r']
dd= 60 * 0.262 /3600 # pix*as/pix converted to deg
newra,newdec=[],[]
for dx,dy in zip([dd,-dd,0,0],[0,0,dd,-dd]):
    #newx,newy=x+dx,y+dy
    #ra,dec= tim.wcs.wcs.pixelxy2radec(newx,newy)
    #print(newx,newy,ra,dec)
    newra.append(ra + dx)
    newdec.append(dec + dy)
print(newra)
print(newdec)
_,newx,newy= tim.wcs.wcs.radec2pixelxy(newra,newdec)
print(newx)
print(newy)

_,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(sliceImage(tims['r'].getImage(),xslice=xslc,yslice=yslc),
                   ax)
# cirlces
plotImage().circles(newx,newy,ax,
                    xslice=xslc,yslice=yslc,
                    r_pixels=5./0.262,color='y')
print('IMAGE SLICE:',xslc,yslc)









    



(28.4194, -16.4362)
[28.423766666666666, 28.415033333333334, 28.4194, 28.4194]
[-16.4362, -16.4362, -16.431833333333334, -16.440566666666665]
[ 918.6980714   918.75985247  859.09916007  978.3570594 ]
[ 477.86830435  363.5315702   420.66734678  420.73262761]






    














    











    



('IMAGE SLICE:', slice(818, 1018, None), slice(320, 520, None))

Now find ccd_x0, ..., ccd_y1 for each ccd



In [125]:

    
# start with coadd/ccds
ccds= fits_table(os.path.join(DATA_DIR,
                 'ccds_0285m165_dr5coadddir.fits'))



In [111]:

    
kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=False, invvar=False)


survey = LegacySurveyData(ccds=ccds)
vals=[]
for band in ['g','r','z']:
    t= ccds[ccds.filter == band]
    for ccd in t:
        im = survey.get_image_object(ccd)
        tim = im.get_tractor_image(**kwargs)
        _,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
        print(band,int(x),int(y))
        vals.append( (band,int(x),int(y)) )









    



Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ooi_g_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390498.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390498/decam-00390498-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.00575166300675
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390498.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390498/decam-00390498-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.10, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.123362
Galaxy norm 0.0972685759985
('g', 1814, 3205)
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ooi_r_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390924.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390924/decam-00390924-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.00801456716987
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390924.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390924/decam-00390924-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.41, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.162455
Galaxy norm 0.120581835731
('r', 1783, 3215)
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ooi_z_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00384/decam-00384128.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00384/00384128/decam-00384128-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.0346027164371
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00384/decam-00384128.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00384/00384128/decam-00384128-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.53, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.137621
Galaxy norm 0.111473202692
('z', 1782, 3263)
[('g', 1814, 3205), ('r', 1783, 3215), ('z', 1782, 3263)]



In [114]:

    
print(vals)
dx=100
from collections import defaultdict
x0,x1,y0,y1= defaultdict(list),defaultdict(list),defaultdict(list),defaultdict(list)
for b,xc,yc in vals:
    x0[b] += [=xc-100,xc+100,yc-100,yc+100)









    



[('g', 1814, 3205), ('r', 1783, 3215), ('z', 1782, 3263)]
('g', 1714, 1914, 3105, 3305)
('r', 1683, 1883, 3115, 3315)
('z', 1682, 1882, 3163, 3363)



In [122]:

    
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz_200x200/')
g=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141223_020302_ooi_g_a1-S10.fits'))[1].read()
r=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141224_033245_ooi_r_a1-S10.fits'))[1].read()
z=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141206_043719_ooi_z_a1-S10.fits'))[1].read()
g.shape









    Out[122]:





(200, 200)



In [123]:

    
fig,ax=plt.subplots(1,3,figsize=(6,3))
plotImage().imshow(g,ax[0])
plotImage().imshow(r,ax[1])
plotImage().imshow(z,ax[2])

It worked!

6'' separation + testcase_DR_z



In [5]:

    
OUT_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'out_testcase_DR5_z/elg/174/1741p242/rs0/')
IN_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'testcase_DR5_z/')

simcat= fits_table(os.path.join(OUT_DIR,'obiwan/simcat-elg-1741p242.fits'))
obitractor= fits_table(os.path.join(OUT_DIR,'tractor/tractor-1741p242.fits'))

blobs= fitsio.FITS(os.path.join(OUT_DIR,'metrics/blobs-1741p242.fits.gz'))[0].read()

img_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-image.jpg'),
                   jpeg=True)
model_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-model.jpg'),
                     jpeg=True)
resid_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-resid.jpg'),
                     jpeg=True)



In [6]:

    
len(simcat),len(obitractor),blobs.shape









    Out[6]:





(4, 3, (200, 200))

allblobs == False

Bug in Tractor, 3/5 sources modeled because they are apparently too close for the algorithm!



In [7]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(blobs,ax[0,0],qs=None)
plotImage().imshow(img_jpg,ax[0,1],qs=None)
plotImage().imshow(model_jpg,ax[1,0],qs=None)
plotImage().imshow(resid_jpg,ax[1,1],qs=None)

allblobs == True

Same Bug in Tractor, even when allblobs (there is only 1 blob anyway)



In [8]:

    
OUT_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'out_testcase_DR5_z_allblobs/elg/174/1741p242/rs0/')
IN_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'testcase_DR5_z_allblobs/')

simcat= fits_table(os.path.join(OUT_DIR,'obiwan/simcat-elg-1741p242.fits'))
obitractor= fits_table(os.path.join(OUT_DIR,'tractor/tractor-1741p242.fits'))

blobs= fitsio.FITS(os.path.join(OUT_DIR,'metrics/blobs-1741p242.fits.gz'))[0].read()

img_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-image.jpg'),
                   jpeg=True)
model_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-model.jpg'),
                     jpeg=True)
resid_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-resid.jpg'),
                     jpeg=True)



In [9]:

    
len(simcat),len(obitractor),blobs.shape









    Out[9]:





(4, 3, (200, 200))



In [10]:

    
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(blobs,ax[0,0],qs=None)
plotImage().imshow(img_jpg,ax[0,1],qs=None)
plotImage().imshow(model_jpg,ax[1,0],qs=None)
plotImage().imshow(resid_jpg,ax[1,1],qs=None)