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:

# 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.11.6
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()
Latitude Spacecraft Longitude Radius SiteCode B_CF_N B_CF_E B_CF_C sigma_CF_N sigma_CF_E sigma_CF_C
Timestamp
2023-10-15 12:00:00 -77.29173 - -165.52934 6861200.0 S77W166 -1686.36131 9417.76708 -46116.62875 4.514341 5.423621 6.418425
2023-10-15 12:00:00 -77.29173 - -114.10077 6861200.0 S77W114 7202.88941 11125.37982 -39981.35583 2.145221 5.339294 7.018290
2023-10-15 12:00:00 -77.29173 - -62.67219 6861200.0 S77W063 13678.82870 5104.79529 -34094.71856 2.472257 3.297114 2.456658
2023-10-15 12:00:00 -77.29173 - -11.24362 6861200.0 S77W011 13768.71708 -4025.04369 -32282.14367 2.650041 4.129388 4.460318
2023-10-15 12:00:00 -89.90000 - 0.00000 6861200.0 S90E000 10212.20632 -6870.00789 -41204.01398 17.570136 9.868229 5.075469

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: 119, SiteCode: 300, NEC: 3)
Coordinates:
  * Timestamp  (Timestamp) datetime64[ns] 2013-12-15T12:00:00 ... 2023-10-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.12e+04
    sigma_CF   (SiteCode, Timestamp, NEC) float64 nan nan nan ... 9.868 5.075
Attributes:
    Sources:         ['SW_OPER_VOBS_1M_2__20131215T000000_20231015T000000_0101']
    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: 119, NEC: 3)
Coordinates:
  * Timestamp  (Timestamp) datetime64[ns] 2013-12-15T12:00:00 ... 2023-10-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 ... -18.38 4.619e+04
    sigma_CF   (Timestamp, NEC) float64 nan nan nan nan ... 8.658 5.88 8.826
Attributes:
    Sources:         ['SW_OPER_VOBS_1M_2__20131215T000000_20231015T000000_0101']
    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");
../_images/116ab8120785ed677bfe3992e9d45c3d0aa1d0135333e99ec335b7eefc35633d.png
ds.sel(SiteCode="N90E000")["B_CF"].plot.line(x="Timestamp", col="NEC", sharey=False);
../_images/720839a4ea0d394df625972f79b7c50bfc3357490cbbfa7c1344fed9d73c6219.png

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");
../_images/3452f3f5e1b017ba783ca6740413ab84f9163c8deba1785026c1804ca728328b.png
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

SW_OPER_VOBS_1M_2_

Swarm 1 month data from all virtual observatories

SW_OPER_VOBS_4M_2_

Swarm 4 month data from all virtual observatories

CH_OPER_VOBS_4M_2_

CHAMP 4 month data from all virtual observatories

CR_OPER_VOBS_4M_2_

Cryosat-2 4 month data from all virtual observatories

OR_OPER_VOBS_4M_2_

Ørsted 4 month data from all virtual observatories

CO_OPER_VOBS_4M_2_

Composite (Ørsted, CHAMP, Cryosat-2, Swarm) 4 month data from all virtual observatories

These collections each contain the variables:

Variable

Unit

Dimension

Description

SiteCode

$\(-\)$

char [7]

virtual observatory identifier

Timestamp

$\(-\)$

scalar

UTC time of observation

Latitude

$\(\text{deg}\)$

scalar

ITRF geocentric latitude

Longitude

$\(\text{deg}\)$

scalar

ITRF geocentric longitude

Radius

$\(\text{m}\)$

scalar

ITRF geocentric radius

B_CF

$\(\text{nT}\)$

vector [3]

Core magnetic field vector in ITRF NEC frame.

B_OB

$\(\text{nT}\)$

vector [3]

Observed magnetic field vector in ITRF NEC frame.

sigma_CF

$\(\text{nT}\)$

vector [3]

Estimated error of the core magnetic field vector in ITRF NEC frame.

sigma_OB

$\(\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

SW_OPER_VOBS_1M_2_:SecularVariation

Swarm 1 month secular variation data from all virtual observatories

SW_OPER_VOBS_4M_2_:SecularVariation

Swarm 4 month secular variation data from all virtual observatories

CH_OPER_VOBS_4M_2_:SecularVariation

CHAMP 4 month secular variation data from all virtual observatories

CR_OPER_VOBS_4M_2_:SecularVariation

Cryosat-2 4 month secular variation data from all virtual observatories

OR_OPER_VOBS_4M_2_:SecularVariation

Ørsted 4 month secular variation data from all virtual observatories

CO_OPER_VOBS_4M_2_:SecularVariation

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

SiteCode

$\(-\)$

char [7]

virtual observatory identifier

Timestamp

$\(-\)$

scalar

UTC time of observation

Latitude

$\(\text{deg}\)$

scalar

ITRF geocentric latitude

Longitude

$\(\text{deg}\)$

scalar

ITRF geocentric longitude

Radius

$\(\text{m}\)$

scalar

ITRF geocentric radius

B_SV

$\(\text{nT}/\text{yr}\)$

vector [3]

Field secular variation vector in ITRF NEC frame.

sigma_SV

$\(\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:
    Spacecraft  (Timestamp) object '-' '-' '-' '-' '-' ... '-' '-' '-' '-' '-'
    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
    Latitude    (Timestamp) float64 65.26 65.26 65.26 ... 65.26 65.26 65.26
    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
    SiteCode    (Timestamp) <U7 'N65W051' 'N65W051' ... 'N65W051' 'N65W051'
    sigma_CF    (Timestamp, NEC) float64 5.858 2.708 5.326 ... 5.858 2.708 5.326
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));
../_images/79a7f11f37a7a43075511201514ec905d91d425f57a3ae6d0ceb33a955770793.png ../_images/7bdffd9dc7f1a14f18ed9ed5ae9c7a42c216933984548513270ae51d1be4b1ac.png

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:04<00:36,  4.05s/it]
 20%|██        | 2/10 [00:07<00:30,  3.85s/it]
 30%|███       | 3/10 [00:11<00:27,  3.87s/it]
 40%|████      | 4/10 [00:15<00:22,  3.78s/it]
 50%|█████     | 5/10 [00:19<00:19,  3.84s/it]
 60%|██████    | 6/10 [00:23<00:16,  4.09s/it]
 70%|███████   | 7/10 [00:27<00:11,  3.94s/it]
 80%|████████  | 8/10 [00:30<00:07,  3.76s/it]
 90%|█████████ | 9/10 [00:34<00:03,  3.63s/it]
100%|██████████| 10/10 [00:38<00:00,  3.70s/it]
100%|██████████| 10/10 [00:38<00:00,  3.80s/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: 27)
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.915e+03
    B_OB          (SiteCode, Timestamp, RTP) float64 5.284e+03 ... -6.913e+03
    sigma_CF      (SiteCode, Timestamp, RTP) float64 -1.379 -3.671 ... 7.107
    sigma_OB      (SiteCode, Timestamp, RTP) float64 -0.9393 -2.205 ... 11.54
    sigma_SV      (SiteCode, Timestamp_SV, RTP) float64 -1.195 -2.934 ... 7.586
    B_SV          (SiteCode, Timestamp_SV, RTP) float64 43.79 -10.41 ... -43.67
Attributes:
    Sources:         ['SW_OPER_VOBS_4M_2__20140301T000000_20231101T000000_0101']
    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");
../_images/09bc8c84064260700d727052c147c444a045adbc0bc529a9f116cf2c1eb4a197.png

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");
../_images/b942c37e266863ba3d54021e91ee3ad8d317f6e527fbec655685854d118f5d46.png

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")
../_images/1b7eda5d3b9a34c44dd083ad885a87478b5a750acdaba7d642f24a6ace5df095.png

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")
../_images/96cc59f7c9bdf9c7318670685baef5b3da6981244360409dfa27f52a88969b3e.png