VOBS (Virtual Observatories)#
Abstract: Geomagnetic Virtual Observatories (GVOs) spatially aggregate magnetic field measurements made by spacecraft at different times, in order to create time series of the magnetic field at 300 fixed locations (the virtual observatories). Data collected near each virtual observatory are combined by fitting to a local potential field in order to produce estimates of the observed field at monthly and four-monthly cadences. The resulting data products are suitable for studying Earth’s core dynamo, being comparable to classical ground observatories but with uniform global coverage. For more information, see the project webpage and Hammer, M.D., et al. Geomagnetic Virtual Observatories: monitoring geomagnetic secular variation with the Swarm satellites. Earth Planets Space 73, 54 (2021). https://doi.org/10.1186/s40623-021-01357-9
See also:
https://nbviewer.jupyter.org/github/pacesm/jupyter_notebooks/blob/master/VOBS/VOBS_data_access.ipynb
http://www.spacecenter.dk/files/magnetic-models/GVO/GVO_Product_Definition.pdf
(link to dashboard TBD)
# Show important version numbers to help debugging
%load_ext watermark
%watermark -i -v -p viresclient,pandas,xarray,matplotlib,numpy
Python implementation: CPython
Python version : 3.11.6
IPython version : 8.18.0
viresclient: 0.12.0
pandas : 2.1.3
xarray : 2023.12.0
matplotlib : 3.8.2
numpy : 1.26.2
import datetime as dt
import numpy as np
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from tqdm import tqdm
from viresclient import SwarmRequest
Quickstart#
Core field predictions ("B_CF" and associated uncertainty "sigma_CF") from the Swarm 1-monthly VOBS product, "SW_OPER_VOBS_1M_2_", can be fetched with:
from viresclient import SwarmRequest
import datetime as dt
request = SwarmRequest()
request.set_collection("SW_OPER_VOBS_1M_2_")
request.set_products(
measurements=["SiteCode", "B_CF", "sigma_CF"]
)
data = request.get_between(
dt.datetime(2000, 1, 1),
dt.datetime.now(),
asynchronous=False, # Process synchronously (suitable for small data)
show_progress=False # Disable the downloading progress bar
)
We can load the data as a pandas dataframe (use expand=True to split the vectors out into columns for each component):
df = data.as_dataframe(expand=True)
df.tail()
| SiteCode | Latitude | Radius | Longitude | B_CF_N | B_CF_E | B_CF_C | sigma_CF_N | sigma_CF_E | sigma_CF_C | |
|---|---|---|---|---|---|---|---|---|---|---|
| Timestamp | ||||||||||
| 2024-04-15 12:00:00 | S77W166 | -77.29173 | 6861200.0 | -165.52934 | -1680.25969 | 9433.38183 | -46080.55981 | 4.344905 | 5.277821 | 6.545478 |
| 2024-04-15 12:00:00 | S77W114 | -77.29173 | 6861200.0 | -114.10077 | 7213.01020 | 11112.62000 | -39943.53297 | 2.244145 | 5.203383 | 6.956875 |
| 2024-04-15 12:00:00 | S77W063 | -77.29173 | 6861200.0 | -62.67219 | 13678.82153 | 5082.40554 | -34070.67029 | 2.224223 | 3.206085 | 2.378352 |
| 2024-04-15 12:00:00 | S77W011 | -77.29173 | 6861200.0 | -11.24362 | 13768.52254 | -4026.03511 | -32268.98151 | 2.650189 | 3.981576 | 4.793556 |
| 2024-04-15 12:00:00 | S90E000 | -89.90000 | 6861200.0 | 0.00000 | 10181.00857 | -6910.78946 | -41192.08502 | 16.890712 | 10.515642 | 4.985520 |
Or as an xarray dataset (use reshape=True to get the data in a higher-dimensional form that more naturally accommodates the shape of the data):
ds = data.as_xarray(reshape=True)
ds
<xarray.Dataset>
Dimensions: (Timestamp: 125, SiteCode: 300, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2013-12-15T12:00:00 ... 2024-04-15T...
Latitude (SiteCode) float64 5.964 5.964 5.964 ... -77.29 -77.29 -89.9
Longitude (SiteCode) float64 1.844 13.46 25.07 36.68 ... -114.1 -165.5 0.0
Radius (SiteCode) float64 6.861e+06 6.861e+06 ... 6.861e+06 6.861e+06
* NEC (NEC) <U1 'N' 'E' 'C'
* SiteCode (SiteCode) <U7 'N06E002' 'N06E013' ... 'S77W166' 'S90E000'
Data variables:
B_CF (SiteCode, Timestamp, NEC) float64 nan nan ... -4.119e+04
sigma_CF (SiteCode, Timestamp, NEC) float64 nan nan nan ... 10.52 4.986
Attributes:
Sources: ['SW_OPER_VOBS_1M_2__20131215T000000_20240415T000000_0201']
MagneticModels: []
AppliedFilters: []The "SiteCode" variable is set as an index of the dataset (only when we specify reshape=True above) so we can quickly extract the data for a particular virtual observatory:
ds.sel(SiteCode="N90E000")
<xarray.Dataset>
Dimensions: (Timestamp: 125, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2013-12-15T12:00:00 ... 2024-04-15T...
Latitude float64 89.9
Longitude float64 0.0
Radius float64 6.861e+06
* NEC (NEC) <U1 'N' 'E' 'C'
SiteCode <U7 'N90E000'
Data variables:
B_CF (Timestamp, NEC) float64 nan nan nan ... 1.113e+03 1.35 4.619e+04
sigma_CF (Timestamp, NEC) float64 nan nan nan nan ... 8.171 6.033 8.298
Attributes:
Sources: ['SW_OPER_VOBS_1M_2__20131215T000000_20240415T000000_0201']
MagneticModels: []
AppliedFilters: []This enables rapid visualisation of the data from just that observatory, for example with:
ds.sel(SiteCode="N90E000").sel(NEC="N")["B_CF"].plot.line(x="Timestamp");
ds.sel(SiteCode="N90E000")["B_CF"].plot.line(x="Timestamp", col="NEC", sharey=False);
A quick graph of the virtual observatory locations, and how we can extract the list of their names in the "SiteCode" variable:
ds.isel(Timestamp=0).plot.scatter(y="Latitude", x="Longitude");
ds["SiteCode"].values[0:10]
array(['N06E002', 'N06E013', 'N06E025', 'N06E037', 'N06E048', 'N06E060',
'N06E072', 'N06E083', 'N06E095', 'N06E106'], dtype='<U7')
Product collection details#
There are five products available, available through VirES under collections of the same name:
Collection Name |
Description |
|---|---|
|
Swarm 1 month data from all virtual observatories |
|
Swarm 4 month data from all virtual observatories |
|
CHAMP 4 month data from all virtual observatories |
|
Cryosat-2 4 month data from all virtual observatories |
|
Ørsted 4 month data from all virtual observatories |
|
Composite (Ørsted, CHAMP, Cryosat-2, Swarm) 4 month data from all virtual observatories |
These collections each contain the variables:
Variable |
Unit |
Dimension |
Description |
|---|---|---|---|
|
$\(-\)$ |
char [7] |
virtual observatory identifier |
|
$\(-\)$ |
scalar |
UTC time of observation |
|
$\(\text{deg}\)$ |
scalar |
ITRF geocentric latitude |
|
$\(\text{deg}\)$ |
scalar |
ITRF geocentric longitude |
|
$\(\text{m}\)$ |
scalar |
ITRF geocentric radius |
|
$\(\text{nT}\)$ |
vector [3] |
Core magnetic field vector in ITRF NEC frame. |
|
$\(\text{nT}\)$ |
vector [3] |
Observed magnetic field vector in ITRF NEC frame. |
|
$\(\text{nT}\)$ |
vector [3] |
Estimated error of the core magnetic field vector in ITRF NEC frame. |
|
$\(\text{nT}\)$ |
vector [3] |
Estimated error of the observed magnetic field vector in ITRF NEC frame. |
The secular variation estimates are available within separate collections (because of the different time sampling points) with :SecularVariation appended:
Collection Name |
Description |
|---|---|
|
Swarm 1 month secular variation data from all virtual observatories |
|
Swarm 4 month secular variation data from all virtual observatories |
|
CHAMP 4 month secular variation data from all virtual observatories |
|
Cryosat-2 4 month secular variation data from all virtual observatories |
|
Ørsted 4 month secular variation data from all virtual observatories |
|
Composite (Ørsted, CHAMP, Cryosat-2, Swarm) 4 month secular variation data from all virtual observatories |
These collections similarly contain the variables:
Variable |
Unit |
Dimension |
Description |
|---|---|---|---|
|
$\(-\)$ |
char [7] |
virtual observatory identifier |
|
$\(-\)$ |
scalar |
UTC time of observation |
|
$\(\text{deg}\)$ |
scalar |
ITRF geocentric latitude |
|
$\(\text{deg}\)$ |
scalar |
ITRF geocentric longitude |
|
$\(\text{m}\)$ |
scalar |
ITRF geocentric radius |
|
$\(\text{nT}/\text{yr}\)$ |
vector [3] |
Field secular variation vector in ITRF NEC frame. |
|
$\(\text{nT}/\text{yr}\)$ |
vector [3] |
Estimated error of the field secular variation vector in ITRF NEC frame. |
Sub-collections are also defined for each virtual observatory, named according to SiteCode, so that data from a specific observatory can be fetched alone by specifying collections named like "SW_OPER_VOBS_1M_2_:N65W051" or "SW_OPER_VOBS_1M_2_:SecularVariation:N65W051".
NB: VirES provides the data in the NEC frame in order to be consistent with the other Swarm products available through VirES. This is in contrast to the source files which are provided in the (Radial, Theta, Phi) frame.
The list of available observatory names (i.e. SiteCode) can be queried with:
request = SwarmRequest()
request.available_observatories("SW_OPER_VOBS_1M_2_")
Magnetic model predictions can also be fetched directly just as with the MAGx_LR products (but it is currently not possible to directly fetch the secular variation predictions of models).
For example, we can fetch the data for a specific observatory, together with IGRF predictions:
request = SwarmRequest()
request.set_collection("SW_OPER_VOBS_1M_2_:N65W051")
request.set_products(
measurements=["SiteCode", "B_OB", "sigma_OB", "B_CF", "sigma_CF"],
models=["IGRF"]
)
data = request.get_between(
dt.datetime(2016, 1, 1),
dt.datetime(2017, 1, 1),
asynchronous=False,
show_progress=False
)
ds = data.as_xarray()
ds
<xarray.Dataset>
Dimensions: (Timestamp: 12, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2016-01-15T12:00:00 ... 2016-12-15...
* NEC (NEC) <U1 'N' 'E' 'C'
Data variables:
B_OB (Timestamp, NEC) float64 7.7e+03 -3.608e+03 ... 4.357e+04
Longitude (Timestamp) float64 -51.43 -51.43 -51.43 ... -51.43 -51.43
Radius (Timestamp) float64 6.861e+06 6.861e+06 ... 6.861e+06 6.861e+06
sigma_CF (Timestamp, NEC) float64 5.94 2.612 5.954 ... 5.94 2.612 5.954
Latitude (Timestamp) float64 65.26 65.26 65.26 ... 65.26 65.26 65.26
SiteCode (Timestamp) <U7 'N65W051' 'N65W051' ... 'N65W051' 'N65W051'
sigma_OB (Timestamp, NEC) float64 19.96 33.28 12.13 ... 23.93 29.31 11.13
B_CF (Timestamp, NEC) float64 7.713e+03 -3.606e+03 ... 4.355e+04
B_NEC_IGRF (Timestamp, NEC) float64 7.718e+03 -3.605e+03 ... 4.356e+04
Attributes:
Sources: ['SW_OPER_AUX_IGR_2__19000101T000000_20241231T235959_010...
MagneticModels: ['IGRF = IGRF(max_degree=13,min_degree=1)']
AppliedFilters: []A quick plot comparing the virtual observatory core field estimates with the IGRF predictions:
ds.plot.scatter(x="Timestamp", y="B_CF", col="NEC", sharey=False, figsize=(15,3))
ds.plot.scatter(x="Timestamp", y="B_NEC_IGRF", col="NEC", sharey=False, figsize=(15,3));
Fetching and merging all available data#
Here is an example of loading data from all the products. We merge secular variation (SV) into datasets containing also the field measurements by defining a Timestamp_SV variable, rotate from the NEC frame to the RTP (Radial, Theta, Phi) frame, and collect the five products into a dictionary of five items.
import datetime as dt
import numpy as np
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from tqdm import tqdm
from viresclient import SwarmRequest
COLLECTIONS = {
'Swarm_1M': 'SW_OPER_VOBS_1M_2_',
'Swarm_4M': 'SW_OPER_VOBS_4M_2_',
'CHAMP_1M': 'CH_OPER_VOBS_1M_2_',
'Orsted_1M': 'OR_OPER_VOBS_1M_2_',
'Cryosat_1M': 'CR_OPER_VOBS_1M_2_',
'Composite_1M': 'CO_OPER_VOBS_1M_2_',
'Orsted_4M': 'OR_OPER_VOBS_4M_2_',
'CHAMP_4M': 'CH_OPER_VOBS_4M_2_',
'Cryosat_4M': 'CR_OPER_VOBS_4M_2_',
'Composite_4M': 'CO_OPER_VOBS_4M_2_'
}
def nec2rtp(ds):
"""Convert NEC coords to RTP. NB: data_var names do not change"""
_ds = ds.copy()
# Convert variables from NEC to RTP frame
for var in _ds.data_vars:
if "NEC" in _ds[var].dims:
_ds[var] = np.array([-1, 1, -1])*_ds[var]
_ds[var] = _ds[var].roll({"NEC": 1}, roll_coords=False)
_ds[var].attrs = ds[var].attrs
# Rename NEC dims & coords to RTP
_ds = _ds.rename({"NEC": "RTP"})
_ds = _ds.assign_coords({"RTP": ["Radial", "Theta", "Phi"]})
_ds["RTP"].attrs = {
"units": "",
"description": "RTP frame - Radial, Theta, Phi [R,T,P] = [-C,-N,E]"
}
return _ds
def fetch_vobs(collection, sv=False, reshape=True, rtp=True):
"""Fetch data from VirES for a given collection"""
collection = f"{collection}:SecularVariation" if sv else collection
if sv:
measurements = ['SiteCode', 'B_SV', 'sigma_SV']
else:
measurements = ['SiteCode', 'B_CF', 'B_OB', 'sigma_CF', 'sigma_OB']
request = SwarmRequest()
request.set_collection(collection)
request.set_products(
measurements=measurements,
)
data = request.get_between(
dt.datetime(1999, 1, 1),
dt.datetime(2024, 1, 1),
asynchronous=False, show_progress=False
)
ds = data.as_xarray(reshape=reshape)
if rtp:
ds = nec2rtp(ds)
return ds
def merge_vobs(ds, ds_sv):
"""Merge in SecularVariation data by using new 'Timestamp_SV' coord"""
ds_sv = ds_sv.rename({"Timestamp": "Timestamp_SV"})
# Copy metadata
attrs = ds.attrs.copy()
for k, v in ds_sv.attrs.items():
attrs[k].extend(v)
attrs[k] = list(set(attrs[k]))
ds = xr.merge((ds, ds_sv))
ds.attrs = attrs
return ds
def fetch_all_vobs():
"""Gives a dictionary containing datasets, one for each VOBS collection"""
ALL_VOBS = {}
for key, collection in tqdm(COLLECTIONS.items()):
ds = fetch_vobs(collection)
ds_sv = fetch_vobs(collection, sv=True)
ds = merge_vobs(ds, ds_sv)
ALL_VOBS[key] = ds
return ALL_VOBS
ALL_VOBS = fetch_all_vobs()
0%| | 0/10 [00:00<?, ?it/s]
10%|█ | 1/10 [00:06<00:55, 6.21s/it]
20%|██ | 2/10 [00:11<00:45, 5.73s/it]
30%|███ | 3/10 [00:17<00:41, 5.96s/it]
40%|████ | 4/10 [00:23<00:35, 5.93s/it]
50%|█████ | 5/10 [00:29<00:29, 5.97s/it]
60%|██████ | 6/10 [00:36<00:25, 6.29s/it]
70%|███████ | 7/10 [00:41<00:17, 5.94s/it]
80%|████████ | 8/10 [00:47<00:11, 5.77s/it]
90%|█████████ | 9/10 [00:52<00:05, 5.62s/it]
100%|██████████| 10/10 [00:58<00:00, 5.64s/it]
100%|██████████| 10/10 [00:58<00:00, 5.83s/it]
Each dataset is now available within the dictionary ALL_VOBS, and can be accessed like:
ALL_VOBS["Swarm_4M"]
<xarray.Dataset>
Dimensions: (Timestamp: 30, SiteCode: 300, RTP: 3, Timestamp_SV: 28)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2014-03-01T12:00:00 ... 2023-11-...
Latitude (SiteCode) float64 5.964 5.964 5.964 ... -77.29 -77.29 -89.9
Longitude (SiteCode) float64 1.844 13.46 25.07 ... -114.1 -165.5 0.0
Radius (SiteCode) float64 6.861e+06 6.861e+06 ... 6.861e+06 6.861e+06
* SiteCode (SiteCode) <U7 'N06E002' 'N06E013' ... 'S77W166' 'S90E000'
* RTP (RTP) <U6 'Radial' 'Theta' 'Phi'
* Timestamp_SV (Timestamp_SV) datetime64[ns] 2014-09-15T12:00:00 ... 2023-...
Data variables:
B_CF (SiteCode, Timestamp, RTP) float64 5.286e+03 ... -6.887e+03
sigma_OB (SiteCode, Timestamp, RTP) float64 -0.9269 -2.192 ... 20.62
B_OB (SiteCode, Timestamp, RTP) float64 5.284e+03 ... -6.877e+03
sigma_CF (SiteCode, Timestamp, RTP) float64 -1.333 -2.317 ... 5.867
sigma_SV (SiteCode, Timestamp_SV, RTP) float64 -1.261 -0.7624 ... 7.043
B_SV (SiteCode, Timestamp_SV, RTP) float64 43.97 -10.39 ... -28.29
Attributes:
Sources: ['SW_OPER_VOBS_4M_2__20140301T000000_20240301T000000_0201']
MagneticModels: []
AppliedFilters: []Note that the altitudes of the virtual observatories defined for each mission are different, so the measurements can not directly be compared. Below, we can see that CHAMP and Swarm VO’s are at lower altitude so the observed field is stronger. The composite data product has VO’s defined at a common altitude of 700km.
rad = {}
for mission in ("Orsted", "CHAMP", "Cryosat", "Swarm", "Composite"):
rad[mission] = int(ALL_VOBS[f"{mission}_4M"]["Radius"][0]/1e3)
plt.figure(figsize=(10,3))
for mission in ("Orsted", "CHAMP", "Cryosat", "Swarm", "Composite"):
missionname = "Ørsted" if mission == "Orsted" else mission
label = f"{missionname} @ {rad[mission]}km"
ALL_VOBS[f"{mission}_4M"].sel(SiteCode="S77W114", RTP="Radial").plot.scatter(
x="Timestamp", y="B_OB", label=label
)
plt.title("Radial component at S77W114 (4-monthly estimate)")
plt.xlabel("Time")
plt.grid()
plt.legend(loc=(1.05, 0), title="Geocentric distance");
Visualising data from one location#
def plot_vobs_data(mission, site):
"""Plot a 3x3 grid of the data from a given mission at a given VO site"""
fig, axes = plt.subplots(nrows=3, ncols=3, sharex="col", figsize=(15, 10))
available_1m = True if f"{mission}_1M" in ALL_VOBS.keys() else False
if available_1m:
ds_1m = ALL_VOBS[f"{mission}_1M"].sel(SiteCode=site)
ds_4m = ALL_VOBS[f"{mission}_4M"].sel(SiteCode=site)
for i, rtp in enumerate(("Radial", "Theta", "Phi")):
if available_1m:
_ds_1m = ds_1m.sel(RTP=rtp)
# Observed field
axes[i, 0].errorbar(
_ds_1m["Timestamp"].values, _ds_1m["B_OB"].values, np.abs(_ds_1m["sigma_OB"].values),
fmt=".", label="1M"
)
# Core field
axes[i, 1].errorbar(
_ds_1m["Timestamp"].values, _ds_1m["B_CF"].values, np.abs(_ds_1m["sigma_CF"].values),
fmt=".",
)
# Secular variation (of core field)
axes[i, 2].errorbar(
_ds_1m["Timestamp_SV"].values, _ds_1m["B_SV"].values, np.abs(_ds_1m["sigma_SV"].values),
fmt=".",
)
_ds_4m = ds_4m.sel(RTP=rtp)
# Observed field
axes[i, 0].errorbar(
_ds_4m["Timestamp"].values, _ds_4m["B_OB"].values, np.abs(_ds_4m["sigma_OB"].values),
fmt=".", label="4M"
)
# Core field
axes[i, 1].errorbar(
_ds_4m["Timestamp"].values, _ds_4m["B_CF"].values, np.abs(_ds_4m["sigma_CF"].values),
fmt=".",
)
axes[i, 1].set_ylim(axes[i, 0].get_ylim())
# Secular variation (of core field)
axes[i, 2].errorbar(
_ds_4m["Timestamp_SV"].values, _ds_4m["B_SV"].values, np.abs(_ds_4m["sigma_SV"].values),
fmt=".",
)
axes[0, 0].set_ylabel("Radial\n[nT]")
axes[1, 0].set_ylabel("Theta\n[nT]")
axes[2, 0].set_ylabel("Phi\n[nT]")
for i in (0, 1, 2):
axes[i, 1].set_ylabel("[nT]")
axes[i, 2].set_ylabel("[nT/yr]")
axes[0, 0].set_title("B_OB (Observed field)")
axes[0, 1].set_title("B_CF (Core field)")
axes[0, 2].set_title("B_SV (Secular variation)")
axes[0, 0].legend()
fig.tight_layout()
fig.suptitle(site, va="bottom", y=1, fontsize=15)
return fig, axes
plot_vobs_data("Swarm", "N65W051");
Visualising data from all locations (using cartopy)#
LINE_COLORS = {
"Radial": mpl.colors.to_hex("tab:blue"),
"Theta": mpl.colors.to_hex("tab:orange"),
"Phi": mpl.colors.to_hex("tab:green"),
}
TITLES = {
"Radial": r"dB$_r$ / dt",
"Theta": r"dB$_\theta$ / dt",
"Phi": r"dB$_\phi$ / dt"
}
def make_sv_map(ds, RTP="Radial", var="B_SV"):
"""Generate overview map of SV from a given dataset"""
# Set up underlying map to plot onto
fig = plt.figure(figsize=(16, 8), dpi=150)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(),extent=[-180, 180, -90, 90])
ax.add_feature(cfeature.LAND, color="lightgrey")
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
GVO_LOCATIONS = np.vstack((ds["Longitude"].values, ds["Latitude"].values)).T
ax.scatter(*GVO_LOCATIONS.T, color="black", s=3, zorder=2)
# Add fractional year to use as x-axis
tvar = "Timestamp_SV" if var == "B_SV" else "Timestamp"
ds["Year"] = ds[tvar].dt.year + ds[tvar].dt.month/12
ds = ds.set_coords("Year")
# Min and max values to use for scaling
min_y = np.nanmin(ds[var].sel(RTP=RTP).data)
max_y = np.nanmax(ds[var].sel(RTP=RTP).data)
min_x = np.nanmin(ds["Year"].data)
max_x = np.nanmax(ds["Year"].data)
scale_x = 10
scale_y = 40
# Loop through each GVO
for i in range(300):
# Extract data for only that GVO and vector component
gvo_record = ds.sel(RTP=RTP).isel(SiteCode=i)
# Get x & y values and scale them, centred around the GVO location
# Scale values between 0 & 1:
x = (gvo_record["Year"].data - min_x) / (max_x - min_x)
y = (gvo_record[var].data - min_y) / (max_y - min_y)
# Shift values to centre around the lat/lon position:
lat = float(gvo_record["Latitude"])
lon = float(gvo_record["Longitude"])
x = lon + scale_x*(x - 0.5)
y = lat + scale_y*(y - np.nanmean(y))
# Plot these points onto the figure
gvo_xy_verts = np.vstack((x, y))
ax.scatter(
*gvo_xy_verts, transform=ccrs.PlateCarree(),
color=LINE_COLORS.get(RTP), alpha=0.8, s=1
)
# Create scale indicator
dx = scale_x
dy = scale_y * 20 / (max_y-min_y)
p_x = 160
p_y = 105
ax.arrow(p_x, p_y, dx, 0, linewidth=2, head_width=0).set_clip_on(False)
ax.arrow(p_x, p_y, 0, dy, linewidth=2, head_width=0).set_clip_on(False)
ax.text(p_x-2, p_y+dx/2, "20nT/yr", va="top", ha="right")
minyr = str(np.round(min_x, 1))
maxyr = str(np.round(max_x, 1))
ax.text(p_x, p_y-2, f"{minyr} - {maxyr}", va="top", ha="left")
fig.suptitle(TITLES[RTP], fontsize=15)
return fig, ax
fig, ax = make_sv_map(ALL_VOBS["Swarm_4M"], RTP="Radial")
You can now repeat this graph, changing the key to the ALL_VOBS dictionary to one of:
'Swarm_1M', 'Swarm_4M', 'CHAMP_1M', 'Orsted_1M', 'Cryosat_1M', 'Composite_1M', 'Orsted_4M', 'CHAMP_4M', 'Cryosat_4M', 'Composite_4M'
and RTP to one of 'Radial', 'Theta', 'Phi'. For example:
fig, ax = make_sv_map(ALL_VOBS["Composite_1M"], RTP="Phi")
