In [1]:
#all examples I have shown so far require you to download data to your machine
#However with most examples just calling information contained in online databases what is the point?
#to do this it is convient to download the module astroquery
#to do this simply type the following command line
#pip install astroquery
import matplotlib.pyplot as plt
import numpy as np
from astroquery.skyview import SkyView
#Grabed this from the internet for image rescaling however it can be left out if you have your own routines
import img_scale
In [2]:
#Now lets see what images we can query from SkyView
#if you want to know what commands get images can take run the following command
#help(SkyView.get_images)
#Lets look at what has been called my favorite cluster (really it is an association)
#(or the second large (>3000 stars) star forming region within 1kpc of the Sun)
img = SkyView.get_images(position='22:57:00,62:38:00',survey=['DSS2 Blue','DSS2 IR','DSS2 Red'],pixels='2400,2400',coordinates='J2000',grid=True,gridlabels=True)
#Another example would be
#img = SkyView.get_images(position='22:57:00,62:38:00',survey=['2MASS-J','2MASS-H','2MASS-K'])
In [3]:
#lets do a quick view to see some idea what we just downloaded
#lets look at all three data
fig,ax = plt.subplots(ncols=3,figsize=(24,8))
plot = ax[0].imshow(img[0][0].data,vmax=np.max(img[0][0].data)*.65,vmin=np.max(img[2][0].data)*.45)
plot1 = ax[1].imshow(img[1][0].data,vmax=np.max(img[1][0].data)*.75,vmin=np.max(img[1][0].data)*.4)
plot2 = ax[2].imshow(img[2][0].data,vmax=np.max(img[2][0].data)*.485,vmin=np.max(img[2][0].data)*.25)
In [4]:
#okay now lets make a 3-D image
#blue will be the DSS blue data
b = img[0][0]
#green will be the DSS Near IR data
g = img[1][0]
#red will be the DSS Red data
r = img[2][0]
#However just stacking everything is urealistic usually, but nice things have been done for us
#as we quickly can see by checking there shapes
print(b.data.shape,g.data.shape,r.data.shape)
In [5]:
#Quickly check what is in the header
b.header
Out[5]:
In [6]:
#first lets find out the coverage in space of the images we have just set up
def coverage(dat):
try:
print(dat.header['CRVAL1'],dat.header['CRVAL2'],dat.header['CRVAL1']+
dat.header['NAXIS1']*dat.header['CD1_1'],dat.header['CRVAL2']+dat.header['NAXIS2']*dat.header['CD2_2'])
#Because people love using different ways to do the same thing
except KeyError:
print(b.header['CRVAL1']-(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL1']+(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL2']+(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2'],
b.header['CRVAL2']-(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2'])
#Again these are already matched so work is avoided
coverage(b)
coverage(g)
coverage(r)
In [7]:
#Do some seting up of the data so we can stack it
bi = (b.data)
ri = (r.data)
gi = (g.data)
#Rescale the data
ri = img_scale.linear(ri,scale_min=.22*np.max(ri),scale_max=.45*np.max(ri))
gi = img_scale.linear(gi,scale_min=.44*np.max(gi),scale_max=.65*np.max(gi))
bi = img_scale.linear(bi,scale_min=.40*np.max(bi),scale_max=.57*np.max(bi))
#Replace any value higher than 1 with 1
replace = bi > 1
bi[replace] = 1
replace = gi > 1
gi[replace] = 1
replace = ri > 1
ri[replace] = 1
#Stack the Images and set up as 8-bit integers (highest value 2^8 = 256)
imgs = (np.dstack((ri,gi,bi))*255.999).astype(np.uint8)
fig, ax = plt.subplots(figsize=(10,10))
ax.imshow(imgs,interpolation='none',aspect='equal',extent=[b.header['CRVAL1']-(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL1']+(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL2']+(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2'],
b.header['CRVAL2']-(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2']])
fig.gca().invert_yaxis()
#Ultimately we would like to display the plot in hours, minutes, and seconds rather than degrees because we are people
ax.set_xticks([344.0,344.25,344.5])
ax.set_xticklabels(['22$^h$:56$^m$','22$^h$:57$^m$','22$^h$:58$^m$'])
ax.set_yticks([62.3833333,62.633333333,62.883333333333])
ax.set_yticklabels(['62:23$^\prime$','62:38$^\prime$','62:53$^\prime$'])
Out[7]:
In [8]:
#Now Let us zoom in on the interesting region and we can see stars blowing away their natal gas
fig, ax = plt.subplots(figsize=(10,10))
ax.imshow(imgs,interpolation='none',aspect='equal',extent=[b.header['CRVAL1']-(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL1']+(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL2']+(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2'],
b.header['CRVAL2']-(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2']])
ax.set_xlim([344.40,344.])
ax.set_ylim([62.5,62.883333333333])
#Ultimately we would like to display the plot in hours, minutes, and seconds rather than degrees because we are people
ax.set_xticks([344.0,344.25])
ax.set_xticklabels(['22$^h$:56$^m$','22$^h$:57$^m$'])
ax.set_yticks([62.5,62.883333333333])
ax.set_yticklabels(['62:30$^\prime$','62:53$^\prime$'])
Out[8]:
In [9]:
#Now all other plotting commands still work in the images so we can overplot some data on the images
from astroquery.simbad import Simbad
result_table = Simbad.query_region("22 57 00 +62 38 00", radius='0d3m0s')
print(result_table)
In [10]:
print(result_table.colnames)
In [11]:
#Or we can query an object we know is in the region
Bstar = Simbad.query_object("HD 217061")
print(Bstar)
print(Bstar.colnames)
print(Bstar['RA'],Bstar['DEC'])
In [12]:
#And another just for fun
Ostar = Simbad.query_object("HD 217086")
print(Ostar['RA'],Ostar['DEC'])
In [13]:
#Now let us put it on the graph
fig, ax = plt.subplots(figsize=(10,10))
ax.imshow(imgs,interpolation='none',aspect='equal',extent=[b.header['CRVAL1']-(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL1']+(b.header['NAXIS1']-b.header['CRPIX1'])*b.header['CDELT1'],
b.header['CRVAL2']+(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2'],
b.header['CRVAL2']-(b.header['NAXIS2']-b.header['CRPIX2'])*b.header['CDELT2']])
fig.gca().invert_yaxis()
#Bstar data convert RA and Dec to degrees
ras = Bstar['RA'].data[0]
decs = Bstar['DEC'].data[0]
rab = (np.float64(ras[:3])+np.float64(ras[3:5])/60.+np.float64(ras[5:])/3600.)*15.
decb = np.float64(decs[:3])+np.float64(decs[4:6])/60.+np.float64(decs[6:])/3600.
#Ostar data convert RA and Dec to degrees
ras = Ostar['RA'].data[0]
decs = Ostar['DEC'].data[0]
rao = (np.float64(ras[:3])+np.float64(ras[3:5])/60.+np.float64(ras[5:])/3600.)*15.
deco = np.float64(decs[:3])+np.float64(decs[4:6])/60.+np.float64(decs[6:])/3600.
#The image is projected Gnomomically so we need to correct for that to get the circles in the right place
def fix_project(ra,dec):
a = np.cos(np.radians(dec))*np.cos(np.radians(ra-b.header['CRVAL1']))
f = 1./b.header['CDELT2']*(180./np.pi)/(np.sin(np.radians(b.header['CRVAL2']))*np.sin(np.radians(dec))+a*np.cos(np.radians(b.header['CRVAL2'])))
decn = -f*(np.cos(np.radians(b.header['CRVAL2']))*np.sin(np.radians(dec))-a*np.sin(np.radians(b.header['CRVAL2'])))
ran = -f*np.cos(np.radians(dec))*np.sin(np.radians(ra-b.header['CRVAL1']))
ran = b.header['CRVAL1']+(ran)*b.header['CDELT1']
decn = b.header['CRVAL2']-(decn)*b.header['CDELT2']
return ran,decn
#Correct pixel placement of circle center for image projection
rab,decb = fix_project(rab,decb)
rao,deco = fix_project(rao,deco)
#Create circle
circleb = plt.Circle((rab,decb),.015,color='green',fill=False,linewidth=3)
circleo = plt.Circle((rao,deco),.015,color='white',fill=False,linewidth=3)
#draw circles on the plot
fig.gca().add_artist(circleb)
fig.gca().add_artist(circleo)
#Ultimately we would like to display the plot in hours, minutes, and seconds rather than degrees because we are people
ax.set_xticks([344.0,344.25,344.5])
ax.set_xticklabels(['22$^h$:56$^m$','22$^h$:57$^m$','22$^h$:58$^m$'])
ax.set_yticks([62.3833333,62.633333333,62.883333333333])
ax.set_yticklabels(['62:23$^\prime$','62:38$^\prime$','62:53$^\prime$'])
Out[13]:
