Image Scaling module

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'])
WARNING: AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead. [astropy.config.configuration]
WARNING:astropy:AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead.
/wandajune_home/jakub/PYTHON/loopy/lib/python2.7/site-packages/astroquery/utils/commons.py:373: UserWarning: FITS files must be read as binaries; error is likely.
  warnings.warn("FITS files must be read as binaries; error is likely.")

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)
((2400, 2400), (2400, 2400), (2400, 2400))

In [5]:
#Quickly check what is in the header
b.header
Out[5]:
SIMPLE  =                    T / Written by SkyView Tue Apr 28 16:07:46 EDT 2015
BITPIX  =                  -32 / 4 byte floating point                          
NAXIS   =                    2 / Two dimensional image                          
NAXIS1  =                 2400 / Width of image                                 
NAXIS2  =                 2400 / Height of image                                
CRVAL1  =               344.25 / Reference longitude                            
CRVAL2  =    62.63333333333333 / Reference latitude                             
RADESYS = 'FK5     '           / Coordinate system                              
EQUINOX =               2000.0 / Epoch of the equinox                           
CTYPE1  = 'RA---TAN'           / Coordinates -- projection                      
CTYPE2  = 'DEC--TAN'           / Coordinates -- projection                      
CRPIX1  =               1200.5 / X reference pixel                              
CRPIX2  =               1200.5 / Y reference pixel                              
CDELT1  =         -2.777777E-4 / X scale                                        
CDELT2  =          2.777777E-4 / Y scale                                        
COMMENT                                                                         
COMMENT  SkyView Survey metadata                                                
COMMENT                                                                         
COMMENT Provenance:  Data taken by ROE, AAO, and CalTech, Compression           
COMMENT                       and distribution by Space Telescope Science Inst  
COMMENT           itute.                                                        
COMMENT Copyright:   There are multiple copyright holders depending         up  
COMMENT           on the source plate or plates.          See full copyright    
COMMENT                  notices The coverage link below describes the          
COMMENT            coverage for each element of the surveys when the original   
COMMENT           plate used is         not included in the SkyView files..     
COMMENT Regime:      Optical                                                    
COMMENT NSurvey:     1                                                          
COMMENT Frequency:   637 THz                                                    
COMMENT Bandpass:    517-967 THz                                                
COMMENT Coverage:    All-sky, but some data not yet be processed.               
COMMENT PixelScale:  1" in the north 1.7" in the south                          
COMMENT PixelUnits:  Pixel values are given as scaled densities                 
COMMENT Resolution:  Depends on plate. Typically better than 2".                
COMMENT Coordinates: Equatorial                                                 
COMMENT Projection:  Schmidt                                                    
COMMENT Equinox:     2000                                                       
COMMENT Epoch:       1950-1998                                                  
COMMENT Reference:   Some information on the DSS2 is given in            McLea  
COMMENT           n, 2000, The Second Generation Guide Star Catalog.            
COMMENT             Characteristics of the DSS surveys are summarized in        
COMMENT                this document.                                           
COMMENT                                                                         
COMMENT Survey specific cards                                                   
COMMENT                                                                         
SURVEY  = 'DSS2BLUE'                                                            
HISTORY                                                                         
HISTORY  Settings used in processing:                                           
HISTORY                                                                         
HISTORY cache = /skyview8/cache/surveys/                                        
HISTORY catalogids = on                                                         
HISTORY coordinates = J2000.0                                                   
HISTORY deedger = skyview.process.Deedger                                       
HISTORY descriptionxslt = cgifiles/description.xsl                              
HISTORY dummy                                                                   
HISTORY equinox = 2000                                                          
HISTORY fileprefix = http://skyview.gsfc.nasa.gov/surveys/dss                   
HISTORY finalpostprocessor = skyview.ij.IJProcessor,skyview.data.BoxSmoother,sky
HISTORY float = on                                                              
HISTORY footertemplate = cgifiles/skyfooter.html                                
HISTORY gallerydir = /userimages                                                
HISTORY galleryxslt = cgifiles/gallerymultipage.xsl                             
HISTORY grid = True                                                             
HISTORY gridlabels                                                              
HISTORY headertemplate = cgifiles/skyheader.html                                
HISTORY htmlwriter = skyview.request.HTMLWriter                                 
HISTORY imagefactory = skyview.survey.DSSImageFactory                           
HISTORY imagesize = 6                                                           
HISTORY localurl = http://skyview.gsfc.nasa.gov/surveys,/skyview/htdocs-VME/surv
HISTORY lut = colortables/b-w-linear.bin                                        
HISTORY lutcbarpath = ../images/colorbars/                                      
HISTORY mosaicker = skyview.process.Mosaicker                                   
HISTORY name = 2nd Digitized Sky Survey (Blue)                                  
HISTORY noexit =                                                                
HISTORY nullimagedir = ../images/nodata                                         
HISTORY output = ../../tempspace/fits/skv4501859425211_1                        
HISTORY outputroot = ../../tempspace/fits                                       
HISTORY pixels = 2400,2400                                                      
HISTORY position = 22:57:00,62:38:00                                            
HISTORY postprocessor = skyview.ij.IJProcessor,skyview.request.HTMLWriter       
HISTORY projection = Tan                                                        
HISTORY quicklook = jpeg                                                        
HISTORY requested_coords = 344.25, 62.63333333333333                            
HISTORY reqxpos = 344.25                                                        
HISTORY reqypos = 62.63333333333333                                             
HISTORY resolver = SIMBAD-NED                                                   
HISTORY rgbtemplate = cgifiles/skyrgb.html                                      
HISTORY rgbwriter = skyview.request.RGBWriter                                   
HISTORY sampler = Default                                                       
HISTORY savebysurvey                                                            
HISTORY scale = 0.0002777777                                                    
HISTORY scaling = Log Sqrt Linear HistEq LogLog                                 
HISTORY settingsupdaters = BatchCompatibility,SettingsFixer,skyview.request.Toas
HISTORY shortname = DSS2B, DSS2 Blue                                            
HISTORY sia_header = ./sia.header                                               
HISTORY siabase = http://skyview.gsfc.nasa.gov/cgi-bin/images?                  
HISTORY siaimagetimeout = 300000                                                
HISTORY size = 0.6666664800000001,0.6666664800000001                            
HISTORY survey = DSS2 Blue,DSS2 IR,DSS2 Red                                     
HISTORY surveyfinder = skyview.survey.XMLSurveyFinder                           
HISTORY surveymanifest = surveys/survey.manifest                                
HISTORY surveyregimes = Radio,Millimeter,Infrared,Optical,Ultraviolet,X-ray,Gamm
HISTORY surveysheader = cgifiles/survey.header                                  
HISTORY surveytemplate = cgifiles/skysurvey.html                                
HISTORY url.heasarcbase = http://heasarc.gsfc.nasa.gov/xamin/vo/cone?showoffsets
HISTORY url.ned = http://nedwww.ipac.caltech.edu/cgi-bin/nph-NEDobjsearch?search
HISTORY url.simbad = http://simbad.u-strasbg.fr/simbad-conesearch.pl?           
HISTORY url.vizierbase = http://vizier.u-strasbg.fr/viz-bin/votable/-dtd/-A?-out
HISTORY urlcoordinates = http://heasarc.gsfc.nasa.gov/cgi-bin/Tools/convcoord/co
HISTORY urllocalhelp = http://skyview.gsfc.nasa.gov/help/help.html              
HISTORY usedsswcs                                                               
HISTORY version = 3.1.12                                                        
HISTORY webrootpath = /skyview/htdocs-VME                                       
HISTORY                                                                         
HISTORY  Map generated at: Tue Apr 28 16:07:46 EDT 2015                         
HISTORY                                                                         
HISTORY  Resampler used: NNSampler                                              
HISTORY                                                                         
HISTORY                                                                         
HISTORY Image mosaicking using skyview.geometry.Mosaicker                       
HISTORY                                                                         
HISTORY   Used image:xj147                                                      
HISTORY   Used image:xj109                                                      
HISTORY                                                                         
HISTORY                                                                         
HISTORY Edge adjustments applied (skyview.geometry.DeedgerList                  
HISTORY                                                                         
HISTORY      Image xj109 offset by -1833.0                                      
HISTORY                                                                         
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)
(344.58319435115, 343.91680564885, 62.96652768448333, 62.300138982183334)
(344.58319435115, 343.91680564885, 62.96652768448333, 62.300138982183334)
(344.58319435115, 343.91680564885, 62.96652768448333, 62.300138982183334)

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$']) 
img_scale : linear
img_scale : linear
img_scale : linear

Out[7]:
[<matplotlib.text.Text at 0x7ff8025eeed0>,
 <matplotlib.text.Text at 0x7ff806a8fd10>,
 <matplotlib.text.Text at 0x7ff806ead6d0>]
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]:
[<matplotlib.text.Text at 0x7ff806b23b90>,
 <matplotlib.text.Text at 0x7ff806bcda90>]
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)
        MAIN_ID              RA      ... COO_WAVELENGTH     COO_BIBCODE    
                          "h:m:s"    ...                                   
----------------------- ------------ ... -------------- -------------------
2MASS J22570070+6238018  22 57 00.70 ...              I 2003yCat.2246....0C
2MASS J22570049+6238116 22 57 00.491 ...              N 2003yCat.2246....0C
2MASS J22570191+6237514  22 57 01.91 ...              I 2003yCat.2246....0C
2MASS J22570133+6237422  22 57 01.34 ...              I 2003yCat.2246....0C
              [TOH95] C   22 57 01.3 ...                1995A&A...303..881T
2MASS J22565915+6238192 22 56 59.151 ...              N 2003yCat.2246....0C
2MASS J22565889+6237402  22 56 58.90 ...              I 2003yCat.2246....0C
2MASS J22570003+6238240  22 57 00.04 ...              I 2003yCat.2246....0C
2MASS J22565633+6238037 22 56 56.337 ...              N 2003yCat.2246....0C
2MASS J22570220+6238251  22 57 02.20 ...              I 2003yCat.2246....0C
                    ...          ... ...            ...                 ...
2MASS J22570466+6240530  22 57 04.66 ...              I 2003yCat.2246....0C
       [MCD93] Cep B 22   22 57 11.6 ...                1993A&A...273..619M
2MASS J22565866+6240559  22 56 58.66 ...              I 2003yCat.2246....0C
2MASS J22572473+6238501  22 57 24.74 ...              I 2003yCat.2246....0C
2MASS J22564566+6235349  22 56 45.66 ...              I 2003yCat.2246....0C
2MASS J22563671+6236498  22 56 36.72 ...              I 2003yCat.2246....0C
2MASS J22570223+6240570 22 57 02.235 ...              I 2003yCat.2246....0C
2MASS J22564369+6240171 22 56 43.690 ...              I 2003yCat.2246....0C
2MASS J22563464+6237303 22 56 34.645 ...              I 2003yCat.2246....0C
               [C51] 99      22 57.0 ...              O 1976A&AS...25...25D
Length = 235 rows

In [10]:
print(result_table.colnames)
['MAIN_ID', 'RA', 'DEC', 'RA_PREC', 'DEC_PREC', 'COO_ERR_MAJA', 'COO_ERR_MINA', 'COO_ERR_ANGLE', 'COO_QUAL', 'COO_WAVELENGTH', 'COO_BIBCODE']

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'])
 MAIN_ID        RA           DEC      ... COO_WAVELENGTH     COO_BIBCODE    
             "h:m:s"       "d:m:s"    ...                                   
--------- ------------- ------------- ... -------------- -------------------
HD 217061 22 56 42.5531 +62 37 29.555 ...              O 2007A&A...474..653V
['MAIN_ID', 'RA', 'DEC', 'RA_PREC', 'DEC_PREC', 'COO_ERR_MAJA', 'COO_ERR_MINA', 'COO_ERR_ANGLE', 'COO_QUAL', 'COO_WAVELENGTH', 'COO_BIBCODE']
(<MaskedColumn name='RA' dtype='unicode416' unit='"h:m:s"' description=u'Right ascension' length=1>
22 56 42.5531, <MaskedColumn name='DEC' dtype='unicode416' unit='"d:m:s"' description=u'Declination' length=1>
+62 37 29.555)

In [12]:
#And another just for fun
Ostar = Simbad.query_object("HD 217086")
print(Ostar['RA'],Ostar['DEC'])
(<MaskedColumn name='RA' dtype='unicode416' unit='"h:m:s"' description=u'Right ascension' length=1>
22 56 47.1903, <MaskedColumn name='DEC' dtype='unicode416' unit='"d:m:s"' description=u'Declination' length=1>
+62 43 37.648)

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]:
[<matplotlib.text.Text at 0x7ff8070cac10>,
 <matplotlib.text.Text at 0x7ff8070dfcd0>,
 <matplotlib.text.Text at 0x7ff806b70090>]