{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conjunctions"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"In this advanced tutorial we will combine a low Earth orbiting satellite ephemeris with ASIs, to calculate the mean auroral intensity at the satellite footprint as a function of time.\n",
"\n",
"The mean auroral intensity is calculated two ways: 1) the nearest pixel, and 2) a (20x20) km area at 110 km emission altitude."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"asilib version: 0.23.2\n"
]
}
],
"source": [
"from datetime import datetime, timedelta\n",
"from IPython.display import Video\n",
"import numpy as np\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import asilib\n",
"import asilib.asi\n",
"import asilib.map\n",
"\n",
"plt.style.use('dark_background')\n",
"\n",
"print(f'asilib version: {asilib.__version__}')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"location_code = 'RANK'\n",
"alt=110\n",
"time_range = (datetime(2017, 9, 15, 2, 32, 0), datetime(2017, 9, 15, 2, 35, 0))\n",
"asi = asilib.asi.themis(location_code, time_range=time_range, alt=alt)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Satellite footprint\n",
"Now we define an orbit path of a low Earth orbiting satellite (i.e. its footprint). This is a north-south satellite track oriented to the east of the THEMIS/RANK imager. In this context, `lla` stands for the `(latitude, longitude, altitude)` coordinates."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"n = int((time_range[1] - time_range[0]).total_seconds() / 3) # 3 second cadence.\n",
"lats = np.linspace(asi.meta[\"lat\"] + 5, asi.meta[\"lat\"] - 5, n)\n",
"lons = (asi.meta[\"lon\"] - 0.5) * np.ones(n)\n",
"alts = alt * np.ones(n)\n",
"sat_lla = np.array([lats, lons, alts]).T\n",
"sat_time = asi.data.time"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Create an `asilib.Conjunction()` object that handles mapping between the satellite and the imager. It takes in an `Imager` instance and arrays of the satellite times and LLA coordinates."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"conjunction_obj = asilib.Conjunction(asi, (sat_time, sat_lla))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, here are two steps that we are ignoring that you'll likely need to implement:\n",
"\n",
"1. Map the satellite's LLA coordinates along the magnetic field line from the satellite altitude down to 110 km (or whatever you chose for the `alt` kwarg.) This is done via `Conjunction.lla_footprint()` that requires the IRBEM library. IRBEM can be hard to install; in the future, I plan to change remove IRBEM in favor of geopack (or a similar package).\n",
"2. Normally the satellite LLA time stamps are not the same as the ASI. In that case you will need to call `Conjunction.interp_sat()` to interpolate the LLA coordinates to the ASI timestamps. Note: this method does not handle interpolation well across the anti-meridian (-180/180 degree longitude). If it detects that you're interpolating over it, it will issue a warning. \n",
"\n",
"### Nearest pixel intensity"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"sat_azel, sat_azel_pixels = conjunction_obj.map_azel()\n",
"nearest_pixel_intensity = conjunction_obj.intensity(box=None)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean pixel intensity in a 20x20 km area. \n",
"\n",
"The mean intensity is calculated with a masked array. It contains `np.nan` outside the 20x20 km area, and 1s inside. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"area_intensity = conjunction_obj.intensity(box=(10, 10))\n",
"\n",
"# You don't need to calculate the area mask if you just need the intensity, but this is useful if you \n",
"# want to animate and visualize the area\n",
"area_mask = conjunction_obj.equal_area(box=(10,10))\n",
"# Need to change masked NaNs to 0s so we can plot the rectangular area contours.\n",
"area_mask[np.where(np.isnan(area_mask))] = 0"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code block contains many steps.\n",
"\n",
"1. We create three subplots and initialize the animation generator.\n",
"2. We loop over each image.\n",
" - In the first subplot plot the entire satellite footprint using the `sat_azel_pixels` array.\n",
" - In the first subplot plot the instantaneous footprint. \n",
" - In the first subplot plot the 20x20 km area contour.\n",
" - In the second and third subplots plot the auroral intensity from a) the nearest pixel, and b) the 20x20 km area.\n",
" - In the second and third subplots plot a vertical line at the current time.\n",
" - Annotate the first subplot."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"20170915_023200_023500_themis_rank_fisheye.mp4: |###################### | 97%\n",
"Animation saved to C:\\Users\\shumkms1\\asilib-data\\animations\\20170915_023200_023500_themis_rank_fisheye.mp4\n"
]
}
],
"source": [
"fig, ax = plt.subplots(\n",
" 3, 1, figsize=(7, 10), gridspec_kw={'height_ratios': [4, 1, 1]}, constrained_layout=True\n",
")\n",
"ax[1].set(ylabel='ASI intensity\\nnearest pixel [counts]')\n",
"ax[2].set(xlabel='Time', ylabel='ASI intensity\\n10x10 km area [counts]')\n",
"\n",
"gen = asi.animate_fisheye_gen(\n",
" ax=ax[0], azel_contours=True, overwrite=True, cardinal_directions='news'\n",
")\n",
"\n",
"for i, (time, image, _, im) in enumerate(gen):\n",
" # Plot the entire satellite track, its current location, and a 20x20 km box\n",
" # around its location.\n",
" ax[0].plot(sat_azel_pixels[:, 0], sat_azel_pixels[:, 1], 'red')\n",
" ax[0].scatter(sat_azel_pixels[i, 0], sat_azel_pixels[i, 1], c='red', marker='o', s=50)\n",
" ax[0].contour(area_mask[i, :, :], levels=[0.99], colors=['yellow'])\n",
"\n",
" if 'vline1' in locals():\n",
" vline1.remove()\n",
" vline2.remove() \n",
" text_obj.remove()\n",
" else:\n",
" # Plot the ASI intensity along the satellite path\n",
" ax[1].plot(sat_time, nearest_pixel_intensity)\n",
" ax[2].plot(sat_time, area_intensity)\n",
" vline1 = ax[1].axvline(time, c='b')\n",
" vline2 = ax[2].axvline(time, c='b')\n",
"\n",
" # Annotate the location_code and satellite info in the top-left corner.\n",
" location_code_str = (\n",
" f'THEMIS/{location_code} '\n",
" f'LLA=({asi.meta[\"lat\"]:.2f}, '\n",
" f'{asi.meta[\"lon\"]:.2f}, {asi.meta[\"alt\"]:.2f})'\n",
" )\n",
" satellite_str = f'Satellite LLA=({sat_lla[i, 0]:.2f}, {sat_lla[i, 1]:.2f}, {sat_lla[i, 2]:.2f})'\n",
" text_obj = ax[0].text(\n",
" 0,\n",
" 1,\n",
" location_code_str + '\\n' + satellite_str,\n",
" va='top',\n",
" transform=ax[0].transAxes,\n",
" color='red',\n",
" )\n",
"plt.close()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# When you run this, you should see the video below in your asilib-data/movies directory.\n",
"Video('https://github.com/mshumko/asilib/raw/main/docs/_static/example_outputs/20170915_023200_023500_themis_rank_fisheye.mp4')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.12 ('asilib_test')",
"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.11.3"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "b9163b3a154006dba32002ac69546f1602817992f65b3d9c15b81834e8ed58d8"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}