Investigation of SDSS image of NGC 450

In this notebook I will import the SDSS image of NGC 450 in the g band and investigate some properties of the object.

TODO: 3D visualisation of flux counts. Compare to Gaussian and Sersic profiles for Galfit value of n. Visualise Galfit fits for one vs two Sersic profiles.


In [41]:
from regphot import git_version
print("This notebook was run with regphot version: \n{}".format(git_version()))


This notebook was run with regphot version: 
62b9b19 (Thu Apr 6 19:19:04 2017 +0100) [with local modifications]

In [42]:
#Import relevant modules
import matplotlib.pyplot as plt
%matplotlib inline
from astropy.io import fits
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import numpy as np
import aplpy
from regphot import galfit

Import the whole SDSS frame


In [43]:
full = fits.open('frame-g-004858-1-0480.fits')
wcs = WCS('frame-g-004858-1-0480.fits')
apfig = aplpy.FITSFigure(full)
apfig.show_grayscale()
apfig.add_scalebar(0.05, "0.05 degrees", color='white', corner='top right')
apfig.set_theme('publication')
apfig.tick_labels.set_xformat('hh:mm:ss.ssss') #ddd.dddddd')
apfig.tick_labels.set_yformat('dd:mm:ss.ss') #ddd.dddddd')
apfig.ticks.show()
#apfig.ticks.set_xspacing(0.05)
apfig.tick_labels.show()
#apfig.add_grid()


WARNING: FITSFixedWarning: RADECSYS= 'ICRS ' / International Celestial Ref. System 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
INFO: Auto-setting vmin to -9.657e-02 [aplpy.core]
INFO: Auto-setting vmax to  3.577e-01 [aplpy.core]
/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/aplpy/wcs_util.py:515: UserWarning: Pixels are not square, using an average pixel scale
  warnings.warn("Pixels are not square, using an average pixel scale")

Plot 1D profile across centre of main galaxy


In [54]:
print(len(full[0].data))
y=full[0].data[1210]
x=np.arange(0,len(y),1)
plt.plot(x, y)
plt.xlim([1500,len(y)])
plt.savefig('1dsimple.pdf')


1489

In [45]:
full[0].header


Out[45]:
SIMPLE  =                    T /                                                
BITPIX  =                  -32 / 32 bit floating point                          
NAXIS   =                    2                                                  
NAXIS1  =                 2048                                                  
NAXIS2  =                 1489                                                  
EXTEND  =                    T /Extensions may be present                       
BZERO   =              0.00000 /Set by MRD_SCALE                                
BSCALE  =              1.00000 /Set by MRD_SCALE                                
TAI     =        4602733435.14 / 1st row - Number of seconds since Nov 17 1858  
RA      =            19.456261 / 1st row - Right ascension of telescope boresigh
DEC     =            0.000000  / 1st row - Declination of telescope boresight (d
SPA     =              90.000  / 1st row - Camera col position angle wrt north (
IPA     =              69.299  / 1st row - Instrument rotator position angle (de
IPARATE =              0.0000  / 1st row - Instrument rotator angular velocity (
AZ      =            335.10436 / 1st row - Azimuth  (encoder) of tele (0=N?) (de
ALT     =            54.749787 / 1st row - Altitude (encoder) of tele        (de
FOCUS   =           -237.60000 / 1st row - Focus piston (microns?)              
DATE-OBS= '2004-09-24'         / 1st row - TAI date                             
TAIHMS  = '09:03:55.14'        / 1st row - TAI time (HH:MM:SS.SS) (TAI-UT = appr
COMMENT  TAI,RA,DEC,SPA,IPA,IPARATE,AZ,ALT,FOCUS at reading of col 0, row 0     
ORIGIN  = 'SDSS    '                                                            
TELESCOP= '2.5m    '                                                            
TIMESYS = 'TAI     '                                                            
RUN     =                 4858 / Run number                                     
FRAME   =                  488 / Frame sequence number within the run           
CCDLOC  =                   51 / Survey location of CCD (e.g., rowCol)          
STRIPE  =                   82 / Stripe index number (23 <--> eta=0)            
STRIP   = 'N       '           / Strip in the stripe being tracked.             
FLAVOR  = 'science '           / Flavor of this run                             
OBSERVER= 'long    '           / Observer                                       
SYS_SCN = 'mean    '           / System of the scan great circle (e.g., mean)   
EQNX_SCN=           2000.00    / Equinox of the scan great circle. (years)      
NODE    =           95.00000   / RA of the great circle's ascending node (deg)  
INCL    =           0.00000    / Great circle's inclination wrt cel. eq. (deg)  
XBORE   =           -22.74     / Boresight x offset from the array center (mm)  
YBORE   =           0.00       / Boresight x offset from the array center (mm)  
OBJECT  = '82 N    '           / e.g., 'stripe 50.6 degrees, north strip'       
EXPTIME = '53.907456'          / Exposure time (seconds)                        
SYSTEM  = 'FK5     '           / System of the TCC coordinates (e.g., mean)     
CCDMODE = 'DRIFT   '           / 'STARING' or 'DRIFT'                           
C_OBS   =                26322 / CCD row clock rate (usec/row)                  
COLBIN  =                    1 / Binning factor perpendicular to the columns    
ROWBIN  =                    1 / Binning factor perpendicular to the rows       
DAVERS  = 'v14_47  '           / Version of DA software                         
SCDMETHD= 'sqrtDynamic'        / scdMethod                                      
SCDWIDTH=                 1280 / scdDisplayWidth                                
SCDDECMF=                    1 / scdDecimateFactor                              
SCDOFSET=                  410 / scdDisplayOffset                               
SCDDYNTH=                  -20 / scdDynamicThresh                               
SCDSTTHL=                   30 / scdStaticThreshL                               
SCDSTTHR=                   30 / scdStaticThreshR                               
SCDREDSZ=                  515 / scdReduceSize                                  
SCDSKYL =                 2128 / scdSkyLeft                                     
SCDSKYR =                 2129 / scdSkyRight                                    
COMMENT  CCD-specific parameters                                                
CAMROW  =                    5 / Row in the imaging camera                      
BADLINES=                    0 / Number of bad lines in frame                   
EQUINOX =              2000.00 /                                                
SOFTBIAS=                 1000 / software "bias" added to all DN                
BUNIT   = 'nanomaggy'          / 1 nanomaggy = 3.631e-6 Jy                      
FILTER  = 'g       '           / filter used                                    
CAMCOL  =                    1 / column in the imaging camera                   
VERSION = 'v5_6_3  '                                                            
DERV_VER= 'NOCVS:v8_23'                                                         
ASTR_VER= 'NOCVS:v5_24'                                                         
ASTRO_ID= '2009-05-17T03:25:23 12613'                                           
BIAS_ID = 'PS      '                                                            
FRAME_ID= '2009-09-21T08:59:15 03168'                                           
KO_VER  = 'devel   '                                                            
PS_ID   = '2009-05-16T18:09:56 03674 camCol 1'                                  
ATVSN   = 'NOCVS:v5_24'        / ASTROTOOLS version tag                         
RADECSYS= 'ICRS    '           / International Celestial Ref. System            
CTYPE1  = 'RA---TAN'           /Coordinate type                                 
CTYPE2  = 'DEC--TAN'           /Coordinate type                                 
CUNIT1  = 'deg     '           /Units                                           
CUNIT2  = 'deg     '           /Units                                           
CRPIX1  =        1025.00000000 /X of reference pixel                            
CRPIX2  =        745.000000000 /Y of reference pixel                            
CRVAL1  =        18.8251174090 /RA of reference pixel (deg)                     
CRVAL2  =      -0.944890439256 /Dec of reference pixel (deg)                    
CD1_1   =   -9.86650666540E-09 /RA deg per column pixel                         
CD1_2   =    0.000110015128778 /RA deg per row pixel                            
CD2_1   =    0.000109968022201 /Dec deg per column pixel                        
CD2_2   =    2.00246429838E-09 /Dec deg per row pixel                           
HISTORY GSSSPUTAST: Sep 21 08:59:21 2009                                        
COMMENT  Calibration parameters                                                 
COMMENT  Floats truncated at 10 binary digits with FLOATCOMPRESS                
NMGY    =           0.00332422 / Calibration factor [nMgy per count]            
NMGYIVAR=            0.0508784 / Calibration factor inverse variance            
VERSIDL = '7.0     '           / Version of IDL                                 
VERSUTIL= 'v5_5_5  '           / Version of idlutils                            
VERSPOP = 'v1_11_1 '           / Version of photoop product                     
PCALIB  = '/clusterfs/riemann/raid006/bosswork/groups/boss/calib/2009-06-14/cal'
PSKY    = '/clusterfs/riemann/raid006/bosswork/groups/boss/photo/sky' / Value of
RERUN   = '301     '           / rerun                                          
HISTORY SDSS_FRAME_ASTROM: Astrometry fixed for dr9 Thu Jun 21 03:27:56 2012    

Extract 1D profile across both galaxies


In [55]:
# from http://stackoverflow.com/questions/7878398/how-to-extract-an-arbitrary-line-of-values-from-a-numpy-array
#-- Extract the line...
# Make a line with "num" points...
x0, y0 = 1650, 1000  # These are in _pixel_ coordinates!!
x1, y1 = 1940, 1440 
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length), np.linspace(y0, y1, length)

# Extract the values along the line (why seemingly wrong way round?)
zi = full[0].data[y.astype(np.int), x.astype(np.int)]

#-- Plot...
fig = plt.plot(zi)

#plt.show()
plt.savefig('1dblended.pdf')


Show line used to generate 1D flux graph


In [56]:
apfig = aplpy.FITSFigure(full)
apfig.show_grayscale()
apfig.add_scalebar(0.05, "0.05 degrees", color='white', corner='top right')
apfig.set_theme('publication')
apfig.tick_labels.set_xformat('hh:mm:ss.ssss') #ddd.dddddd')
apfig.tick_labels.set_yformat('dd:mm:ss.ss') #ddd.dddddd')
apfig.ticks.show()
#apfig.ticks.set_xspacing(0.05)
apfig.tick_labels.show()
#apfig.add_grid()
ra0, dec0 = wcs.all_pix2world(x0, y0, 0) #, 1500 # These are in _pixel_ coordinates!!
ra1, dec1 = wcs.all_pix2world(x1,y1, 0 )#, 2045
iline = np.array([[ra0, ra1],[dec0,dec1]])
apfig.show_lines([iline], color = 'r')
apfig.savefig('lineonmap.pdf')


WARNING: FITSFixedWarning: RADECSYS= 'ICRS ' / International Celestial Ref. System 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
INFO: Auto-setting vmin to -9.657e-02 [aplpy.core]
INFO: Auto-setting vmax to  3.577e-01 [aplpy.core]
/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/aplpy/wcs_util.py:515: UserWarning: Pixels are not square, using an average pixel scale
  warnings.warn("Pixels are not square, using an average pixel scale")
INFO: Auto-setting resolution to 264.258 dpi [aplpy.core]

Find all 'peak' pixels

For every pixel in the map mark it true if every surrounding pixel is lower in value and false otherwise. This will show us the peak pixels. We can then investigate the positions of these pixels in colour space.


In [ ]: