%load_ext autoreload
%autoreload 2
import subprocess
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 as hv
import cartopy.crs as ccrs
import geoviews as gv
import geoviews.feature as gf
gv.extension('bokeh', 'matplotlib')
from bokeh.models import PrintfTickFormatter
import eoxmagmod
from chaosmagpy.plot_utils import nio_colormap
from src import build_model_xarray
# hv.help(hv.Contours)
# TODO: find gufm1 and stitch it with IGRF13
times = np.array([dt.datetime(year, 1, 1) for year in range(1900, 2026, 5)])
ds = build_model_xarray("IGRF13.shc", times)
# Some configuration to use for all 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=200,
hover=False, dynamic=False
)
def quadmesh(ds, component, **opts):
"""Plot X, Y, Z, F, H with nT"""
_fig_kwargs = fig_kwargs.copy()
if component in ["B_N", "B_E", "B_C"]:
_fig_kwargs.update(dict(
cmap = nio_colormap(), clim=(-70000, 70000)))
elif component in ["F", "H"]:
_fig_kwargs.update(dict(
cmap = "viridis", clim=(20000, 70000)))
_fig_kwargs.update(dict(
alpha=0.8, clabel="nT", coastline=True,
rasterize=True))#, datashade=True, aggregator="mean",))
_fig_kwargs.update(opts)
return ds.hvplot.quadmesh(
'lon', 'lat', component,
**_fig_kwargs
)
def contourf(ds, component, unit="nT", **opts):
"""Currently crashes on DEC/INC ..."""
ds = ds.copy()
_fig_kwargs = fig_kwargs.copy()
if component in ["B_N", "B_E", "B_C"]:
_fig_kwargs.update(dict(
clabel="nT",
cmap=nio_colormap(), levels=np.linspace(-70000, 70000, 51)))
elif component in ["F", "H"]:
_fig_kwargs.update(dict(
clabel="nT",
cmap="viridis", levels=np.linspace(20000, 70000, 51)))
elif component in ["DEC", "INC"]:
fig_kwargs.update(dict(
clabel="deg",
cmap=nio_colormap(), levels=np.linspace(-90, 90, 37)))
# NB need to set the contours explicitly in order to fix them across timeslices
_fig_kwargs.update(dict(
alpha=0.8, coastline=True))
if unit == "nT":
pass
elif unit == "deg":
pass
elif unit == "μT":
ds[component] = _ds[component]/1e3
_fig_kwargs["levels"] = _fig_kwargs["levels"]/1e3
_fig_kwargs["clabel"] ="μT"
_fig_kwargs.update(opts)
return ds.hvplot.contourf(
'lon', 'lat', component,
**_fig_kwargs
)
def contour(ds, component, color_continents=False, **opts):
"""Plot DEC or INC in degrees"""
_fig_kwargs = fig_kwargs.copy()
_fig_kwargs.update(dict(
levels=30, clim=(-90, 90), cmap=nio_colormap(), line_width=2, hover=True,
coastline=True,
# colorbar=False, legend=False
))
_fig_kwargs.update(opts)
ax = (
ds.hvplot.contour(
'lon', 'lat', component,
**_fig_kwargs)
# set degree formatter on colorbar
.opts({"Contours": {"colorbar_opts": {"formatter": PrintfTickFormatter(format="%i°")}}}))
if color_continents:
ax = (gf.land(fill_color="darkseagreen", alpha=0.6) *
gf.ocean(fill_color="steelblue", alpha=0.6) *
ax)
return ax
# _ds = ds.isel(time=[0]).squeeze("time")
# contourf(_ds, "DEC")
# # identify integer contour options
# for i in range(55):
# a = np.linspace(20000, 70000, i)
# if all([(num%1 == 0) for num in a]):
# print(i, a)
_ds = ds.isel(x=slice(0, -1, 4), y=slice(0, -1, 4), time=slice(0, -1, 4))
ax1 = quadmesh(_ds, "B_C", alpha=0.5)
# NB should be more accurate to use datashader, but seems to glitch and fill in the background of the figure
# hv.save(ax1, "test1.html")
# ax1
_ds = ds.isel(time=slice(0, -1, 4))
ax2 = contourf(_ds, "F", alpha=0.6)
ax2
# hv.save(ax2, "test4.html")
# hv.save(ax2, "test6.html", fmt="scrubber")
# contourf plots are ~10x smaller in filesize than quadmesh !!!
# BUT GLITCHES AT SOME TIME STEPS. try Martin's suggestion of projecting *before* plotting
# NB there are also glitches at the antemeridian but not really noticeable when using high sampling
# https://github.com/pacesm/jupyter_notebooks/blob/master/examples/CHAOS-6_Cartopy_Contours.ipynb
# # how to add labels to the contours?
# hv.help(hv.Contours)
_ds = ds.isel(time=[-2])
# ax3 = contour(_ds, "DEC", hover=True, colorbar=False, legend=False, color_continents=True)
# ax3
ax = (
contourf(_ds, "F", levels=np.linspace(20, 70, 11), unit="μT", clabel="Intensity [μT]",
frame_height=300, title="IGRF-13 evaluated at 01/01/2020")
.opts({"Polygons": {"colorbar_position": "right"}})#, "colorbar_opts":{"title_standoff": 10}}})#, "title_text_font_size":12}}})
*
contour(_ds, "DEC", hover=True, color_continents=False, clabel="Declination [deg]",
colorbar=True, legend=True,
)
.opts({"Contours": {"colorbar_position": "left"}})
)
# hv.save(ax, "hvplot_igrf13_2020b.html")
hvplot.save(ax.opts(toolbar=None), "igrf13_2020_FD.png")
ax
# Increase spacing of colorbar tick labels from the bar
formatter = PrintfTickFormatter(format=" %i")
for i in range(ds.time.size):
# Get a time string to add to the plots...
t = pd.to_datetime(ds.time.data[i]).strftime('%Y-%m-%d %H:%M:%S')
_ds = ds.isel(time=[i]).squeeze("time")
ax = (
quadmesh(_ds, "B_C", colorbar=True, frame_height=300)
.opts(toolbar=None, title=t)
.opts({"Image":
{"colorbar_opts": {"formatter": formatter}}})
)
hvplot.save(ax, f'tmp/frame_{i:03}.png')
subprocess.run("ffmpeg -framerate 8 -i tmp/frame_%3d.png core_igrf_Z.webm -y".split(" "), check=True)
for i in range(ds.time.size):
# Get a time string to add to the plots...
# t = pd.to_datetime(ds.time.data[i]).strftime('%Y-%m-%d %H:%M:%S') + " (IGRF-13)"
t = pd.to_datetime(ds.time.data[i]).strftime('%Y-%m-%d') + " (IGRF-13)"
_ds = ds.isel(time=[i]).squeeze("time")
# Convert to uT
for var in ["B_N", "B_E", "B_C", "F"]:
_ds[var] = _ds[var]/1e3
ax = (
quadmesh(_ds, "B_N", title="Northward (X)", clabel="μT", clim=(-70, 70)) +
quadmesh(_ds, "B_E", title="Eastward (Y)", clabel="μT", clim=(-70, 70)) +
quadmesh(_ds, "B_C", title="Downward (Z)", clabel="μT", clim=(-70, 70)) +
contour(_ds, "DEC", title="Declination (D)", color_continents=True) +
contour(_ds, "INC", title="Inclination (I)", color_continents=True) +
quadmesh(_ds, "F", title="Intensity (F)", clabel="μT", clim=(20, 70))
).options(toolbar=None, title=t).cols(3)
hvplot.save(ax, f'tmp/frame_{i:03}.png')
# remove unintended background transparency (requires imagemagick)
subprocess.run("mogrify -flatten tmp/*.png".split(" "))
subprocess.run("ffmpeg -framerate 8 -i tmp/frame_%3d.png core_igrf_fullvector.webm -y".split(" "), check=True)
# subprocess.run("rm tmp/frame_*".split(" "))
# !rm tmp/frame_*
for i in range(ds.time.size):
t = pd.to_datetime(ds.time.data[i]).strftime('%Y-%m-%d') + " (IGRF-13)"
_ds = ds.isel(time=[i]).squeeze("time")
ax = (
contourf(_ds, "B_N", title="Northward (X)", unit="μT") +
contourf(_ds, "B_E", title="Eastward (Y)", unit="μT") +
contourf(_ds, "B_C", title="Downward (Z)", unit="μT") +
contour(_ds, "DEC", title="Declination (D)") +
contour(_ds, "INC", title="Inclination (I)") +
contourf(_ds, "F", title="Intensity (F)", unit="μT")
).options(toolbar=None, title=t).cols(3)
hvplot.save(ax, f'tmp/frame_{i:03}.png')
# remove unintended background transparency (requires imagemagick)
subprocess.run("mogrify -flatten tmp/*.png".split(" "))
subprocess.run("ffmpeg -framerate 8 -i tmp/frame_%3d.png core_igrf_fullvector_2.webm -y".split(" "), check=True)