This notebook shows an example of rotating a feature across the solar disk and mapping it to an image coordinate system.

feature-rot-to-pix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2da73849-a330-4676-8107-67fdc9a65da1",
   "metadata": {},
   "source": [
    "# Rotating Features across the Solar Disk\n",
    "\n",
    "This notebook shows an example of a loop rotating across the solar disk and how the position at the loop apex maps to an image coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f34c1fd2-b68f-44a4-87ac-76a991c4d665",
   "metadata": {},
   "outputs": [],
   "source": [
    "from synthesizAR.models.geometry import semi_circular_loop\n",
    "from sunpy.map import make_fitswcs_header, Map\n",
    "from sunpy.coordinates import propagate_with_solar_surface, Helioprojective\n",
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.visualization import quantity_support, time_support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d1676e2-45b5-4a6f-ad10-3d32204c6905",
   "metadata": {},
   "outputs": [],
   "source": [
    "loop = semi_circular_loop(\n",
    "    250*u.Mm,\n",
    "    observer=SkyCoord(lon=0*u.deg, lat=20*u.deg, frame='heliographic_stonyhurst', obstime='2020-01-01'),\n",
    "    gamma=90*u.deg,\n",
    "    n_points=100,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac0be0d6-398a-4573-b2a1-f9e92f321b3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.arange(-13.5, 13.5, 0.5) * u.day + loop.obstime\n",
    "observer = SkyCoord(0*u.deg, 0*u.deg, 1*u.AU, frame='heliographic_stonyhurst', obstime=loop.obstime)\n",
    "shape = (500,500)\n",
    "scale = [5,5]*u.arcsec/u.pix\n",
    "frames = [Helioprojective(obstime=t, observer=observer) for t in times]\n",
    "headers = [make_fitswcs_header(shape, SkyCoord(0,0,unit='arcsec', frame=f), scale=scale, rotation_angle=25*u.deg) for f in frames]\n",
    "data = np.nan*np.ones(shape)\n",
    "maps = [Map(data, h) for h in headers]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f565e425-3f44-48d5-87ea-d474269b4e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "projected_coords = []\n",
    "with propagate_with_solar_surface(rotation_model='snodgrass'):\n",
    "    for m in maps:\n",
    "        projected_coords.append(loop.transform_to(m.coordinate_frame))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ec208c3-e351-4a96-b26b-31a1c456b3f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_coords = u.Quantity([u.Quantity(m.world_to_pixel(c[50]) )for c, m in zip(projected_coords,maps)])\n",
    "visibility = np.array([c[50].is_visible() for c in projected_coords])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d26dc0b2-8c75-4a2e-a5fe-0c7a882af0ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_coords_visible = np.where(np.tile(visibility[:,np.newaxis], 2), pixel_coords, np.nan)\n",
    "pixel_coords_invisible = np.where(~np.tile(visibility[:,np.newaxis], 2), pixel_coords, np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "416f81d1-622c-4df5-9b22-e10f3a5b910c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for frame_index in range(len(maps)):\n",
    "    fig = plt.figure(figsize=(10,5))\n",
    "    ax1 = fig.add_subplot(121, projection=maps[frame_index])\n",
    "    maps[frame_index].plot(axes=ax1)\n",
    "    maps[frame_index].draw_grid(axes=ax1,color='k')\n",
    "    \n",
    "    with propagate_with_solar_surface():\n",
    "        is_visible = loop.transform_to(maps[frame_index].coordinate_frame).is_visible()\n",
    "        if is_visible.any():\n",
    "            ax1.plot_coord(loop[is_visible], color='k')\n",
    "            ax1.plot_coord(loop[50], marker='.', color='C3', label='Loop Apex')\n",
    "            ax1.legend(frameon=False, loc=1)\n",
    "    \n",
    "    ax2 = fig.add_subplot(122)\n",
    "    with quantity_support():\n",
    "        with time_support():\n",
    "            ax2.plot(times[:frame_index+1], pixel_coords_visible[:frame_index+1,0], color='C0', ls='-', label='x, visible')\n",
    "            ax2.plot(times[:frame_index+1], pixel_coords_invisible[:frame_index+1,0], color='C0', ls='--', label='x, behind limb')\n",
    "            ax2.plot(times[:frame_index+1], pixel_coords_visible[:frame_index+1,1], color='C1', ls='-', label='y, visible')\n",
    "            ax2.plot(times[:frame_index+1], pixel_coords_invisible[:frame_index+1,1], color='C1', ls='--', label='y, behind limb')\n",
    "            ax2.set_xlim(times[[0,-1]])\n",
    "    \n",
    "    ax2.set_ylim(0,500)\n",
    "    ax2.set_ylabel('pixel position')\n",
    "    ax2.set_title('Loop Apex Position')\n",
    "    ax2.legend(loc=4, frameon=False)\n",
    "    fig.savefig(f'loop-rotation-animation/frame_{frame_index:06d}')\n",
    "    plt.close(fig=fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "decbbb35-31e9-4edb-b992-93fb446d63f4",
   "metadata": {},
   "source": [
    "To make an animation, then run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e122e173-31bd-4ff5-97e5-72994fdce2f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "!ffmpeg -framerate 15 -pattern_type glob -i 'loop-rotation-animation/*.png' -c:v libx264 -pix_fmt yuv420p loop-rotation-animation/animation.mp4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "761d366b-3147-4d44-bf3f-27d25fefd5aa",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:qpp-synthesis]",
   "language": "python",
   "name": "conda-env-qpp-synthesis-py"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}