{
"cells": [
{
"cell_type": "markdown",
"id": "e175dfa5-b401-416b-b8b3-642382f6a8d4",
"metadata": {},
"source": [
"# Modelling with Swarm data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "731f9d08-0da4-4218-9ad6-d2b73ee197ef",
"metadata": {},
"outputs": [],
"source": [
"import datetime as dt\n",
"import numpy as np\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"import holoviews as hv\n",
"import hvplot.xarray\n",
"\n",
"from chaosmagpy.model_utils import design_gauss, synth_values, power_spectrum\n",
"from chaosmagpy.coordinate_utils import geo_to_gg, gg_to_geo\n",
"from chaosmagpy.plot_utils import plot_power_spectrum\n",
"from viresclient import SwarmRequest"
]
},
{
"cell_type": "markdown",
"id": "07e20d67-6deb-496c-9ee5-775e16403999",
"metadata": {},
"source": [
"## Fetch Swarm data to use"
]
},
{
"cell_type": "markdown",
"id": "3b04caa1-8dca-4f0a-8eb4-9fa331ce9489",
"metadata": {},
"source": [
"Following the data selection demonstrated in the previous notebook, but choosing a time period near the equinox:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f2af9967-baab-4e4e-9664-68cc8844378b",
"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(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 = fetch_selected_data()\n",
"ds"
]
},
{
"cell_type": "markdown",
"id": "15d64efe-4361-49cd-b4c4-b165721494d9",
"metadata": {},
"source": [
"## Constructing a model"
]
},
{
"cell_type": "markdown",
"id": "569356cf-7ca9-4cb7-a97b-4d2db208d8d4",
"metadata": {},
"source": [
"We can perform a least-squares inversion as we did before with the small number of ground observatory measurements. In this case, we can easily go to higher spherical harmonic degree because we have many more data and they are globally distributed.\n",
"\n",
"We'll do one thing differently: use only the radial component of the measurements."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6960a893-ae23-4f1d-81f0-49d742507561",
"metadata": {},
"outputs": [],
"source": [
"def build_model(ds, nmax=13):\n",
" \"\"\"Use the contents of dataset, ds, to build a SH model\"\"\"\n",
" # Get the positions and measurements from ds\n",
" radius = (ds[\"Radius\"]/1e3).values\n",
" theta = (90-ds[\"Latitude\"]).values\n",
" phi = (ds[\"Longitude\"]).values\n",
" B_radius = -ds[\"B_NEC\"].sel(NEC=\"C\").values\n",
" B_theta = -ds[\"B_NEC\"].sel(NEC=\"N\").values\n",
" B_phi = ds[\"B_NEC\"].sel(NEC=\"E\").values\n",
" # Use ChaosMagPy to build the design matrix, G\n",
" # https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.model_utils.design_gauss.html\n",
" G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax)\n",
" # Using only the radial component, perform the inversion:\n",
" G = np.concatenate((G_radius, ))\n",
" d = np.concatenate((B_radius, ))\n",
" m = np.linalg.inv(G.T @ G) @ (G.T @ d)\n",
" return m\n",
"\n",
"model_coeffs = build_model(ds, nmax=13)\n",
"print(f\"First 10 coefficients of the model:\\n{model_coeffs[:10]}\")"
]
},
{
"cell_type": "markdown",
"id": "d547a0a6-e2ae-4001-8b59-0756f60f3921",
"metadata": {},
"source": [
"Let's sample our model and visualise it with contours:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6f39740-54cb-4501-b4ae-adf44f7c6695",
"metadata": {},
"outputs": [],
"source": [
"def sample_model_on_grid(m, radius=6371.200):\n",
" \"\"\"Evaluate Gauss coefficients over a regular grid over Earth\"\"\"\n",
" theta = np.arange(1, 180, 1)\n",
" phi = np.arange(-180, 180, 1)\n",
" theta, phi = np.meshgrid(theta, phi)\n",
" B_r, B_t, B_p = synth_values(m, radius, theta, phi)\n",
" ds_grid = xr.Dataset(\n",
" data_vars={\"B_NEC_model\": ((\"y\", \"x\", \"NEC\"), np.stack((-B_t, B_p, -B_r), axis=2))},\n",
" coords={\"Latitude\": ((\"y\", \"x\"), 90-theta), \"Longitude\": ((\"y\", \"x\"), phi), \"NEC\": np.array([\"N\", \"E\", \"C\"])}\n",
" )\n",
" return ds_grid\n",
"\n",
"# Evaluate at the mean radius of the satellite measurements\n",
"mean_radius_km = float(ds[\"Radius\"].mean())/1e3\n",
"ds_grid = sample_model_on_grid(model_coeffs, radius=mean_radius_km)\n",
"# Generate contour plot of vertical component\n",
"ds_grid.sel(NEC=\"C\").hvplot.contourf(\n",
" x=\"Longitude\", y=\"Latitude\", c=\"B_NEC_model\",\n",
" levels=30, coastline=True, projection=ccrs.PlateCarree(), global_extent=True,\n",
" clabel=\"Model: vertical (downwards) magnetic field (nT)\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "e38a570a-4386-4ec4-8f2c-4e5d1862189a",
"metadata": {},
"source": [
"### Exercise"
]
},
{
"cell_type": "markdown",
"id": "49fa9d81-1303-4bfd-a8f6-77c2081b18c7",
"metadata": {},
"source": [
"1. We used `synth_values()` above to evaluate the model over the whole surface of Earth. Evaluate the model just in Lisbon (38.72° N, -9.14° W) and calculate the magnetic declination. \n",
" Hints:\n",
" - Use `synth_values(model_coeffs, ...)` to find the vector components ($B_r, B_\\theta, B_\\phi$)\n",
" - Refer to the [geomagnetic components](https://intermagnet.org/faq/10.geomagnetic-comp) to calculate D\n",
" - NB: ($B_r, B_\\theta, B_\\phi$) corresponds to (-Z, -X, Y). You may neglect the difference between geodetic and geocentric coordinates.\n",
"1. Plot the power spectrum of the model\n",
"1. Try pushing the inversion to higher spherical harmonic degree"
]
},
{
"cell_type": "markdown",
"id": "d5253e24-f63c-4f5d-9bea-d35526c5aa6b",
"metadata": {},
"source": [
"## Data-model residuals"
]
},
{
"cell_type": "markdown",
"id": "b02552fa-768a-43a6-bc3b-af8bbcba4633",
"metadata": {},
"source": [
"Now let's dig a little deeper, comparing the input data and the model...\n",
"\n",
"Let's plot the difference between the data and the model, the \"data-model residuals\":"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7ef550e0-4d83-424d-a420-ff8a6186df36",
"metadata": {},
"outputs": [],
"source": [
"# Evaluate the model at the same points as the measurements\n",
"def append_model_evaluations(ds, m):\n",
" ds = ds.copy()\n",
" B_r, B_t, B_p = synth_values(m, ds[\"Radius\"]/1e3, 90-ds[\"Latitude\"], ds[\"Longitude\"])\n",
" ds[\"B_NEC_model\"] = (\"Timestamp\", \"NEC\"), np.stack((-B_t, B_p, -B_r), axis=1)\n",
" ds[\"B_NEC_res_model\"] = ds[\"B_NEC\"] - ds[\"B_NEC_model\"]\n",
" return ds\n",
"ds = append_model_evaluations(ds, model_coeffs)\n",
"\n",
"# Plot the vertical component residual\n",
"_ds = ds.drop_vars(\"Timestamp\").sel(NEC=\"C\")\n",
"_ds.hvplot.scatter(\n",
" x=\"Longitude\", y=\"Latitude\", color=\"B_NEC_res_model\",\n",
" rasterize=True, x_sampling=2, y_sampling=2,\n",
" colorbar=True, cmap=\"coolwarm\", clim=(-40, 40), clabel=\"Vertical (B_C) data-model residuals (nT)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "45efa58f-a871-4329-91b3-eb59eb20c8e9",
"metadata": {},
"source": [
"(Above): The difference between the model and the data is quite small except over the poles.\n",
"\n",
"Let's look at the Eastward component:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "81ab6c67-8b12-408c-8a9a-7a9c76e7e194",
"metadata": {},
"outputs": [],
"source": [
"_ds = ds.drop_vars(\"Timestamp\").sel(NEC=\"E\")\n",
"_ds.hvplot.scatter(\n",
" x=\"Longitude\", y=\"Latitude\", color=\"B_NEC_res_model\",\n",
" rasterize=True, x_sampling=2, y_sampling=2,\n",
" colorbar=True, cmap=\"coolwarm\", clim=(-40, 40), clabel=\"Eastward (B_E) data-model residuals (nT)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "7df261a0-3b2b-4ce3-964e-27901b59b35a",
"metadata": {},
"source": [
"(Above): The disturbance over the poles is much stronger in the Eastward component... this is related to the geometry of the field-aligned currents producing a larger magnetic perturbation in that direction.\n",
"\n",
"Next let's plot against magnetic local time (MLT) instead, as the relevant coordinate system in which these currents are organised:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c614d8ef-3a39-434c-892c-e93641e9942a",
"metadata": {},
"outputs": [],
"source": [
"_ds = ds.drop_vars(\"Timestamp\").sel(NEC=\"E\")\n",
"_ds.hvplot.scatter(\n",
" x=\"MLT\", y=\"Latitude\", color=\"B_NEC_res_model\",\n",
" rasterize=True, x_sampling=2/15, y_sampling=1,\n",
" colorbar=True, cmap=\"coolwarm\", clim=(-40, 40), clabel=\"Eastward (B_E) data-model residuals (nT)\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c010873-3d3b-4866-8bd2-8a00b4180e35",
"metadata": {},
"outputs": [],
"source": [
"ds.drop_vars(\"Timestamp\").hvplot.scatter(x=\"Latitude\", y=\"B_NEC_res_model\", col=\"NEC\", rasterize=True, colorbar=False, cmap=\"bwr\")"
]
},
{
"cell_type": "markdown",
"id": "a09df52b-1122-4e84-a6b9-4f866c74e68e",
"metadata": {},
"source": [
"(Above and below): Both the Northward (N) and Eastward (E) components are more disturbed over the poles than the downward (C) component. This disturbance is due to Field-Aligned Currents (FACs) and Auroral Electrojets (AEJs) which produce strong magnetic fields over the poles."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "59882eb4-6ff5-4c3c-a2f0-3eae173d7150",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"(ds.groupby_bins(\"Latitude\", 90)\n",
" .apply(lambda x: x[\"B_NEC_res_model\"].std(axis=0))\n",
" .plot.line(x=\"Latitude_bins\", ax=ax)\n",
")\n",
"ax.set_title(\"Standard deviations\");"
]
},
{
"cell_type": "markdown",
"id": "7d7794eb-0ae5-4a85-9591-5f2039821ef7",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"1. Alter the inversion to use all the components (instead of only the vertical component), and inspect the resulting model by plotting the residuals using the code above.\n",
"\n",
" Hint: Change the design matrix, `G`, and data, `d` to: \n",
" ```\n",
" G = np.concatenate((G_radius, G_theta, G_phi))\n",
" d = np.concatenate((B_radius, B_theta, B_phi))\n",
" ```\n",
" \n",
" How is the resultant model affected? Why (and how) might we treat different components differently in the inversion?\n",
"\n",
"1. Plot histograms of the residuals"
]
},
{
"cell_type": "markdown",
"id": "e3059399-4887-48ec-936d-06a33abce4ee",
"metadata": {},
"source": [
"## Exploring further"
]
},
{
"cell_type": "markdown",
"id": "02b77a2b-1a7b-4769-8072-8e33d2f4f54a",
"metadata": {},
"source": [
"Take a look at [Swarm Notebooks](https://notebooks.vires.services/notebooks/04a1_geomag-models-vires) if you want more examples. There you will also find recipes for accessing other Swarm products."
]
}
],
"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
}