{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Images and Coordinates in Astronomy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lesson we are going to look at aspects of processing and viewing images specific to Astronomy and Solar Astronomy. By the end of this lesson you should understand:\n", "\n", "* Projected Coordinate Systems in Images\n", "* World Coordinate Systems\n", "* Using WCS to calculate coordinates in images\n", "* Plotting images with WCS in images\n", "* Using SunPy Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Projected Coordinate Systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When taking images of the sky, we are projecting the spherical celestial coordinate system onto a 2-dimensional plane, which means that there is no simple linear relation between pixel coordinates and celestial coordinates\n", "\n", "There are multiple coordinate systems used to describe the locations in 2D and 3D space for both Astronomy and Solar Physics. We shall use a couple of these systems here as examples but if you want to know more about them there are many of resources avalible.\n", "\n", "![Spherical Coordinates](coord_inset.png)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Physical Coordinate Systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Projected coordinate systems are one type of physical coordinate systems, they are the one we will be focusing on in this lesson due to their applicability to imaging data. Astropy and SunPy support representing point in many different physical coordinate systems, both projected and fully 3D, such as ICRS or Helioprojective." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import astropy.units as u\n", "import sunpy.coordinates\n", "from astropy.coordinates import SkyCoord" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "icrs = SkyCoord(10*u.deg, 86*u.deg, frame='icrs')\n", "icrs" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hpc = SkyCoord(100*u.arcsec, 700*u.arcsec, obstime=\"now\", frame='helioprojective')\n", "hpc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hpc.transform_to('heliographic_stonyhurst')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hpc.transform_to(\"icrs\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### World Coordinate System\n", "\n", "#### From pixels to physical coordinates" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The FITS files have a standard for describing the physical coordinate system associated with imaging data, this is called the world coordinate system or WCS, sometimes the specific FITS version of this is referred to as FITS-WCS.\n", "\n", "There are multiple papers describing the FITS-WCS standard for various types of data, there is a list here: http://fits.gsfc.nasa.gov/fits_wcs.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you learned in the previous lesson we can load FITS files with Astropy. To demonstrate a simple example of a FITS file with FITS-WCS information in the header we shall use an image from SunPy:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sunpy.data.sample import AIA_171_IMAGE\n", "from astropy.io import fits" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hdulist = fits.open(AIA_171_IMAGE)\n", "hdulist.verify('silentfix')\n", "hdulist[0].header" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see there are lots of keys in this and most other real world FITS headers. The ones we need to understand for FITS-WCS are:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reference Pixel and Coordinate:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "header = hdulist[0].header\n", "\n", "print(header['CRVAL1'], header['CRVAL2'])\n", "print(header['CRPIX1'], header['CRPIX2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pixel resolution (at the reference pixel):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(header['CDELT1'], header['CDELT2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rotation angle, in degress (at the reference pixel):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(header['CROTA2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Coordinate System and Projection:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(header['CTYPE1'], header['CTYPE2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Keyword Extraction

\n", "
\n", "\n", "\n", "
\n", "\n", "

Extract and print out the TELESCOP value from the header.

\n", "

Next, extract the WAVELNTH and WAVEUNIT values, use these to construct an astropy Quantity object for the wavelength of this image.

\n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Solution

\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "header['TELESCOP']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import astropy.units as u\n", "u.Quantity(header['WAVELNTH'], unit=header['WAVEUNIT'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could now sit down and work out how to convert from a pixel coordinate to a physical coordinate described by this header (Helioprojective).\n", "\n", "However, we can cheat and just use Astropy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.wcs import WCS" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs = WCS(header)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can convert from pixel to world coordinate:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs.wcs_pix2world(((100, 100),), 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or back again:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs.wcs_world2pix([[ 3.59725669e+02, -2.74328093e-01]], 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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).\n", "Here we are using 0 as we want to convert to and from numpy indexes of the array." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Where are my Pixels?\n", "\n", "```\n", "[-500, 0]\n", "[500, 500]\n", "[0, 0]\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Solution

\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(wcs.wcs_pix2world(((-500, 0),), 1))\n", "print(wcs.wcs_pix2world(((500, 500),), 1))\n", "print(wcs.wcs_pix2world(((0, 0),), 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### WCS and SkyCoord" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A common usecase for WCS + Coordinates is to store or transform from pixel coordinates to one or more different physical coordinates. Combining Astropy WCS and Coordinates makes this easy.\n", "\n", "Assuming we have the WCS object we created from the FITS header above we can get an astropy Coordinate Frame:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.wcs.utils import wcs_to_celestial_frame" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hpc_frame = wcs_to_celestial_frame(wcs)\n", "hpc_frame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then use this when creating `SkyCoord` objects:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SkyCoord(100*u.arcsec, -500*u.arcsec, frame=hpc_frame)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Note: For solar data you should always use `sunpy.map` for this*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coordinate Aware Plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this section we are going to use the `astropy.visualization.wcsaxes` subpackage to make WCS aware image plots." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example we are going to use a Hubble image." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.io import fits" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hdulist = fits.open('h_n4571_f555_mosaic.fits.gz')\n", "hdulist" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs = WCS(hdulist[0].header)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# HIDDEN\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax = plt.subplot(111, projection=wcs)\n", "ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(hdulist[0].header['CTYPE1'], hdulist[0].header['CTYPE2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can tell that this is in the FK5 coordinate system by the presence of a 'equinox' entry in the header:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hdulist[0].header['equinox']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Add some labels

\n", "
\n", "\n", "\n", "
\n", "\n", "

Now we have a nice plot with physically meaningful ticks, we should label our axes.

\n", "

Add labels to the axes saying \"Right Ascension [degrees]\" and \"Declination [degrees]\"

\n", "

Also overlay a coordinate grid using:\n", "ax.coords.grid()\n", "Look up the documentation for this method to see what parameters you can specify.

\n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Solution

\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "ax = plt.subplot(111, projection=wcs)\n", "ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower')\n", "ax.set_xlabel(\"Right Ascension [degrees]\")\n", "ax.set_ylabel(\"Declination [degrees]\")\n", "ax.coords.grid(color='white', alpha=0.5, linestyle='solid')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have a nice plot, we can do a couple of things to plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Overplotting in Pixel Coordinates" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "ax = plt.subplot(111, projection=wcs)\n", "\n", "ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower')\n", "ax.set_xlabel(\"Right Ascension [degrees]\")\n", "ax.set_ylabel(\"Declination [degrees]\")\n", "ax.coords.grid(color='white', alpha=0.5, linestyle='solid')\n", "\n", "ax.plot(3000, 3000, 'o')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Overplotting in World Coordinates" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111, projection=wcs)\n", "\n", "ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower')\n", "ax.set_xlabel(\"Right Ascension [degrees]\")\n", "ax.set_ylabel(\"Declination [degrees]\")\n", "ax.coords.grid(color='white', alpha=0.5, linestyle='solid')\n", "\n", "ax.set_autoscale_on(False)\n", "\n", "ax.plot(3000, 3000, 'o')\n", "# Overplot in FK5\n", "ax.plot_coord(SkyCoord(189.25*u.deg, 14.23*u.deg, frame=\"fk5\"), \"o\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Overplot in a different coordinate system

\n", "
\n", "\n", "\n", "
\n", "\n", "

Using ax.get_transform() you can specify coordinates in any system that astropy coordinates can transform to, try overplotting a point in the 'galactic' system.

\n", "

Also overlay a coordinate grid using:\n", "ax.coords.grid()\n", "Look up the documentation for this method to see what parameters you can specify.

\n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Overplotting Another Coordinate System" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111, projection=wcs)\n", "ax.imshow(hdulist[0].data, cmap='gray', vmax=1000, interpolation=None, origin='lower')\n", "ax.set_xlabel(\"Right Ascension [degrees]\")\n", "ax.set_ylabel(\"Declination [degrees]\")\n", "ax.coords.grid(color='white', alpha=0.5, linestyle='solid')\n", "\n", "overlay = ax.get_coords_overlay('galactic')\n", "overlay.grid(color='orange', alpha=1, linestyle='solid')\n", "overlay['l'].set_axislabel(\"Galactic Longitude [degrees]\")\n", "overlay['b'].set_axislabel(\"Galactic Latitude [degrees]\")\n", "\n", "ax.plot_coord(SkyCoord(287.5*u.deg, 76.65*u.deg, frame=\"galactic\"), \"o\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## SunPy Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The SunPy Map class is a wrapper for solar images which makes some of the above opertations easier." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sunpy.map\n", "from sunpy.data.sample import AIA_171_ROLL_IMAGE\n", "\n", "amap = sunpy.map.Map(AIA_171_ROLL_IMAGE)\n", "amap.peek()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Coordinate Systems\n", "\n", "The SunPy map will calculate some of the things we did earlier automatically, and add some extra information." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amap.coordinate_system" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amap.coordinate_frame" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amap.wcs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SunPy map uses these attributes to do it's plotting with WCSAxes like we did manually earlier. We can therefore use the overplotting and grid techniques to improve the plot:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import astropy.units as u" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", "amap = sunpy.map.Map(AIA_171_ROLL_IMAGE)\n", "\n", "im = amap.plot()\n", "ax = plt.gca()\n", "\n", "x = 500*u.arcsec\n", "y = -300*u.arcsec\n", "\n", "ax.plot_coord(SkyCoord(x, y, frame=amap.coordinate_frame), \"o\")\n", "\n", "overlay = ax.get_coords_overlay('heliographic_stonyhurst')\n", "overlay.grid(color='cyan', alpha=1, linestyle='solid')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amap.pixel_to_data(100*u.pix, 200*u.pix)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amap.data_to_pixel(SkyCoord(0*u.arcsec, 0*u.arcsec, frame=amap.coordinate_frame))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "
\n", "

Rotate your Owl

\n", "
\n", "\n", "\n", "
\n", "\n", "

Why is the Sun wonky?

\n", "

Use the rotate() method of SunPy Map to align the coordinate grid to the pixel grid in this sample image.

\n", "

Once you have run rotate, plot the resulting image, and compare with the one above.

\n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mr = amap.rotate()\n", "mr.peek()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" }, "name": "Astro_Images.ipynb" }, "nbformat": 4, "nbformat_minor": 1 }