2019_spd/sunpy-1.0-spd2019-coordinates.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.visualization import quantity_support\n",
    "quantity_support()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from astropy.time import Time\n",
    "from astropy.coordinates import get_body_barycentric, SkyCoord\n",
    "import astropy.units as u\n",
    "\n",
    "from sunpy.coordinates import get_body_heliographic_stonyhurst\n",
    "import sunpy.map\n",
    "\n",
    "import astropy.constants as const\n",
    "import scipy\n",
    "from sunpy.coordinates import frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = (16,8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!rm ~sunpy/data/secchi_l0_a_img_euvi_20100819_20100819_000715_n4eua.fts\n",
    "#!rm ~/sunpy/data/aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# SunPy & Coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SunPy coordinates allows us to deal with points in physical space, and the SunPy coordinates subpackage provides definitions of and transformations between several reference frames commonly used in solar physics. This allows us to manipulate Maps and take advantage of WCS Axes for plotting\n",
    "\n",
    "These reference frames and their associated transformations are implemented using the `astropy.coordinates` subpackage and extend Astropy’s coordinate frame transformation graph to include solar coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# STEREO & SDO loop \n",
    "In this example we will see how powerful the coordinates functionality in SunPy is for performing SDO-STEREO data analysis \n",
    "\n",
    "Lets first download some AIA and STEREO data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sunpy.net import Fido, attrs as a\n",
    "\n",
    "\n",
    "stereo = (a.vso.Source('STEREO_A') &\n",
    "          a.Instrument('EUVI') &\n",
    "          a.Time('2010-08-19', '2010-08-19T00:10:00'))\n",
    "\n",
    "aia = (a.Instrument('AIA') &\n",
    "       a.vso.Sample(24 * u.hour) &\n",
    "       a.Time('2010-08-19', '2010-08-19T00:10:00'))\n",
    "\n",
    "wave = a.Wavelength(17 * u.nm, 18 * u.nm)\n",
    "\n",
    "\n",
    "res = Fido.search(wave, aia | stereo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "files = Fido.fetch(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now lets make maps for both AIA and EUVI "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_aia, map_stereo = sunpy.map.Map(sorted(files))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Compare the Locations of STEREO-A and AIA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun_pos = get_body_heliographic_stonyhurst('sun', time = map_aia.date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(projection='polar')\n",
    "\n",
    "plt.polar(map_stereo.observer_coordinate.lon.to(u.rad), map_stereo.observer_coordinate.radius.to(u.AU), \n",
    "          marker='o',  ms=10, label='STEREO_A')\n",
    "plt.polar(map_aia.observer_coordinate.lon.to(u.rad), map_aia.observer_coordinate.radius.to(u.AU), \n",
    "          marker='o', ms=10, label='AIA')\n",
    "\n",
    "plt.polar(sun_pos.lon.to(u.rad), sun_pos.radius.to(u.AU), 'o', ms=20, label='Sun', color='yellow')\n",
    "ax.set_theta_zero_location(\"S\")\n",
    "\n",
    "plt.title('Position of the Sun, AIA and Stereo')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And plot the two images together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(1,2,1, projection = map_aia)\n",
    "map_aia.plot(axes = ax1)\n",
    "ax2 = fig.add_subplot(1,2,2, projection = map_stereo)\n",
    "map_stereo.plot(axes = ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lets now focus on the active region that is at ~ (500'', 480'') in AIA that is also seen in EUVI\n",
    "To do this we can define a box in AIA and then transform the corners of the box to the coordinate system of EUVI and overplot a rectange on both. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bl_aia = SkyCoord(350*u.arcsec, 300*u.arcsec, frame=map_aia.coordinate_frame)\n",
    "width = 300*u.arcsec\n",
    "height = 280*u.arcsec\n",
    "corners_aia = ([bl_aia.Tx, bl_aia.Ty],\n",
    "               [bl_aia.Tx + width, bl_aia.Ty],\n",
    "               [bl_aia.Ty, bl_aia.Ty + height],\n",
    "               [bl_aia.Tx + width, bl_aia.Ty + height])\n",
    "\n",
    "corners_aia = SkyCoord(corners_aia, frame=map_aia.coordinate_frame)\n",
    "\n",
    "#the coordinates of the corners now in STEREO/EUVI coordinate system\n",
    "corners_stereo = corners_aia.transform_to(map_stereo.coordinate_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(1,2,1, projection=map_aia)\n",
    "map_aia.plot(axes=ax)\n",
    "map_aia.draw_rectangle(corners_aia[0], \n",
    "                       corners_aia[1].Tx - corners_aia[0].Tx, \n",
    "                       corners_aia[2].Ty - corners_aia[0].Ty)\n",
    "\n",
    "ax2 = plt.subplot(1,2,2, projection=map_stereo)\n",
    "map_stereo.plot(axes=ax2)\n",
    "map_stereo.draw_rectangle(corners_stereo[0], \n",
    "                          corners_stereo[1].Tx - corners_stereo[0].Tx, \n",
    "                          corners_stereo[2].Ty - corners_stereo[0].Ty)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lets define some coordinates of a semicirclar loop of the active region and then plot this in both fields of view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "@u.quantity_input\n",
    "def semi_circular_loop(length: u.m, lon: u.deg=0*u.deg, lat: u.deg=0*u.deg):\n",
    "    \"\"\"\n",
    "    Return a Heliographic Stonyhurst coordinate object with points of a semi circular loop in it.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    length : `~astropy.Quantity`\n",
    "        Distance between footpoints of loops\n",
    "    lon : ~astropy.Quantity`\n",
    "        Longitude value of center of the loop\n",
    "    lat : `~astropy.Quantity`\n",
    "        Latitude value of center of loop\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    `astopy.coordinate.SkyCoord`\n",
    "        A SkyCoord of the coordinates of a semi-circular loop\n",
    "    \n",
    "    \"\"\"\n",
    "    r_sun = const.R_sun\n",
    "\n",
    "    def r_2_func(x):\n",
    "        return np.arccos(0.5 * x / r_sun.to(u.cm).value) - np.pi + length.to(u.cm).value / 2. / x\n",
    "\n",
    "    r_2 = scipy.optimize.bisect(r_2_func,\n",
    "                                length.to(u.cm).value / (2 * np.pi),\n",
    "                                length.to(u.cm).value / np.pi) * u.cm\n",
    "\n",
    "    alpha = np.arccos(0.5 * (r_2 / r_sun).decompose())\n",
    "    phi = np.linspace(-np.pi * u.rad + alpha, np.pi * u.rad - alpha, 2000)\n",
    "\n",
    "    # Quadratic formula to find r\n",
    "    a = 1.\n",
    "    b = -2 * (r_sun.to(u.cm) * np.cos(phi.to(u.radian)))\n",
    "    c = r_sun.to(u.cm)**2 - r_2.to(u.cm)**2\n",
    "    r = (-b + np.sqrt(b**2 - 4 * a * c)) / 2 / a\n",
    "    # Choose only points above the surface\n",
    "    i_r = np.where(r > r_sun)\n",
    "    r = r[i_r]\n",
    "    phi = phi[i_r]\n",
    "    hcc_frame = frames.Heliocentric(\n",
    "        observer=SkyCoord(lon=lon, lat=lat, radius=r_sun, frame='heliographic_stonyhurst'))\n",
    "\n",
    "    return SkyCoord(\n",
    "        x=r.to(u.cm) * np.sin(phi.to(u.radian)),\n",
    "        y=u.Quantity(r.shape[0] * [0 * u.cm]),\n",
    "        z=r.to(u.cm) * np.cos(phi.to(u.radian)),\n",
    "        frame=hcc_frame).transform_to('heliographic_stonyhurst')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can create a loop structure that fits over the active region of choice and plot it on AIA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we are defining the loop coordinates\n",
    "loop_coords = semi_circular_loop(200*u.Mm, lon=38*u.deg, lat=31*u.deg)\n",
    "\n",
    "# Make a map the size of our rectangle\n",
    "sub_aia = map_aia.submap(corners_aia[0], corners_aia[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax = plt.subplot(1,2,1, projection=map_aia)\n",
    "map_aia.plot(axes=ax)\n",
    "map_aia.draw_rectangle(corners_aia[0], \n",
    "                       corners_aia[1].Tx - corners_aia[0].Tx, \n",
    "                       corners_aia[2].Ty - corners_aia[0].Ty)\n",
    "\n",
    "ax.plot_coord(loop_coords.transform_to(map_aia.coordinate_frame), color='r', lw=2)\n",
    "\n",
    "ax1 = plt.subplot(1,2,2, projection=sub_aia)\n",
    "\n",
    "sub_aia.plot(axes=ax1)\n",
    "ax1.plot_coord(loop_coords.transform_to(sub_aia.coordinate_frame), color='r', lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plot the Loop on Both submaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_stereo = map_stereo.submap(corners_stereo[0], corners_stereo[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ax = plt.subplot(1,2,1, projection = sub_aia)\n",
    "sub_aia.plot(axes = ax)\n",
    "ax.plot_coord(loop_coords.transform_to(sub_aia.coordinate_frame), color = 'red', lw = 2)\n",
    "\n",
    "\n",
    "ax2 = plt.subplot(1,2,2, projection = sub_stereo)\n",
    "sub_stereo.plot(axes = ax2)\n",
    "ax2.plot_coord(loop_coords.transform_to(sub_stereo.coordinate_frame), color = 'red', lw = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Extract the Intensity Values along the Loop "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.int64(u.Quantity(map_aia.world_to_pixel(loop_coords)))\n",
    "\n",
    "aia_loop_int = map_aia.data[y, x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.int64(u.Quantity(map_stereo.world_to_pixel(loop_coords)))\n",
    "\n",
    "stereo_loop_int = map_stereo.data[y, x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance = np.cumsum(loop_coords[:-1].separation_3d(loop_coords[1:])).insert(0, 0*u.km).to(u.Mm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "l_stereo, = plt.plot(distance, stereo_loop_int)\n",
    "plt.ylabel(\"Intensity (STEREO)\")\n",
    "plt.xlabel(\"Distance Along Loop [Mm]\")\n",
    "plt.twinx()\n",
    "\n",
    "l_aia, = plt.plot(distance, aia_loop_int, color=\"C1\")\n",
    "plt.ylabel(\"Intensity (AIA)\")\n",
    "\n",
    "plt.legend((l_stereo, l_aia), (\"STEREO\", \"AIA\"))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "virtualenv_sunpy-dev",
   "language": "python",
   "name": "virtualenv_sunpy-dev"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}