This notebook demonstrates how to use these tools together to rapidly build nice visualisations of geomagnetic models. A lot of things are used here without explanation (including xarray and cartopy), and you will encounter bugs and inconsistent behaviour with holoviews. Beware!
Probably refactor into separate notebooks. It is too long.
conda install -c pyviz holoviz
?conda install -c conda-forge holoviews geoviews hvplot ...
jupyter labextension install @pyviz/jupyterlab_pyviz
conda install xarray
pip install chaosmagpy
conda install selenium phantomjs
sudo apt install ffmpeg
assumed running on Linux
import datetime as dt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
import hvplot.xarray
import holoviews
import cartopy.crs as ccrs
import geoviews as gv
import geoviews.feature as gf
gv.extension('bokeh', 'matplotlib')
import eoxmagmod
from chaosmagpy.plot_utils import nio_colormap
def datetime_to_mjd2000(t: dt.datetime) -> float:
"Convert a datetime object to MJD2000."
# Convert to datetime64 ns
t = np.datetime64(t).astype("M8[ns]")
# Get offset to year 2000
t = (t - np.datetime64('2000')).astype('int64')
# Convert from ns to days
NS2DAYS = 1.0/(24*60*60*1e9)
return t * NS2DAYS
def generate_latlon_grid(resolution=2, min_lat=-90, max_lat=90):
"Generate a grid of positions over the Earth at a given degree resolution."
lat, lon = np.meshgrid(
np.arange(min_lat, max_lat, resolution),
np.arange(-180, 180, resolution))
REFRAD_KM = 6371.200
rad = np.ones_like(lat)*REFRAD_KM
return lat, lon, rad
def eval_model_on_grid(
lat: np.ndarray, lon: np.ndarray, rad: np.ndarray,
times: np.ndarray, model=None, shc_model=eoxmagmod.data.IGRF12,
**opts):
"""Use eoxmagmod to evaluate a model over a grid.
Evaluate the B_NEC vector at each point in a grid, for a range of times.
Args:
lat (ndarray): Latitude in degrees (Spherical geocentric)
lon (ndarray): Longitude in degrees
rad (ndarray): Geocentric radius in kilometres
times (ndarray): 1-D array of datetimes
model (magnetic_model): Model loaded with eoxmagmod.load_model_<x>
shc_model (str): Path to a shc model
Returns:
ndarray: B_NEC values at each point
"""
if shc_model and not model:
model = eoxmagmod.load_model_shc(shc_model)
times_mjd2000 = [datetime_to_mjd2000(t) for t in times]
# Reshape the input coordinates to use eoxmagmod
orig_shape = lat.shape
_lat, _lon, _rad = map(lambda x: x.flatten(), (lat, lon, rad))
coords = np.stack((_lat, _lon, _rad), axis=1)
coords = np.stack([coords for i in range(len(times_mjd2000))])
timestack = np.stack([np.ones_like(_lat)*t for t in times_mjd2000])
# Use the model and do the computation
b_nec = model.eval(timestack, coords, scale=[1, 1, -1], **opts)
# Reshape output back to original grid
b_nec = b_nec.reshape(times.shape + orig_shape + (3, ))
return b_nec
# The list of times to sample and the grid to use
times = np.array([dt.datetime(year, 1, 1) for year in range(1900, 2021, 10)])
lat, lon, rad = generate_latlon_grid(resolution=2)
# Evaluate the model
b_nec = eval_model_on_grid(lat, lon, rad, times, shc_model=eoxmagmod.data.IGRF13)
# b_nec = eval_model_on_grid(lat, lon, rad, times, shc_model="IGRF13.shc")
# Assign to an xarray.Dataset
ds = xr.Dataset({'B_NEC': (['time', 'x', 'y', 'NEC'], b_nec)},
coords={'lon': (['x', 'y'], lon),
'lat': (['x', 'y'], lat),
'time': times})
# Add some columns
# Intensity
ds["F"] = np.sqrt(np.sum(ds["B_NEC"]**2, axis=3))
# Declination, arctan(B_E / B_N)
ds["DEC"] = np.rad2deg(np.arctan(
ds["B_NEC"][:, :, :, 1] / ds["B_NEC"][:, :, :, 0]))
# # Alternative "expanded" form:
# ds = xr.Dataset(
# {'B_NEC_N': (['time', 'x', 'y'], b_nec[:, :, :, 0]),
# 'B_NEC_E': (['time', 'x', 'y'], b_nec[:, :, :, 1]),
# 'B_NEC_C': (['time', 'x', 'y'], b_nec[:, :, :, 2]),},
# coords={'lon': (['x', 'y'], lon),
# 'lat': (['x', 'y'], lat),
# 'time': times})
ds
ds.sel({"time": "1900"})["B_NEC"][:, :, :, 2].plot(x="lon", y="lat")
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mollweide())
ds.sel({"time": "1900"})["B_NEC"][:, :, :, 2].plot(x="lon", y="lat", ax=ax, transform=ccrs.PlateCarree())
ax.coastlines()
ax
ds["B_NEC"][:, :, :, 2].hvplot.contourf('lon', 'lat', levels=10)
ds["time"][-1].values
# Squeeze out just one time slice to simplify it
ds_sub = ds.sel({"time": "2020"}).squeeze("time")
# Subsample
ds_sub = ds_sub.isel({"x":slice(None, None, 2), "y":slice(None, None,2)})
# Some configuration to use for both plots
fig_kwargs = dict(
projection=ccrs.Mollweide(), global_extent=True, frame_height=400
)
title = """IGRF-13 evaluated at 01/01/2020"""
# Background colour: field intensity\n
# Red/blue contours: declination
# """
# Create a filled contour plot of intensity
ax1 = ds_sub.hvplot.contourf(
'lon', 'lat', 'F', levels=30, cmap='viridis', alpha=0.8,
coastline=True, title=title, hover=False, clabel="Intensity / nT",
**fig_kwargs
)
# .. and overlay with declination contours
ax2 = ds_sub.hvplot.contour(
'lon', 'lat', 'DEC', levels=30, cmap="seismic", line_width=2, **fig_kwargs
)
ax3 = ax1 * ax2
# ax3.opts(tools="box_zoom,reset,tap".split(","), toolbar="right")
print(ax3)
ax3
hvplot.save(ax3, "outputs/hvplot_igrf13_2020.html")
# Some configuration to use for both plots
# NB must set dynamic=False so that the bokeh plot does not rely on a server
fig_kwargs = dict(
projection=ccrs.Mollweide(), global_extent=True, frame_height=300, dynamic=False,
)
# Create a filled contour plot of intensity
# NB need to set the contours explicitly in order to fix them across timeslices
contours = np.arange(30000, 70000, 5000)
ax1 = ds.hvplot.contourf(
'lon', 'lat', 'F', levels=contours, cmap='viridis', alpha=0.8,
coastline=True, hover=False, clabel="Intensity / nT", clim=(30000, 70000),
**fig_kwargs
)
# ax2 = ds.hvplot.contour(
# 'lon', 'lat', 'DEC', levels=40, cmap="seismic", **fig_kwargs
# )
# ax1 * ax2
ax1
ds_sub = ds.isel(x=slice(0, -1, 2), y=slice(0, -1, 2), time=slice(0, -1, 2))
# Some configuration to use for both plots
# NB must set dynamic=False so that the bokeh plot does not rely on a server
fig_kwargs = dict(
projection=ccrs.Mollweide(), global_extent=True, frame_height=300,
dynamic=False
)
# Create a filled contour plot of intensity
ax1 = ds_sub.hvplot.quadmesh(
'lon', 'lat', "F", cmap='viridis', alpha=0.8,
coastline=True, hover=False, clabel="Intensity / nT", clim=(30000, 70000),
rasterize=True, datashade=True, aggregator='mean',
**fig_kwargs
)
ax2 = ds_sub.hvplot.contour(
'lon', 'lat', 'DEC', levels=40, cmap="seismic", colorbar=False, legend=False,
**fig_kwargs
)
ax1 * ax2
# hvplot.save(ax1*ax2, "hvplot_igrf_timeslider.html")
# see bug: https://github.com/holoviz/hvplot/issues/305
renderer = holoviews.renderer('bokeh')
renderer.save(ax1*ax2, 'hvplot_igrf_timeslider')
!du -h hvplot_igrf_timeslider.html
renderer = holoviews.renderer('bokeh')
renderer.save(ax1*ax2, 'hvplot_igrf_timeslider_b', fmt="scrubber")
holoviews.save(ax, 'hvplot_igrf_timeslider_b.html', backend='bokeh', fmt="scrubber")
# The list of times to sample and the grid to use
times = np.array([dt.datetime(2000, 1, 1)])
lat, lon, rad = generate_latlon_grid(resolution=0.4)
# Evaluate the model
b_nec = eval_model_on_grid(lat, lon, rad, times, model=eoxmagmod.data.LCS1)
# Assign to an xarray.Dataset
ds_lcs = xr.Dataset(
{'B_NEC_N': (['time', 'x', 'y'], b_nec[:, :, :, 0]),
'B_NEC_E': (['time', 'x', 'y'], b_nec[:, :, :, 1]),
'B_NEC_C': (['time', 'x', 'y'], b_nec[:, :, :, 2]),},
coords={'lon': (['x', 'y'], lon),
'lat': (['x', 'y'], lat),
'time': times}).squeeze("time")
ds_lcs
fig_kwargs = dict(
projection=ccrs.Mollweide(), global_extent=True, frame_height=300,
rasterize=True, project=True, dynamic=False
#datashade=True, #aggregator="mean"#
)
ax1 = ds.hvplot.quadmesh(
x='lon', y='lat', z="B_NEC_C", cmap=nio_colormap(), alpha=0.4,
# hover=True, hover_cols=["B_NEC_N", "B_NEC_E", "B_NEC_C"],
hover_cols=["B_NEC_C", "lat", "lon"],
coastline=True, clabel="Vertical component / nT",
clim=(-200, 200), colorbar=True, legend=True,
**fig_kwargs
)
ax1
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mollweide())
ds["B_NEC_C"].plot(
x="lon", y="lat", ax=ax, transform=ccrs.PlateCarree(),
cmap=nio_colormap(), vmin=-200, vmax=200)
ax.coastlines()
ax.set_title("LCS-1 vertical component / nT", fontdict={"fontsize": 25})
ax
import urllib.request as request
import zipfile
from tempfile import NamedTemporaryFile
import shutil
def fetch_zipped_file(url, file_name):
"Fetch a given file from an online zip file"
output_file = NamedTemporaryFile()
zip_file, _ = request.urlretrieve(url)
with zipfile.ZipFile(zip_file, 'r') as zip_ref:
with zip_ref.open(file_name) as f:
shutil.copyfileobj(f, output_file)
output_file.seek(0)
return output_file
def load_mio():
url = 'ftp://swarm-diss.eo.esa.int/Level2longterm/MIO/SW_OPER_MIO_SHA_2C_20131201T000000_20170101T000000_0201.ZIP'
file_name = 'SW_OPER_MIO_SHA_2C_20131201T000000_20170101T000000_0201.txt'
# url = 'ftp://swarm-diss.eo.esa.int/Level2longterm/MIO/SW_OPER_MIO_SHA_2D_20131201T000000_20171231T235959_0402.ZIP'
# file_name = 'SW_OPER_MIO_SHA_2D_20131201T000000_20171231T235959_0402.txt'
mio_file = fetch_zipped_file(url, file_name)
mio_model = eoxmagmod.load_model_swarm_mio_external(mio_file.name)
return mio_model
mio_model = load_mio()
def eval_mio_model(times=None, **kwargs):
lat, lon, rad = generate_latlon_grid(**kwargs)
# Evaluate the model
# NB need to pass a value for F107 for the MIO models
b_nec = eval_model_on_grid(lat, lon, rad, times, model=mio_model, f107=70)
# Assign to an xarray.Dataset
ds = xr.Dataset({'B_NEC': (['time', 'x', 'y', 'NEC'], b_nec)},
coords={'lon': (['x', 'y'], lon),
'lat': (['x', 'y'], lat),
'time': times})
return ds
# Create the diurnal variation by sampling hours through one day
# t0 = dt.datetime(2017, 1, 1)
# times = np.array([t0 + dt.timedelta(hours=i) for i in np.linspace(0, 25, 6)])
times = np.array([dt.datetime(2017, 1, 1, hour) for hour in range(0, 24, 4)])
ds_day = eval_mio_model(times)
# Create the seasonal variation by sampling noon-times through one year
# times = np.array([dt.datetime(2017, int(i), 1, 12) for i in np.linspace(1, 12, 6)])
times = np.array([dt.datetime(2017, month, 1, 12) for month in range(1, 12, 2)])
ds_year = eval_mio_model(times)
# ds_day_sub = ds_day.isel(x=slice(0, -1, 2), y=slice(0, -1, 2))#, time=slice(0, -1, 2))
# ds_day_sub = ds_day.isel(time=slice(0, -1, 2))
fig_kwargs = dict(
projection=ccrs.Mollweide(), global_extent=True, frame_height=300,
coastline=True, x='lon', y='lat',
cmap=nio_colormap(), alpha=0.4, clim=(-50, 50),# levels=10,
clabel="nT", colorbar=True, legend=True,
rasterize=True, project=True, dynamic=False, hover=False,
# datashade=True, aggregator="mean"#
)
# ax1 = ds_day["B_NEC"][:, :, :, 2].hvplot.quadmesh(**fig_kwargs).options(toolbar=None)
# ax2 = ds_year["B_NEC"][:, :, :, 2].hvplot.quadmesh(**fig_kwargs)
# ax1
times = np.array([dt.datetime(2017, 1, 1, hour) for hour in range(0, 24, 1)])
ds_day = eval_mio_model(times, resolution=2, min_lat=-60, max_lat=60)
times = np.sort(np.array([dt.datetime(2017, month, 1, 12) for month in range(1, 13, 1)]
+ [dt.datetime(2017, month, 15, 12) for month in range(1, 13, 1)]))
ds_year = eval_mio_model(times, resolution=2, min_lat=-60, max_lat=60)
# # # this crashes the kernel:
# holoviews.renderer('matplotlib').save(ax1, 'test', fmt='gif')
ds_day.isel(time=[0]).squeeze("time")["B_NEC"][:, :, 2].hvplot.contourf(levels=10, **fig_kwargs,).options(toolbar=None)
## this was behaving weirdly before.. (ConnectioRefusedError...)
## seems to work now when running notebook without any previous holoviews visible
for i in range(ds_day.time.size):
# Get a time string to add to the plots...
t = pd.to_datetime(ds_day.time.data[i]).strftime('%Y-%m-%d %H:%M:%S')
ax = ds_day.isel(time=[i]).squeeze("time")["B_NEC"][:, :, 2].hvplot.quadmesh(**fig_kwargs, title=t).options(toolbar=None)
hvplot.save(ax, f'mio_day_{i:02}.png')
# holoviews.save(ax, f'mio_day_{i:02}.png')
# # using imagemagick
# # convert doesn't have webm?
# # gifs are much larger
# !convert -delay 25 -loop 0 mio_day_*.png mio_day.gif
!ffmpeg -framerate 8 -i mio_day_%2d.png mio_day.webm -y
# !rm mio_day_*
for i in range(ds_year.time.size):
t = pd.to_datetime(ds_year.time.data[i]).strftime('%Y-%m-%d %H:%M:%S')
ax = ds_year.isel(time=[i]).squeeze("time")["B_NEC"][:, :, 2].hvplot.quadmesh(**fig_kwargs, title=t).options(toolbar=None)
hvplot.save(ax, f'mio_year_{i:02}.png')
!ffmpeg -framerate 8 -i mio_year_%2d.png mio_year.webm -y
!rm mio_day_*
!rm mio_year_*
def load_mma():
url = 'ftp://swarm-diss.eo.esa.int/Level2longterm/MMA/SW_OPER_MMA_SHA_2C_20131201T000000_20180101T000000_0401.ZIP'
file_name = 'SW_OPER_MMA_SHA_2C_20131201T000000_20180101T000000_0401.cdf'
mma_file = fetch_zipped_file(url, file_name)
mma_model = eoxmagmod.load_model_swarm_mma_2c_external(mma_file.name)
return mma_model
mma_model = load_mma()
mma_model
def eval_mma_model(times=None, **kwargs):
lat, lon, rad = generate_latlon_grid(**kwargs)
# Evaluate the model
# NB need to pass a value for F107 for the MIO models
b_nec = eval_model_on_grid(lat, lon, rad, times, model=mma_model)
# Assign to an xarray.Dataset
ds = xr.Dataset({'B_NEC': (['time', 'x', 'y', 'NEC'], b_nec)},
coords={'lon': (['x', 'y'], lon),
'lat': (['x', 'y'], lat),
'time': times})
return ds
# times = np.array([dt.datetime(2017, 1, 1, hour) for hour in range(0, 24, 4)])
# ds_day = eval_mma_model(times, resolution=4)
# ds_day["B_NEC"][:, :, :, 2].hvplot.quadmesh(**fig_kwargs).options(toolbar=None)
t0 = dt.datetime(2017, 1, 1)
times = np.array([t0 + dt.timedelta(hours=i) for i in range(0, 24*3, 1)])
ds_day = eval_mma_model(times, resolution=2)
## this was behaving weirdly before.. (ConnectioRefusedError...)
## seems to work now when running notebook without any previous holoviews visible
for i in range(ds_day.time.size):
# Get a time string to add to the plots...
t = pd.to_datetime(ds_day.time.data[i]).strftime('%Y-%m-%d %H:%M:%S')
ax = ds_day.isel(time=[i]).squeeze("time")["B_NEC"][:, :, 2].hvplot.quadmesh(**fig_kwargs, title=t).options(toolbar=None)
hvplot.save(ax, f'mma_day_{i:02}.png')
# holoviews.save(ax, f'mio_day_{i:02}.png')
!ffmpeg -framerate 8 -i mma_day_%2d.png mma_days.webm -y
!rm mma_day_*