Modelling with Swarm data#
import datetime as dt
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import holoviews as hv
import hvplot.xarray
from chaosmagpy.model_utils import design_gauss, synth_values, power_spectrum
from chaosmagpy.coordinate_utils import geo_to_gg, gg_to_geo
from chaosmagpy.plot_utils import plot_power_spectrum
from viresclient import SwarmRequest
ERROR 1: PROJ: proj_create_from_database: Open of /opt/conda/share/proj failed
Fetch Swarm data to use#
Following the data selection demonstrated in the previous notebook, but choosing a time period near the equinox:
def fetch_selected_data():
"""Fetch Swarm data, using filters to select the "clean" measurements"""
request = SwarmRequest()
request.set_collection("SW_OPER_MAGA_LR_1B")
request.set_products(
measurements=["B_NEC"],
sampling_step="PT1M",
auxiliaries=["Kp", "dDst", "MLT", "SunZenithAngle"],
)
# Reject bad data according to quality flags
request.add_filter("Flags_B != 255")
request.add_filter("Flags_q != 255")
# Reject geomagnetically active times
request.add_filter("Kp <= 3")
request.add_filter("dDst <= 3")
# Reject sunlit data (retain data where Sun is 10° below horizon)
request.add_filter("SunZenithAngle >= 80")
data = request.get_between(
start_time=dt.datetime(2018, 9, 14),
end_time=dt.datetime(2018, 9, 28)
)
ds = data.as_xarray()
ds.attrs.pop("Sources")
return ds
ds = fetch_selected_data()
ds
<xarray.Dataset>
Dimensions: (Timestamp: 7502, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2018-09-14 ... 2018-09-27T22:2...
* NEC (NEC) <U1 'N' 'E' 'C'
Data variables:
Spacecraft (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A'
MLT (Timestamp) float64 22.08 22.43 22.72 ... 7.747 8.778 17.05
Kp (Timestamp) float64 3.0 3.0 3.0 3.0 ... 0.333 0.333 0.333
dDst (Timestamp) float64 2.0 2.0 2.0 2.0 2.0 ... 3.0 3.0 3.0 3.0
Longitude (Timestamp) float64 0.9939 3.279 4.676 ... -167.7 -145.2
Latitude (Timestamp) float64 -76.42 -72.64 -68.83 ... 79.32 -81.83
SunZenithAngle (Timestamp) float64 107.1 110.8 114.5 ... 85.62 81.77 80.2
Radius (Timestamp) float64 6.824e+06 6.824e+06 ... 6.823e+06
B_NEC (Timestamp, NEC) float64 1.326e+04 -6.031e+03 ... -4.455e+04
Attributes:
MagneticModels: []
AppliedFilters: ['Flags_B != 255', 'Flags_q != 255', 'Kp <= 3', 'SunZeni...Constructing a model#
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.
We’ll do one thing differently: use only the radial component of the measurements.
def build_model(ds, nmax=13):
"""Use the contents of dataset, ds, to build a SH model"""
# Get the positions and measurements from ds
radius = (ds["Radius"]/1e3).values
theta = (90-ds["Latitude"]).values
phi = (ds["Longitude"]).values
B_radius = -ds["B_NEC"].sel(NEC="C").values
B_theta = -ds["B_NEC"].sel(NEC="N").values
B_phi = ds["B_NEC"].sel(NEC="E").values
# Use ChaosMagPy to build the design matrix, G
# https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.model_utils.design_gauss.html
G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax)
# Using only the radial component, perform the inversion:
G = np.concatenate((G_radius, ))
d = np.concatenate((B_radius, ))
m = np.linalg.inv(G.T @ G) @ (G.T @ d)
return m
model_coeffs = build_model(ds, nmax=13)
print(f"First 10 coefficients of the model:\n{model_coeffs[:10]}")
First 10 coefficients of the model:
[-29426.5267484 -1461.20812388 4687.2639744 -2485.94283138
2991.3833521 -2952.8286684 1679.36197956 -705.07457418
1360.15221396 -2373.48205529]
Let’s sample our model and visualise it with contours:
def sample_model_on_grid(m, radius=6371.200):
"""Evaluate Gauss coefficients over a regular grid over Earth"""
theta = np.arange(1, 180, 1)
phi = np.arange(-180, 180, 1)
theta, phi = np.meshgrid(theta, phi)
B_r, B_t, B_p = synth_values(m, radius, theta, phi)
ds_grid = xr.Dataset(
data_vars={"B_NEC_model": (("y", "x", "NEC"), np.stack((-B_t, B_p, -B_r), axis=2))},
coords={"Latitude": (("y", "x"), 90-theta), "Longitude": (("y", "x"), phi), "NEC": np.array(["N", "E", "C"])}
)
return ds_grid
# Evaluate at the mean radius of the satellite measurements
mean_radius_km = float(ds["Radius"].mean())/1e3
ds_grid = sample_model_on_grid(model_coeffs, radius=mean_radius_km)
# Generate contour plot of vertical component
ds_grid.sel(NEC="C").hvplot.contourf(
x="Longitude", y="Latitude", c="B_NEC_model",
levels=30, coastline=True, projection=ccrs.PlateCarree(), global_extent=True,
clabel="Model: vertical (downwards) magnetic field (nT)"
)
Exercise#
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.
Hints:Use
synth_values(model_coeffs, ...)to find the vector components (Refer to the geomagnetic components to calculate D
NB: (
Plot the power spectrum of the model
Try pushing the inversion to higher spherical harmonic degree
Data-model residuals#
Now let’s dig a little deeper, comparing the input data and the model…
Let’s plot the difference between the data and the model, the “data-model residuals”:
# Evaluate the model at the same points as the measurements
def append_model_evaluations(ds, m):
ds = ds.copy()
B_r, B_t, B_p = synth_values(m, ds["Radius"]/1e3, 90-ds["Latitude"], ds["Longitude"])
ds["B_NEC_model"] = ("Timestamp", "NEC"), np.stack((-B_t, B_p, -B_r), axis=1)
ds["B_NEC_res_model"] = ds["B_NEC"] - ds["B_NEC_model"]
return ds
ds = append_model_evaluations(ds, model_coeffs)
# Plot the vertical component residual
_ds = ds.drop_vars("Timestamp").sel(NEC="C")
_ds.hvplot.scatter(
x="Longitude", y="Latitude", color="B_NEC_res_model",
rasterize=True, x_sampling=2, y_sampling=2,
colorbar=True, cmap="coolwarm", clim=(-40, 40), clabel="Vertical (B_C) data-model residuals (nT)",
)
(Above): The difference between the model and the data is quite small except over the poles.
Let’s look at the Eastward component:
_ds = ds.drop_vars("Timestamp").sel(NEC="E")
_ds.hvplot.scatter(
x="Longitude", y="Latitude", color="B_NEC_res_model",
rasterize=True, x_sampling=2, y_sampling=2,
colorbar=True, cmap="coolwarm", clim=(-40, 40), clabel="Eastward (B_E) data-model residuals (nT)",
)
(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.
Next let’s plot against magnetic local time (MLT) instead, as the relevant coordinate system in which these currents are organised:
_ds = ds.drop_vars("Timestamp").sel(NEC="E")
_ds.hvplot.scatter(
x="MLT", y="Latitude", color="B_NEC_res_model",
rasterize=True, x_sampling=2/15, y_sampling=1,
colorbar=True, cmap="coolwarm", clim=(-40, 40), clabel="Eastward (B_E) data-model residuals (nT)",
)
ds.drop_vars("Timestamp").hvplot.scatter(x="Latitude", y="B_NEC_res_model", col="NEC", rasterize=True, colorbar=False, cmap="bwr")
(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.
fig, ax = plt.subplots(1, 1)
(ds.groupby_bins("Latitude", 90)
.apply(lambda x: x["B_NEC_res_model"].std(axis=0))
.plot.line(x="Latitude_bins", ax=ax)
)
ax.set_title("Standard deviations");
Exercise#
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.
Hint: Change the design matrix,
G, and data,dto:G = np.concatenate((G_radius, G_theta, G_phi)) d = np.concatenate((B_radius, B_theta, B_phi))
How is the resultant model affected? Why (and how) might we treat different components differently in the inversion?
Plot histograms of the residuals
Exploring further#
Take a look at Swarm Notebooks if you want more examples. There you will also find recipes for accessing other Swarm products.