[1]:

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#

The easygems package provides some convenience functions for plotting with matplotlib, so let’s use it as well. You might remember seeing nnshow in this place: egh.healpix_show is it’s successor.

[2]:

import easygems.healpix as egh

of course, we also need some model output data:

[3]:

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:

[4]:

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

    egh.healpix_show(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:

[5]:

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

CPU times: user 1.08 s, sys: 260 ms, total: 1.34 s
Wall time: 1.66 s

../../_images/Processing_healpix_healpix_cartopy_8_1.png

regional plots#

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

[6]:

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

egh.healpix_show(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 191 ms, sys: 73.6 ms, total: 265 ms
Wall time: 249 ms

[6]:

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

../../_images/Processing_healpix_healpix_cartopy_10_2.png

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

[7]:

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

egh.healpix_show(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 139 ms, sys: 18.3 ms, total: 157 ms
Wall time: 159 ms

[7]:

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

../../_images/Processing_healpix_healpix_cartopy_12_2.png

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

If you like to work with lower resolution data (or don’t have anything else), but still like smooth plots, you can use the interpolating variant of healpix_show by passing the method="linear" argument.

[8]:

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

egh.healpix_show(data.tas.isel(time=0), ax=ax, cmap=cmocean.cm.thermal, method="linear")
ax.add_feature(cf.COASTLINE, linewidth=0.8)
ax.add_feature(cf.BORDERS, linewidth=0.4)

CPU times: user 158 ms, sys: 16.2 ms, total: 174 ms
Wall time: 181 ms

[8]:

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

../../_images/Processing_healpix_healpix_cartopy_14_2.png

On top, the easygems package also provides a method to show contour lines directly from healpix data:

[9]:

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

projection = ccrs.PlateCarree()

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

egh.healpix_show(
    data.prw.isel(time=0), ax=ax, cmap="Blues", vmin=45, vmax=75, method="linear"
)
contour_lines = egh.healpix_contour(
    data.prw.isel(time=0), ax=ax, levels=[45, 50], colors="grey", linewidths=1, alpha=1
)
plt.clabel(contour_lines, inline=True, fontsize=8, colors="grey", fmt="%d")


ax.add_feature(cf.COASTLINE, linewidth=0.8)
ax.add_feature(cf.BORDERS, linewidth=0.4)

CPU times: user 181 ms, sys: 32.8 ms, total: 214 ms
Wall time: 216 ms

../../_images/Processing_healpix_healpix_cartopy_16_1.png

Now how about the ice sheet on the north pole?

[10]:

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

map_extent = [-65, -5, -10, 25]
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())

egh.healpix_show(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)

[10]:

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

../../_images/Processing_healpix_healpix_cartopy_18_1.png

And a few other variables…

[11]:

%%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 214 ms, sys: 124 ms, total: 338 ms
Wall time: 259 ms

../../_images/Processing_healpix_healpix_cartopy_20_1.png

[12]:

%%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 182 ms, sys: 106 ms, total: 288 ms
Wall time: 235 ms

../../_images/Processing_healpix_healpix_cartopy_21_1.png

[13]:

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

CPU times: user 131 ms, sys: 9.04 ms, total: 140 ms
Wall time: 145 ms

../../_images/Processing_healpix_healpix_cartopy_22_1.png