Add cartopy Greenland map to plotted data

[yoniv1]

Notebook description

Plotting Greenland ice flow velocity

Notebook link or upload

https://github.com/bopen/c3s-eqc-toolbox-template/blob/main/notebooks/wp5/greenland_ice_velocity.ipynb

Anything else we need to know?

Hi Mattia,

Is it possible to include a cartopy background in the plotting function of this notebook?

Right now, it is:

def plot_maps(da, suptitle=None, **kwargs):
    kwargs.setdefault("cmap", "rainbow")
    kwargs.setdefault("col", "period" if "period" in da.dims else None)
    if kwargs["col"]:
        kwargs.setdefault("subplot_kws", {})
        kwargs["subplot_kws"].setdefault("aspect", "equal")
    plot_obj = da.plot.imshow(**kwargs)
    if kwargs["col"]:
        for ax in plot_obj.axs.flatten():
            ax.axis("off")
    else:
        plt.axis("equal")
        plt.axis("off")
    return plot_obj

da = mean_std
_ = plot_maps(
    da,
    norm=matplotlib.colors.LogNorm(),
    vmin=1e-3,
    vmax=1e-1,
    suptitle="Horizontal velocity standard deviation of the Greenland Ice Sheet.",
)

which results in:

is it possible to integrate the cartopy background (also note that the x,y coordinates in the dataset are given in EPSG 3413):

projection = ccrs.Stereographic(central_latitude=70, central_longitude=-42)

fig, ax = plt.subplots(
    figsize=(6, 12), facecolor="w",
    subplot_kw=dict(projection=projection),
)

ax.coastlines()
ax.add_feature(cfeature.LAND,color='w')
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, linewidth=0.25, alpha=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
ax.set_extent((-65, -15, 58, 84), crs=ccrs.PlateCarree())

Environment

[yoniv1]

It can also be another basemap, if it is easier to keep working in the projected coordinates. I dont know if that exists though. But I would anyway like some kind of background with an indication of the land borders and the ocean as a background. And the coordinates. Is it possible to plot each 4 years (2017-18,2018-19,2019-20,2020-21) in a 2x2 subplot? The year 2020-21 is available from version 1.4 on the CDS.

[yoniv1]

Hi Mattia,

I tried something on my own:

def plot_maps_and_print_stats(da, suptitle=None, **kwargs):
    print(f"Mean {da.attrs['long_name']} [{da.attrs['units']}]")
    print(da.mean(("x", "y")).to_pandas())

    kwargs.setdefault("cmap", "rainbow")

    # Create 2x2 subplots
    fig, axs = plt.subplots(2, 2, figsize=(12, 10), subplot_kw={'projection': ccrs.Stereographic()})
    axs = axs.flatten()  # Flatten the array of axes for easy iteration

    # Plot each period on a different subplot
    periods = da.coords['period'].values if 'period' in da.coords else [None]

    for i, (ax, period) in enumerate(zip(axs, periods)):
        if period is not None:
            subset_da = da.sel(period=period)
        else:
            subset_da = da
        im=subset_da.plot.imshow(ax=ax, **kwargs, add_colorbar=False,transform=ccrs.epsg(3413))
        ax.set_title(f'Period: {period}' if period is not None else 'Full Data')
        ax.axis("off")
        ax.add_feature(cfeature.LAND, edgecolor='black')
        ax.add_feature(cfeature.OCEAN, color='lightblue')
        ax.coastlines()
        ax.gridlines(draw_labels=True)
        ax.set_extent([-638000, 860000, -3450000, -656000], ccrs.Stereographic(
        central_longitude=-45,
        central_latitude=90,
        true_scale_latitude=70))
        fig.colorbar(im,ax=ax,extend='both',shrink=0.49,label=f"{da.attrs['long_name']} [{da.attrs['units']}]")

    plt.tight_layout()  
    plt.show()

But the resolution of the plots lowered. I don't know what might have caused that...

[malmans2]

Hi @yoniv1,

I've never worked with this ESPG images and I'm not a Cartopy expert, so I'm not sure if this is a bug or the expected behaviour. The resampling is done when you set the extent, so if you move the imshow command at the bottom you don't get this issue.

This works fine for me:

def plot_maps_and_print_stats(da, suptitle=None, **kwargs):
    print(f"Mean {da.attrs['long_name']} [{da.attrs['units']}]")
    print(da.mean(("x", "y")).to_pandas())

    kwargs.setdefault("cmap", "rainbow")

    # Create 2x2 subplots
    fig, axs = plt.subplots(
        2, 2, figsize=(12, 10), subplot_kw={"projection": ccrs.Stereographic()}
    )
    axs = axs.flatten()  # Flatten the array of axes for easy iteration

    # Plot each period on a different subplot
    periods = da.coords["period"].values if "period" in da.coords else [None]

    for i, (ax, period) in enumerate(zip(axs, periods)):
        if period is not None:
            subset_da = da.sel(period=period)
        else:
            subset_da = da
        ax.set_title(f"Period: {period}" if period is not None else "Full Data")
        ax.axis("off")
        ax.add_feature(cfeature.LAND, edgecolor="black")
        ax.add_feature(cfeature.OCEAN, color="lightblue")
        ax.coastlines()
        ax.gridlines(draw_labels=True)
        ax.set_extent(
            [-638000, 860000, -3450000, -656000],
            ccrs.Stereographic(
                central_longitude=-45, central_latitude=90, true_scale_latitude=70
            ),
        )
        im = subset_da.plot.imshow(
            ax=ax, **kwargs, add_colorbar=False, transform=ccrs.epsg(3413)
        )
        fig.colorbar(
            im,
            ax=ax,
            extend="both",
            shrink=0.49,
            label=f"{da.attrs['long_name']} [{da.attrs['units']}]",
        )

    plt.tight_layout()
    plt.show()

[yoniv1]

Thanks, Mattia! it looks good.