import intake
import cartopy.crs as ccrs
import cartopy.feature as cf
import cmocean
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np

Plotting with cartopy

This function will ask the plot axes for their projection and use this information to sample (nearest-neighbor) variable in healpix nested ordering at a given resolution.

def nnshow(var, nx=1000, ny=1000, ax=None, **kwargs):
    """
    var: variable on healpix coordinates (array-like)
    nx: image resolution in x-direction
    ny: image resolution in y-direction
    ax: axis to plot on
    kwargs: additional arguments to imshow
    """
    if ax is None:
        ax = plt.gca()

    xlims = ax.get_xlim()
    ylims = ax.get_ylim()
    # NOTE: we want the center coordinate of each pixel, thus we have to
    # compute the linspace over halve a pixel size less than the plot's limits
    dx = (xlims[1] - xlims[0]) / nx
    dy = (ylims[1] - ylims[0]) / ny
    xvals = np.linspace(xlims[0] + dx / 2, xlims[1] - dx / 2, nx)
    yvals = np.linspace(ylims[0] + dy / 2, ylims[1] - dy / 2, ny)
    xvals2, yvals2 = np.meshgrid(xvals, yvals)
    latlon = ccrs.PlateCarree().transform_points(
        ax.projection, xvals2, yvals2, np.zeros_like(xvals2)
    )
    valid = np.all(np.isfinite(latlon), axis=-1)
    points = latlon[valid].T
    pix = hp.ang2pix(
        hp.npix2nside(len(var)), theta=points[0], phi=points[1], nest=True, lonlat=True
    )
    res = np.full(latlon.shape[:-1], np.nan, dtype=var.dtype)
    res[valid] = var[pix]
    return ax.imshow(res, extent=xlims + ylims, origin="lower", **kwargs)

of course, we also need some model output data:

cat = intake.open_catalog("https://data.nextgems-h2020.eu/catalog.yaml")
experiment = cat.ICON["ngc3028"]

For convenience, let’s also define a worldmap plot for later use:

def worldmap(var, **kwargs):
    projection = ccrs.Robinson(central_longitude=-135.5808361)
    fig, ax = plt.subplots(
        figsize=(8, 4), subplot_kw={"projection": projection}, constrained_layout=True
    )
    ax.set_global()

    nnshow(var, ax=ax, **kwargs)
    ax.add_feature(cf.COASTLINE, linewidth=0.8)
    ax.add_feature(cf.BORDERS, linewidth=0.4)

Let’s first have a look at global temperatures:

%%time
data = experiment(time="P1D", zoom=7).to_dask()
worldmap(data.tas.isel(time=0), cmap=cmocean.cm.thermal)

CPU times: user 970 ms, sys: 190 ms, total: 1.16 s
Wall time: 1.32 s

If we do the same for Europe, we might want to zoom in a little further:

%%time
data = experiment(time="PT3H", zoom=10).to_dask()

projection = ccrs.Robinson(central_longitude=10)
fig, ax = plt.subplots(
    figsize=(8, 4), subplot_kw={"projection": projection}, constrained_layout=True
)
ax.set_extent([-10, 30, 30, 60], crs=ccrs.PlateCarree())

nnshow(data.tas.isel(time=0), ax=ax, cmap=cmocean.cm.thermal)
ax.add_feature(cf.COASTLINE, linewidth=0.8)
ax.add_feature(cf.BORDERS, linewidth=0.4)

CPU times: user 358 ms, sys: 56.4 ms, total: 414 ms
Wall time: 408 ms

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fff87b6b5e0>

Just for reference, let’s also check zoom=7 from the world map, but applied to this region:

%%time
data = experiment(time="PT3H", zoom=7).to_dask()

projection = ccrs.Robinson(central_longitude=10)
fig, ax = plt.subplots(
    figsize=(8, 4), subplot_kw={"projection": projection}, constrained_layout=True
)
ax.set_extent([-10, 30, 30, 60], crs=ccrs.PlateCarree())

nnshow(data.tas.isel(time=0), ax=ax, cmap=cmocean.cm.thermal)
ax.add_feature(cf.COASTLINE, linewidth=0.8)
ax.add_feature(cf.BORDERS, linewidth=0.4)

CPU times: user 323 ms, sys: 27.4 ms, total: 350 ms
Wall time: 351 ms

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fff87883070>

We still get a good glimpse (might be enough for a quick check) but of course, we now can see the resolution limit.

Now how about the ice sheet on the north pole?

%%time
data = experiment(time="P1D", zoom=7).to_dask()

projection = ccrs.NorthPolarStereo()

fig, ax = plt.subplots(
    figsize=(8, 4), subplot_kw={"projection": projection}, constrained_layout=True
)
ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())

nnshow(data.conc.isel(time=0), ax=ax, cmap=cmocean.cm.ice)
ax.add_feature(cf.COASTLINE, linewidth=0.8)
ax.add_feature(cf.BORDERS, linewidth=0.4)

CPU times: user 750 ms, sys: 20.1 ms, total: 770 ms
Wall time: 773 ms

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fff85f56080>

And a few other variables…

%%time
worldmap(
    experiment(time="P1D", zoom=8).to_dask().u.isel(time=0, depth_full=0),
    cmap=cmocean.cm.balance,
    vmin=-1.5,
    vmax=1.5,
)

CPU times: user 366 ms, sys: 80.7 ms, total: 447 ms
Wall time: 401 ms

%%time
worldmap(
    experiment(time="P1D", zoom=8).to_dask().ua.isel(time=0, level_full=-1),
    cmap=cmocean.cm.balance,
    vmin=-15,
    vmax=15,
)

CPU times: user 367 ms, sys: 81.1 ms, total: 449 ms
Wall time: 405 ms

%%time
worldmap(
    experiment(time="P1D", zoom=8).to_dask().rsds.isel(time=0),
    cmap=cmocean.cm.thermal,
)

CPU times: user 312 ms, sys: 21 ms, total: 333 ms
Wall time: 335 ms