#
#

## Keyword Extraction

#
# # Extract and print out the `TELESCOP` value from the header. # # Next, extract the `WAVELNTH` and `WAVEUNIT` values, use these to construct an astropy Quantity object for the wavelength of this image. # # # In[ ]: header['TELESCOP'] # In[ ]: import astropy.units as u u.Quantity(header['WAVELNTH'], unit=header['WAVEUNIT']) # We could now sit down and work out how to convert from a pixel coordinate to a physical coordinate described by this header (Helioprojective). # # However, we can cheat and just use Astropy. # In[ ]: from astropy.wcs import WCS wcs = WCS(header) # We can convert from pixel to world coordinate: # In[ ]: wcs.wcs_pix2world(((100, 100),), 0) # Or back again: # In[ ]: wcs.wcs_world2pix([[ 3.59725669e+02, -2.74328093e-01]], 0) # The last parameter to the two above examples is the 'origin' parameter. It is a flag that tells WCS if you indexes should be 0-based (like numpy) or 1-based (like FITS). # Here we are using 0 as we want to convert to and from numpy indexes of the array. #
#
#

## How large is the image?

#
#
# To get a little practise using Astropy's WCS calculate the world coordinates of the following pixels: # ``` # [-500, 0] # [500, 500] # [0, 0] # ``` #
#
# In[ ]: print(wcs.wcs_pix2world(((-500, 0),), 1)) print(wcs.wcs_pix2world(((500, 500),), 1)) print(wcs.wcs_pix2world(((0, 0),), 1)) # ## Plotting with wcsaxes # In this section we are going to use the wcsaxes package to make WCS aware image plots. # In[ ]: import wcsaxes # For this example we are going to use a Hubble image. # In[ ]: from astropy.io import fits hdulist = fits.open('./h_n4571_f555_mosaic.fits.gz') hdulist # In[ ]: wcs = WCS(hdulist[0].header) # In[ ]: get_ipython().run_line_magic('matplotlib', 'inline') import matplotlib.pyplot as plt # In[ ]: ax = plt.subplot(111, projection=wcs) ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower') # This image now has physcial labels in the native coordinate system of the image. We can see what the coordinate system and projection of this image is using the 'CTYPE' header entries we saw earlier. # In[ ]: print(hdulist[0].header['CTYPE1'], hdulist[0].header['CTYPE2']) # We can tell that this is in the FK5 coordinate system by the presence of a 'equinox' entry in the header: # In[ ]: hdulist[0].header['equinox'] # There is also a quick way to generate an Astropy coordinate frame from a WCS object, which confirms this diagnosis. # In[ ]: from astropy.wcs.utils import wcs_to_celestial_frame wcs_to_celestial_frame(wcs) # for more information on the very useful `astropy.coordinates` module see http://docs.astropy.org/en/stable/coordinates/ #
#
#

#
#
# Now we have a nice plot with physically meaningful ticks, we should label our axes. #
# # Add labels to the axes saying "Right Ascension [degrees]" and "Declination [degrees]" #
# # Also overlay a coordinate grid using: # `ax.coords.grid()` # Look up the documentation for this method to see what parameters you can specify. #
#
# In[ ]: ax = plt.subplot(111, projection=wcs) ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower') ax.set_xlabel("Right Ascension [degrees]") ax.set_ylabel("Declination [degrees]") ax.coords.grid(color='white', alpha=0.5, linestyle='solid') # Now we have a nice plot, we can do a couple of things to plot. # ### Overplotting in Pixel Coordinates # In[ ]: ax = plt.subplot(111, projection=wcs) ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower') ax.set_xlabel("Right Ascension [degrees]") ax.set_ylabel("Declination [degrees]") ax.coords.grid(color='white', alpha=0.5, linestyle='solid') ax.plot(3000, 3000, 'o') # ### Overplotting in World Coordinates # In[ ]: ax = plt.subplot(111, projection=wcs) ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower') ax.set_xlabel("Right Ascension [degrees]") ax.set_ylabel("Declination [degrees]") ax.coords.grid(color='white', alpha=0.5, linestyle='solid') ax.set_autoscale_on(False) ax.plot(3000, 3000, 'o') # Overplot in FK5 in Degrees ax.plot(189.25, 14.23, 'o', transform=ax.get_transform('world')) #
#
#