Thanks to Ani Thakar for SQL
Request: improve the astronomy!
# standard first block for defining the token and makinhg it available as a system variable for the session
# token must be replaced with new one once it has expired
with open('/home/idies/keystone.token', 'r') as f:
token = f.read().rstrip('\n')
import SciServer.CasJobs
import pandas
import tables
import numpy as np
import astropy
from astropy.io import fits
from astropy import wcs
import skimage.io
import urllib
import os
import matplotlib
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
For doing some scaling of FITS images we use functions from the img_scale.py python file, found at http://dept.astro.lsa.umich.edu/~msshin/science/code/Python_fits_image/img_scale.py. Updated to python 3 (by putting parentheses in print statements).
# python scripts can be uploaded to container and imported
# downloaded from
import img_scale
sql="""
SELECT a.objid as head, c.objid2 as match, b.matchcount,
p.fieldid as head_field, d.fieldid as match_field,
dbo.fGetUrlFitsCFrame(d.fieldid, 'g') as fits_g,
dbo.fGetUrlFitsCFrame(d.fieldid, 'r') as fits_r,
dbo.fGetUrlFitsCFrame(d.fieldid, 'z') as fits_z,
p.ra, d.ra as match_ra, p.dec, d.dec as match_dec
, p.petror90_r
from (select top 1 * from galaxy where objId=8658194378960928809) a
join matchhead b on a.objid=b.objid -- join with matchhead
join photoobj p on a.objid=p.objid -- get matchhead photoobj
join match c on c.objid1=b.objid -- join with all the matches
join photoobjall d on c.objid2=d.objid -- get match photoobj
order by d.fieldid -- order by matchhead objid
"""
queryResponse = SciServer.CasJobs.executeQuery(sql, "Stripe82",token=token)
obss = pandas.read_csv(queryResponse,index_col=None)
# show first 10 rows
obss[:10]
def cutout(ff,ra,dec,hw):
hdulist = fits.open(ff)
w = wcs.WCS(hdulist[0].header)
crd = np.array([[ra,dec]], np.float_)
pixcoords = np.around(list(w.wcs_world2pix(crd,1)))
sh=hdulist[0].data.shape
xfrom=max(0,pixcoords[0][1]-hw)
xto=min(pixcoords[0][1]+hw-1,sh[0])
yfrom=max(pixcoords[0][0]-hw,0)
yto=min(pixcoords[0][0]+hw-1,sh[1])
scidata = hdulist[0].data[xfrom:xto,yfrom:yto]
_img=np.zeros((hw*2-1,hw*2-1))
_mask=np.zeros(_img.shape)
x0=hw-(pixcoords[0][1]-xfrom)
x1=hw+xto-pixcoords[0][1]
y0=hw-(pixcoords[0][0]-yfrom)
y1=hw+yto-pixcoords[0][0]
_img[x0:x1,y0:y1]=scidata
_mask[x0:x1,y0:y1]=1
_img=np.fliplr(_img.T)
_mask=np.fliplr(_mask.T)
return _img,_mask
width=200
height=200
pixelsize=0.396
_cmap='afmhot'
#_cmap='cool'
_vmin=1000
_vmax=10000
nx=4
ny=2
nmax=nx*ny
plt.figure(figsize=(15, 15*nx/ny))
sp = 1
_min=np.inf
_max=-np.inf
for ix in range(0,nmax-2):
if ix == 0:
ra=obss.ra[ix]
dec=obss.dec[ix]
r90=obss.petror90_r[ix]
size=2*r90/pixelsize
scale=size/width
url="http://skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx?ra="+str(ra)
url+="&dec="+str(dec)+"&scale="""+str(scale)+"&width="+str(width)
url+="&height="+str(height)
fits_r=obss.fits_r[ix]
# IFF sdss_das is mounted read from file
# NOTE relative path should be changed to conform to user's folder structure
# OR use absolute /home/idies/workspace/sdss_das/das2
if os.path.isdir('../../sdss_das/'):
fits_r=fits_r.replace('http://das.sdss.org','../../sdss_das/das2')
print ('downloading "local" file',fits_r)
else:
print('downloading remote file',fits_r)
fimg,_mask=cutout(fits_r,ra,dec,np.floor(size))
if ix == 0:
stack = np.zeros(fimg.shape, dtype=float)
mask = np.zeros(stack.shape, dtype=float)
stack=stack+fimg
mask=mask+_mask
plt.subplot(nx,ny,sp)
# use img_scale
scimg=img_scale.asinh(fimg,scale_min=np.min(fimg[_mask>0]))
im=plt.imshow(scimg, cmap=_cmap,origin='lower')
sp+=1
ix=np.where(mask > 0.1)
stack[ix]=stack[ix]/mask[ix]
stack[mask==0]=0
plt.subplot(nx,ny,sp)
#scimg=loggray(stack,mask)
scimg=img_scale.asinh(stack,scale_min=np.min(stack[mask>0]))
im=plt.imshow(scimg,cmap=_cmap,origin='lower')
plt.title('stack')
img=skimage.io.imread(url)
sp+=1
plt.subplot(nx,ny,sp)
plt.imshow(img)
plt.title("JPEG")
NBINS=100
histogram = plt.hist(scimg.flat, NBINS)
plt.yscale('log')
Artefact visible in final co-added image is probably feature of too-simplistic algorithm. Fix this (please!).