Abstract: We Demonstrate geomagnetic ground observatory data access through VirES - this is the AUX_OBS product distributed by BGS to support the Swarm mission, and contains data from INTERMAGNET and the World Data Centre (WDC) for Geomagnetism. Data are available as three collections: 1 second and 1 minute cadences (INTERMAGNET definitive & quasi-definitive data), as well as specially derived hourly means over the past century (WDC).
See also:
The data are also available from the BGS FTP server (ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS/
). If that is more useful to you, you can refer to this older notebook demonstration of access to the FTP server.
Please note the data are under different usage terms than the Swarm data:
Caution: The magnetic vector components have been rotated into the geocentric (NEC) frame rather than the geodetic frame, so that they are consistent with the Swarm data. This is in contrast with the data provided directly from observatories.
# Display important package versions used
%load_ext watermark
%watermark -i -v -p viresclient,pandas,xarray,matplotlib
from viresclient import SwarmRequest
import matplotlib.pyplot as plt
Data are organised into AUX_OBSH (hour), AUX_OBSM (minute), AUX_OBSS (second) types. For example, to access the hourly data, use the collection name SW_OPER_AUX_OBSH2_
.
request = SwarmRequest()
print(request.available_collections("AUX_OBSH", details=False))
print(request.available_collections("AUX_OBSM", details=False))
print(request.available_collections("AUX_OBSS", details=False))
Within each collection, the following variables are available:
print(request.available_measurements("SW_OPER_AUX_OBSH2_"))
print(request.available_measurements("SW_OPER_AUX_OBSM2_"))
print(request.available_measurements("SW_OPER_AUX_OBSS2_"))
B_NEC
and F
are the magnetic field vector and intensityIAGA_code
gives the official three-letter IAGA codes that identify each observatoryQuality
is either "D" or "Q" to indicate whether data is definitive (D) or quasi-definitive (Q)ObsIndex
is an increasing integer (0, 1, 2...) attached to the hourly data - this indicates a change in the observatory (e.g. of precise location) while the 3-letter IAGA code remained the sameNote that the IAGA_code
variable is necessary in order to distinguish records from each observatory
Note that there is a special message issued regarding the data terms
Let's fetch all the variables available within the 1-minute data, from two days:
request = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSM2_")
request.set_products(["IAGA_code", "B_NEC", "F", "Quality"])
data = request.get_between("2016-01-01", "2016-01-03")
ds = data.as_xarray()
ds
Above, we loaded the data as an xarray Dataset
, but we could also load the data as a pandas DataFrame
- note that we should use expand=True
to separate the vector components of B_NEC
into distinct columns:
df = data.as_dataframe(expand=True)
df
available_observatories
to find possible IAGA codes¶We can get a dataframe containing the availability times of data from all the available observatories for a given collection:
request.available_observatories("SW_OPER_AUX_OBSM2_", details=True)
We can also get a list of only the available observatories during a given time window:
print(request.available_observatories("SW_OPER_AUX_OBSM2_", '2016-01-01', '2016-01-02'))
IAGA_code
to specify a particular observatory¶Subset the collection with a special collection name like "SW_OPER_AUX_OBSM2_:<IAGA_code>"
to get data from only that observatory:
request = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSM2_:ABK")
request.set_products(["IAGA_code", "B_NEC", "F", "Quality"])
data = request.get_between("2016-01-01", "2016-01-03")
ds = data.as_xarray()
ds
The VirES API treats these data similarly to the Swarm MAG products, and so all the same model handling behaviour applies. For example, we can directly remove the CHAOS core and crustal model predictions:
request = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSM2_:ABK")
request.set_products(
measurements=["B_NEC"],
models=["'CHAOS-internal' = 'CHAOS-Core' + 'CHAOS-Static'"],
residuals=True
)
data = request.get_between("2016-01-01", "2016-01-03")
ds = data.as_xarray()
ds["B_NEC_res_CHAOS-internal"].plot.line(x="Timestamp", col="NEC")
(This roughly shows the disturbance sensed by the observatory due to the magnetospheric and ionospheric sources)
Let's run through a visualisation of one year of hourly means from three observatories.
First fetch the data from our chosen observatories across the UK: LER (Lerwick), ESK (Eskdalemuir), HAD (Hartland). We can apply a few optional settings to reduce unnecessary output:
verbose=False
to disable the data terms messageasynchronous=False
to enable synchronous processing on the server - it will be slightly faster but only works for smaller data requestsshow_progress=False
to hide the progress barsSpacecraft
variable when we load the datarequest = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSH2_:LER", verbose=False)
request.set_products(measurements=["B_NEC"])
data = request.get_between("2013-01-01", "2014-01-01", asynchronous=False, show_progress=False)
ds_ler = data.as_xarray().drop("Spacecraft")
request = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSH2_:ESK", verbose=False)
request.set_products(measurements=["B_NEC"])
data = request.get_between("2013-01-01", "2014-01-01", asynchronous=False, show_progress=False)
ds_esk = data.as_xarray().drop("Spacecraft")
request = SwarmRequest()
request.set_collection("SW_OPER_AUX_OBSH2_:HAD", verbose=False)
request.set_products(measurements=["B_NEC"])
data = request.get_between("2013-01-01", "2014-01-01", asynchronous=False, show_progress=False)
ds_had = data.as_xarray().drop("Spacecraft")
Now our data is in three objects which look like this:
ds_ler
We can quickly preview the data using the xarray plotting tools:
ds_ler["B_NEC"].plot.line(x="Timestamp", col="NEC", sharey=False)
Let's make a more complex figure to display data from all three observatories together. We can use matplotlib directly now to create the figure and pass the xarray objects to it to fill the contents. Note that we slice out a particular vector component with e.g. ds_ler["B_NEC"].sel(NEC="N")
.
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(15, 5), sharex="all", sharey="row")
for i, NEC in enumerate("NEC"):
axes[i, 0].plot(ds_ler["Timestamp"], ds_ler["B_NEC"].sel(NEC=NEC))
axes[i, 1].plot(ds_esk["Timestamp"], ds_esk["B_NEC"].sel(NEC=NEC))
axes[i, 2].plot(ds_had["Timestamp"], ds_had["B_NEC"].sel(NEC=NEC))
axes[i, 0].set_ylabel(f"B ({NEC}) [nT]")
axes[0, 0].set_title("LER: Lerwick (60.0°N)")
axes[0, 1].set_title("ESK: Eskdalemuir (55.1°N)")
axes[0, 2].set_title("HAD: Hartland (50.8°N)")
fig.tight_layout()
This shows us the difference in the main field between these locations - further North (Lerwick), the field is pointing more downwards so the vertical component (C) is stronger. We can also see a small secular variation over the year as the field changes.