:

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
"""
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 - xlims) / nx
dy = (ylims - ylims) / ny
xvals = np.linspace(xlims + dx / 2, xlims - dx / 2, nx)
yvals = np.linspace(ylims + dy / 2, ylims - 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, phi=points, 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)


Let’s first have a look at global temperatures:

:

%%time
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

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)

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

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)

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

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)

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

CPU times: user 312 ms, sys: 21 ms, total: 333 ms 