{
"cells": [
{
"cell_type": "markdown",
"id": "8bc9a443-dec8-44a4-bec7-4d3ebc2c9655",
"metadata": {},
"source": [
"# Satellite data selection"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0f6eec56-3613-431a-ad6d-f736329d078a",
"metadata": {},
"outputs": [],
"source": [
"import datetime as dt\n",
"import numpy as np\n",
"import xarray as xr\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"\n",
"from viresclient import SwarmRequest"
]
},
{
"cell_type": "markdown",
"id": "1c4595be-0709-41c2-ac8c-3056fbf352e3",
"metadata": {},
"source": [
"## Fetching Swarm data using VirES"
]
},
{
"cell_type": "markdown",
"id": "15190996-9744-4ec1-a033-c909a14de83e",
"metadata": {},
"source": [
"Magnetic field measurements from Swarm Alpha are available within the [SW_OPER_MAGA_LR_1B](https://swarmhandbook.earth.esa.int/catalogue/sw_magx_lr_1b) dataset. Here we use the [viresclient](https://viresclient.readthedocs.io/) Python package to fetch these measurements."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a9436ac-f4b0-4f1a-b271-77a5c9f15500",
"metadata": {},
"outputs": [],
"source": [
"def fetch_data():\n",
" \"\"\"Fetch Swarm data using the VirES service\"\"\"\n",
" \n",
" request = SwarmRequest()\n",
" # Specify dataset to use\n",
" request.set_collection(\"SW_OPER_MAGA_LR_1B\")\n",
" request.set_products(\n",
" measurements=[\"B_NEC\"],\n",
" sampling_step=\"PT1M\", # Sample at 1-minute resolution\n",
" )\n",
" # Fetch 1 month of data\n",
" data = request.get_between(\n",
" start_time=dt.datetime(2025,1,1),\n",
" end_time=dt.datetime(2025,2,1)\n",
" )\n",
" # Load data using xarray\n",
" ds = data.as_xarray()\n",
" ds.attrs.pop(\"Sources\")\n",
" return ds\n",
"\n",
"ds_1 = fetch_data()\n",
"ds_1"
]
},
{
"cell_type": "markdown",
"id": "b32e48b8-af96-4b04-b482-b07f3a8dd61e",
"metadata": {},
"source": [
"## Inspecting data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d1524390-79f3-40a4-9160-ef9039105695",
"metadata": {},
"outputs": [],
"source": [
"def plot_data_preview(ds):\n",
" \"\"\"Plot magnetic field vector against latitude and time\"\"\"\n",
"\n",
" # Configure figure layout\n",
" fig = plt.figure(figsize=(20, 5))\n",
" gs = fig.add_gridspec(3, 2)\n",
" ax_L0 = fig.add_subplot(gs[:, 0])\n",
" ax_R0 = fig.add_subplot(gs[0, 1])\n",
" ax_R1 = fig.add_subplot(gs[1, 1])\n",
" ax_R2 = fig.add_subplot(gs[2, 1])\n",
" # Left: Plot B_NEC vector components against latitude, on one figure\n",
" ax_L0.scatter(x=ds[\"Latitude\"], y=ds[\"B_NEC\"].sel(NEC=\"N\"), s=.1, label=\"$B_N$ (North)\")\n",
" ax_L0.scatter(x=ds[\"Latitude\"], y=ds[\"B_NEC\"].sel(NEC=\"E\"), s=.1, label=\"$B_E$ (East)\")\n",
" ax_L0.scatter(x=ds[\"Latitude\"], y=ds[\"B_NEC\"].sel(NEC=\"C\"), s=.1, label=\"$B_C$ (Centre)\")\n",
" ax_L0.set_ylabel(\"$B_{NEC}$ [nT]\")\n",
" ax_L0.set_xlabel(\"Latitude [deg]\")\n",
" ax_L0.grid()\n",
" ax_L0.legend(markerscale=10)\n",
" # Right: Plot B_NEC against time, over three figures\n",
" ax_R0.scatter(x=ds[\"Timestamp\"], y=ds[\"B_NEC\"].sel(NEC=\"N\"), s=.1, c=\"tab:blue\")\n",
" ax_R1.scatter(x=ds[\"Timestamp\"], y=ds[\"B_NEC\"].sel(NEC=\"E\"), s=.1, c=\"tab:green\")\n",
" ax_R2.scatter(x=ds[\"Timestamp\"], y=ds[\"B_NEC\"].sel(NEC=\"C\"), s=.1, c=\"tab:orange\")\n",
" # Adjust axes etc\n",
" fig.subplots_adjust(wspace=.1, hspace=0)\n",
" ax_R0.set_xticklabels([])\n",
" ax_R1.set_xticklabels([])\n",
" plt.setp(ax_R2.get_xticklabels(), rotation=45, ha='right')\n",
" N_samples = len(ds[\"B_NEC\"])\n",
" fig.suptitle(f\"Exploratory plots of vector field components vs latitude and time\\nNumber of data = {N_samples}\");\n",
"\n",
" return fig\n",
"\n",
"fig_1 = plot_data_preview(ds_1)"
]
},
{
"cell_type": "markdown",
"id": "ba276b8e-b19a-48b7-a648-eb17c78b4822",
"metadata": {},
"source": [
"## Selecting data to use for modelling"
]
},
{
"cell_type": "markdown",
"id": "44eaa9c6-99bf-4867-bd03-d354834f79a0",
"metadata": {},
"source": [
"Here we apply three types of data selection:\n",
"\n",
"- Data quality: Instrumental problems can occur and should(!) be indicated by the quality flags (`Flags_B`, `Flags_F`, `Flags_q`)\n",
"- Geomagnetic activity: Select quiet periods based on Kp and rate of change of Dst\n",
"- Sunlight: The dayside ionosphere is more active so we choose to use only nightside data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c6789204-d38e-418d-9c8b-0d49875ecdf2",
"metadata": {},
"outputs": [],
"source": [
"def fetch_selected_data():\n",
" \"\"\"Fetch Swarm data, using filters to select the \"clean\" measurements\"\"\"\n",
" \n",
" request = SwarmRequest()\n",
" request.set_collection(\"SW_OPER_MAGA_LR_1B\")\n",
" request.set_products(\n",
" measurements=[\"B_NEC\"],\n",
" sampling_step=\"PT1M\",\n",
" # auxiliaries=[\"Kp\", \"dDst\", \"MLT\", \"SunZenithAngle\"],\n",
" )\n",
" # Reject bad data according to quality flags\n",
" request.add_filter(\"Flags_B != 255\")\n",
" request.add_filter(\"Flags_q != 255\")\n",
" # Reject geomagnetically active times\n",
" request.add_filter(\"Kp <= 3\")\n",
" request.add_filter(\"dDst <= 3\")\n",
" # Reject sunlit data (retain data where Sun is 10° below horizon)\n",
" request.add_filter(\"SunZenithAngle >= 80\")\n",
" data = request.get_between(\n",
" start_time=dt.datetime(2025,1,1),\n",
" end_time=dt.datetime(2025,2,1)\n",
" # start_time=dt.datetime(2018, 9, 14),\n",
" # end_time=dt.datetime(2018, 9, 28)\n",
" )\n",
" ds = data.as_xarray()\n",
" ds.attrs.pop(\"Sources\")\n",
" return ds\n",
"\n",
"ds_2 = fetch_selected_data()\n",
"ds_2"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0cbe2715-fd97-45c1-ba07-cf3628b98acd",
"metadata": {},
"outputs": [],
"source": [
"fig_2 = plot_data_preview(ds_2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f184d3d7-cabd-4cf7-a73d-956652240a90",
"metadata": {},
"outputs": [],
"source": [
"def global_plot(ds):\n",
" \"\"\"Plot intensity over the glob\"\"\"\n",
" \n",
" # Append intensity, F, to the dataset\n",
" ds[\"F\"] = np.sqrt((ds[\"B_NEC\"]**2).sum(axis=1))\n",
" \n",
" fig = plt.figure(figsize=(10, 10))\n",
" gs = fig.add_gridspec(2, 2)\n",
" ax_N = fig.add_subplot(gs[1, 0], projection=ccrs.Orthographic(0,90))\n",
" ax_S = fig.add_subplot(gs[1, 1], projection=ccrs.Orthographic(180,-90))\n",
" ax_G = fig.add_subplot(gs[0, :], projection=ccrs.EqualEarth())\n",
" \n",
" for ax in (ax_N, ax_S, ax_G):\n",
" pc = ax.scatter(\n",
" ds[\"Longitude\"], y=ds[\"Latitude\"], c=ds[\"F\"],\n",
" s=5, edgecolors=\"none\", cmap=\"inferno\", transform = ccrs.PlateCarree()\n",
" )\n",
" ax.coastlines(linewidth=1.1)\n",
" # Add colorbar inset axes into global map and move upwards\n",
" cax = inset_axes(ax_G, width=\"2%\", height=\"100%\", loc=\"right\", borderpad=-5)\n",
" clb = plt.colorbar(pc, cax=cax)\n",
" clb.set_label(\"Field intensity [nT]\", rotation=270, va=\"bottom\")\n",
"\n",
"fig_3 = global_plot(ds_2)"
]
},
{
"cell_type": "markdown",
"id": "5189f2fe-7d8a-44a6-9a6c-a83f39e162eb",
"metadata": {},
"source": [
"## Exercises\n",
"\n",
"1. Experiment with changing the selection criteria for Kp, dDst, solar zenith angle\n",
" (hint: alter the line `request.add_filter(\"Kp <= 3\")`)\n",
"1. Experiment using different months of data\n",
"1. How might the selection criteria bias the measurements?\n",
"1. Change the time period to:\n",
" ```\n",
" start_time=dt.datetime(2018, 9, 14),\n",
" end_time=dt.datetime(2018, 9, 28)\n",
" ```\n",
" (which we will use in the following notebook)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}