OverplottingViSPFOV-ARD-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "710792f5",
   "metadata": {},
   "source": [
    "# Showing the Field of View of VISP on AIA / HMI #\n",
    "\n",
    "<div class=\"alert-info\">This example requires sunpy>=6.1 for the sunpy.visualization.extent function.</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "62536102",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dkist\n",
    "import dkist.net\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.units as u\n",
    "from astropy.time import Time\n",
    "from astropy.visualization import ImageNormalize, SqrtStretch\n",
    "\n",
    "import sunpy.map\n",
    "from sunpy.net import Fido, attrs as a\n",
    "from sunpy.visualization import drawing\n",
    "import os\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9343a36b",
   "metadata": {},
   "source": [
    "## Obtaining some data\n",
    "\n",
    "Let's load the data with dkist.load_dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fb635ebd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c4185b3dffa9406fa5751d5beb367f06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "33406d4182c64cb4923e001932875f6a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VISP_L1_20240903T190933_AKLQK_metadata.asdf:   0%|          | 0.00/1.79M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "86c85365e8cc467a85bf1fcd1cfb3a39",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBI_L1_20240903T190930_QZXUJQ_metadata.asdf:   0%|          | 0.00/4.99M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# results = Fido.search(a.dkist.Dataset('AJZRR'))\n",
    "results = Fido.search((a.dkist.Dataset('AKLQK')) | (a.dkist.Dataset('QZXUJQ')))\n",
    "files = Fido.fetch(results, path=\"~/DKISTDownloads/{primary_proposal_id}/{dataset_id}/\")\n",
    "\n",
    "ds = dkist.load_dataset(files)\n",
    "vbi_ds = ds[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c299935-4af5-4e42-b0be-da84ae8969cd",
   "metadata": {},
   "source": [
    "Need to download the VBI file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "78a8cbd7-7699-4351-9c9e-cf69382ff80d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Refreshing dataset inventory for dataset QZXUJQ [dkist.io.file_manager]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd90913db1364f34a1e426b3f38c3b94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "|                                                                                                             …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PENDING: Starting Transfer\n",
      "Task completed with SUCCEEDED status.\n"
     ]
    }
   ],
   "source": [
    "vbi_ds[:11].files.download()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7d5d8c7",
   "metadata": {},
   "source": [
    "Extract the two pixel axes which are correlated to the spatial dimensions from the ViSP dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "19876730",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<dkist.dataset.dataset.Dataset object at 0x7fb48453a650>\n",
       "This VISP Dataset consists of 20 frames.\n",
       "Files are stored in /home/stuart/DKISTDownloads/pid_2_70/AKLQK\n",
       "\n",
       "This calibration has Dataset ID AKLQK.\n",
       "The unique identifier for the input observe frames (Product ID) is L1-RFJNL.\n",
       "\n",
       "This Dataset has 2 pixel and 3 world dimensions.\n",
       "\n",
       "The data are represented by a <class 'dask.array.core.Array'> object:\n",
       "dask.array<getitem, shape=(20, 2555), dtype=float32, chunksize=(1, 2555), chunktype=numpy.ndarray>\n",
       "\n",
       "Array Dim  Axis Name                Data size  Bounds\n",
       "        0  raster scan step number         20  None\n",
       "        1  spatial along slit            2555  None\n",
       "\n",
       "World Dim  Axis Name                  Physical Type                   Units\n",
       "        2  time                       time                            s\n",
       "        1  helioprojective latitude   custom:pos.helioprojective.lat  arcsec\n",
       "        0  helioprojective longitude  custom:pos.helioprojective.lon  arcsec\n",
       "\n",
       "Correlation between pixel and world axes:\n",
       "\n",
       "                          |  PIXEL DIMENSIONS\n",
       "                          | spatial |  raster\n",
       "                          |  along  |   scan\n",
       "                          |   slit  |   step\n",
       "         WORLD DIMENSIONS |         |  number\n",
       "------------------------- | ------- | -------\n",
       "helioprojective longitude |    x    |    x\n",
       " helioprojective latitude |    x    |    x\n",
       "                     time |         |    x"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spatial_slice = ds[0][0, 0, :, 0, :]\n",
    "spatial_slice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bff01477",
   "metadata": {},
   "source": [
    "Extract the times of all the raster positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f1642568",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Time object: scale='utc' format='isot' value=['2024-09-03T19:09:33.922' '2024-09-03T19:10:03.081'\n",
       " '2024-09-03T19:10:32.240' '2024-09-03T19:11:01.400'\n",
       " '2024-09-03T19:11:30.559' '2024-09-03T19:11:59.718'\n",
       " '2024-09-03T19:12:28.877' '2024-09-03T19:12:58.037'\n",
       " '2024-09-03T19:13:27.196' '2024-09-03T19:13:56.355'\n",
       " '2024-09-03T19:14:25.514' '2024-09-03T19:14:54.674'\n",
       " '2024-09-03T19:15:23.833' '2024-09-03T19:15:52.992'\n",
       " '2024-09-03T19:16:22.151' '2024-09-03T19:16:51.311'\n",
       " '2024-09-03T19:17:20.470' '2024-09-03T19:17:49.629'\n",
       " '2024-09-03T19:18:18.788' '2024-09-03T19:18:47.948']>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "times = spatial_slice.axis_world_coords(\"time\")[0]\n",
    "times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77e35c66",
   "metadata": {},
   "source": [
    "And then find and download the closest AIA and HMI images to the time of the first raster position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "53a088b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e513685e57264867afe0c017535e04ba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "06bf96146c674464aadbcfccf4e11c11",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "hmi.m_45s.2024.09.03_19_10_30_TAI.magnetogram.fits:   0%|          | 0.00/16.3M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "df466c19892a436687efee2aca50af1f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "aia.lev1.171A_2024_09_03T19_09_33.35Z.image_lev1.fits:   0%|          | 0.00/12.2M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#results = Fido.search(a.Instrument.aia, a.Wavelength(171*u.AA), a.Time(times[0], times[-1], times[0]))\n",
    "query = Fido.search(\n",
    "    a.Time(times[0], times[-1], times[0]),\n",
    "    (a.Instrument(\"HMI\") & a.Physobs(\"los_magnetic_field\")) |\n",
    "    (a.Instrument(\"AIA\") & a.Wavelength(171 * u.angstrom))\n",
    ")\n",
    "files = Fido.fetch(query, site=\"NSO\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1733c4f0",
   "metadata": {},
   "source": [
    "Now we load the downloaded files into a sunpy.map objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9c6f86a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = sunpy.map.Map(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7627c716",
   "metadata": {},
   "outputs": [],
   "source": [
    "hmi = maps[0].rotate(order=3)\n",
    "aia = maps[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "510f0075",
   "metadata": {},
   "source": [
    "Finally we can use a combination of sunpy.map.GenericMap.plot and sunpy.visualization.drawing.extent to plot the extent of each tile on the disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "4815e845-bf0c-4496-a934-68e2c64defed",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "from astropy.visualization import AsinhStretch, ImageNormalize, PercentileInterval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "b6077bcf-280e-4405-b7d9-bdcc9f32f09e",
   "metadata": {},
   "outputs": [],
   "source": [
    "vbi_0_norm = ImageNormalize(vbi_ds[0].data, stretch=AsinhStretch(), interval=PercentileInterval(99))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "1e0c51cb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/stuart/.virtualenvs/dkist/lib/python3.13/site-packages/erfa/core.py:16909: RuntimeWarning: invalid value encountered in taiutc\n",
      "  utc1, utc2, c_retval = ufunc.taiutc(tai1, tai2)\n",
      "/home/stuart/.virtualenvs/dkist/lib/python3.13/site-packages/erfa/core.py:133: ErfaWarning: ERFA function \"taiutc\" yielded 2556 of \"dubious year (Note 4)\"\n",
      "  warn(f'ERFA function \"{func_name}\" yielded {wmsg}', ErfaWarning)\n",
      "/home/stuart/.virtualenvs/dkist/lib/python3.13/site-packages/erfa/core.py:16909: RuntimeWarning: invalid value encountered in taiutc\n",
      "  utc1, utc2, c_retval = ufunc.taiutc(tai1, tai2)\n",
      "/home/stuart/.virtualenvs/dkist/lib/python3.13/site-packages/erfa/core.py:133: ErfaWarning: ERFA function \"taiutc\" yielded 2556 of \"dubious year (Note 4)\"\n",
      "  warn(f'ERFA function \"{func_name}\" yielded {wmsg}', ErfaWarning)\n"
     ]
    },
    {
     "data": {