DDT_2023_co-observe.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "60a84fd2-f352-49a8-be20-ba6266908a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, EarthLocation\n",
    "from astropy.time import Time\n",
    "\n",
    "import dkist\n",
    "import dkist.net\n",
    "import sunpy_soar\n",
    "import sunpy.map\n",
    "import sunraster.instr.spice\n",
    "from sunpy.net import Fido, attrs as a\n",
    "from sunkit_instruments import iris"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a56dc10f-3b27-4195-b953-c3b3413d205f",
   "metadata": {},
   "source": [
    "## Download some Data\n",
    "\n",
    "First we do a big multi-instrument Fido search and download a selection of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "61bcf0d7-9be6-4462-a33e-c07c35e836a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = Path(\"/data/sunpy/ddt_2023\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dc69755d-731d-4486-b38f-2f056edb6705",
   "metadata": {},
   "outputs": [],
   "source": [
    "dkist_query = a.Instrument.visp | a.Instrument.vbi\n",
    "solo_query = ((a.soar.Product('EUI-FSI174-IMAGE') | a.soar.Product('SPICE-N-RAS')) & a.Level(2))\n",
    "hmi_query = (a.Instrument.hmi & a.Physobs.intensity)\n",
    "iris_query = a.Instrument.iris"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7f714a21-37d5-46c4-add4-9f957c01f93e",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_res = Fido.search(\n",
    "    a.Time(\"2023/10/16 18:10\", \"2023/10/17 00:15\"),\n",
    "    dkist_query | solo_query | hmi_query | iris_query,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c6359a97-7ab2-467e-8478-1d8342789ffd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cbeed794cd36422da8cb489d3f357068",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Files Downloaded:   0%|          | 0/6 [00:00<?, ?file/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "solo_L2_spice-n-ras_20231016T182011_V01_218104099-071.fits:   0%|          | 0.00/116M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "solo_L2_eui-fsi174-image_20231016T181055250_V01.fits:   0%|          | 0.00/5.31M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VISP_L1_20231016T220247_ALDLJ.asdf:   0%|          | 0.00/9.08M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBI_L1_20231016T181613_AZYQR.asdf:   0%|          | 0.00/845k [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "hmi_ic_45s_2023_10_16_18_11_15_tai_continuum.fits:   0%|          | 0.00/15.7M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "iris_l2_20231016_180036_3460259102_sji_1330_t000_fits.gz:   0%|          | 0.00/30.5M [00:00<?, ?B/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<parfive.results.Results object at 0x7f7c9a9e7f10>\n",
       "['/data/sunpy/ddt_2023/VISP/VISP_L1_20231016T220247_ALDLJ.asdf', '/data/sunpy/ddt_2023/VBI/VBI_L1_20231016T181613_AZYQR.asdf', '/data/sunpy/ddt_2023/EUI/solo_L2_eui-fsi174-image_20231016T181055250_V01.fits', '/data/sunpy/ddt_2023/SPICE/solo_L2_spice-n-ras_20231016T182011_V01_218104099-071.fits', '/data/sunpy/ddt_2023/HMI/hmi_ic_45s_2023_10_16_18_11_15_tai_continuum.fits', '/data/sunpy/ddt_2023/IRIS/iris_l2_20231016_180036_3460259102_sji_1330_t000_fits.gz']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fido.fetch(all_res[:, 0], path=data_dir / \"{instrument}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e2639bea-9b02-4a02-aee0-bdcf31fe500b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;38;5;63m/data/sunpy/ddt_2023\u001b[0m\n",
      "├── \u001b[1;38;5;63mEUI\u001b[0m\n",
      "│   └── \u001b[1;38;5;40msolo_L2_eui-fsi174-image_20231016T181055250_V01.fits\u001b[0m\n",
      "├── \u001b[1;38;5;63mHMI\u001b[0m\n",
      "│   └── \u001b[1;38;5;40mhmi_ic_45s_2023_10_16_18_11_15_tai_continuum.fits\u001b[0m\n",
      "├── \u001b[1;38;5;63mIRIS\u001b[0m\n",
      "│   └── \u001b[1;38;5;40miris_l2_20231016_180036_3460259102_sji_1330_t000_fits.gz\u001b[0m\n",
      "├── \u001b[1;38;5;63mSPICE\u001b[0m\n",
      "│   └── \u001b[1;38;5;40msolo_L2_spice-n-ras_20231016T182011_V01_218104099-071.fits\u001b[0m\n",
      "├── \u001b[1;38;5;63mVBI\u001b[0m\n",
      "│   └── \u001b[1;38;5;40mVBI_L1_20231016T181613_AZYQR.asdf\u001b[0m\n",
      "└── \u001b[1;38;5;63mVISP\u001b[0m\n",
      "    └── \u001b[1;38;5;40mVISP_L1_20231016T220247_ALDLJ.asdf\u001b[0m\n",
      "\n",
      "7 directories, 6 files\n"
     ]
    }
   ],
   "source": [
    "!tree /data/sunpy/ddt_2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af25e74e-7465-46aa-a409-5dd32ee1959f",
   "metadata": {},
   "source": [
    "## Load some Data\n",
    "\n",
    "Now we load a bunch of this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d0a965e4-72a4-4360-98e8-673301506dec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# HMI for context\n",
    "m_hmi = sunpy.map.Map(data_dir / \"HMI\").rotate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f86537c5-d38f-4bca-9e4f-ead530fcc447",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: CROTA = 6.01664515285 / [deg] S/C counter-clockwise roll rel to Solar N\n",
      "keyword looks very much like CROTAn but isn't. [astropy.wcs.wcs]\n",
      "WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJDREF to 60233.764017 from DATEREF.\n",
      "Set MJD-OBS to 60233.764017 from DATE-OBS.\n",
      "Set MJD-BEG to 60233.764017 from DATE-BEG.\n",
      "Set MJD-AVG to 60233.767943 from DATE-AVG.\n",
      "Set MJD-END to 60233.771865 from DATE-END'. [astropy.wcs.wcs]\n",
      "WARNING: UnitsWarning: 'W/m2/sr/nm' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]\n"
     ]
    }
   ],
   "source": [
    "# Turn one part of the SPICE data into a map to extract a bounding box later.\n",
    "raster_spice = sunraster.instr.spice.read_spice_l2_fits(list((data_dir / \"SPICE\").glob(\"*\")))\n",
    "spice_window = raster_spice['N IV 765 - SH - Comp 8 ... Ne VIII 770 - LH - Comp 8 (Merged)']\n",
    "m_spice = sunpy.map.Map(spice_window[0,52,:,:].data, spice_window[0,52,:,:].meta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fd5e25a5-6497-4bc5-b945-84553a1a6325",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load all the steps of a SJI sequence into a MapSequence\n",
    "iris_maps = iris.SJI_to_sequence(list((data_dir / \"IRIS\").glob(\"*\"))[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "202f1f68-a6a9-4cfb-ba4a-87369dc8d86d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dkist_fido_convert_bb(res):\n",
    "    \"\"\"\n",
    "    Convert the BoundingBox column of the Fido results into a SkyCoord\n",
    "    \"\"\"\n",
    "    res = res.copy()\n",
    "    coords = np.array(list(map(literal_eval, res[\"Bounding Box\"])))\n",
    "    observer = EarthLocation.of_site(\"DKIST\").get_itrs(res[\"Start Time\"])\n",
    "    observer = observer.transform_to(sunpy.coordinates.HeliographicStonyhurst(obstime=res[\"Start Time\"]))\n",
    "    res[\"Bounding Box\"] = SkyCoord(coords[..., 0][:, ::-1], coords[..., 1][:, ::-1], unit=u.arcsec, frame=\"helioprojective\", obstime=res[\"Start Time\"][:, None], observer=observer[:, None])\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "777e1834-1f4f-4782-aade-24b329d2ab78",
   "metadata": {},
   "outputs": [],
   "source": [
    "visp_res, vbi_res = all_res[\"dkist\"]\n",
    "visp_res = dkist_fido_convert_bb(visp_res)\n",
    "vbi_res = dkist_fido_convert_bb(vbi_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f995f95e-f6e0-454d-afb2-7f8113311a6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Zoom in on the context image\n",
    "bl_hmi = SkyCoord(-400*u.arcsec, -180*u.arcsec, frame=m_hmi.coordinate_frame)\n",
    "tr_hmi = SkyCoord(100*u.arcsec, 240*u.arcsec, frame=m_hmi.coordinate_frame)\n",
    "hmi_sub = m_hmi.submap(bl_hmi, top_right=tr_hmi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "87b3db50-5e52-4bb4-a77d-0c2c6e8d1606",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.\n",
      "For frame 'heliographic_stonyhurst' the following metadata is missing: hglt_obs,hgln_obs\n",
      "For frame 'heliographic_carrington' the following metadata is missing: crln_obs,crlt_obs\n",
      " [sunpy.map.mapbase]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Bounding Boxes of all DKIST datasets from pid_2_114 (and some other instruments)')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {