Quickstart: Access Data from Disc at DKRZ

Browsing the STAC catalog

We create a Client from the pystac_client package. For this example we assume, that we already know the collection ID and item ID. This allows us to navigate directly through the catalog. Browsing the catalog is a topic for another guide.

To open the catalog, we need its base URL, here https://wwestac.cloud.dkrz.de/stac-fastapi-es/. Using a different URL, grants access to a different catalog.

import pystac_client
catalog = pystac_client.Client.open("https://wwestac.cloud.dkrz.de/stac-fastapi-es/")
collection = catalog.get_collection("ngc3026")
items = asset = collection.get_items()
item = collection.get_item('ngc3026_P1D_7')
print(item)

<Item id=ngc3026_P1D_7>

Opening the Dataset

As the last step, we select an asset, in other words an way to open our desired dataset. Here, we select disk. Printing the path to would allow us using this path for the data analysis. Each item can have several assets with different options to open the desired dataset.

zarr_path=item.assets['disk'].href
print(zarr_path)

file:///work/bm1235/k203123/nextgems_cycle3/experiments/ngc3026/outdata/ngc3026_P1D_7.zarr

Now we can open the dataset and inspect it. The default output gives already several information on data axis, variables etc. In this example we use the xarray package, which gives a common interface and supports a wide variety of file formats. xarray also gives extensive features for data processing. For the usage, you need to install the python zarr package. The required import is handeled by xarray.

import xarray as xr
ds = xr.open_zarr(zarr_path)
ds

<xarray.Dataset>
Dimensions:                              (time: 370, depth_half: 129,
                                          cell: 196608, level_full: 90, crs: 1,
                                          depth_full: 128,
                                          soil_depth_water_level: 5,
                                          level_half: 91,
                                          soil_depth_energy_level: 5)
Coordinates:
  * crs                                  (crs) float32 nan
  * depth_full                           (depth_full) float32 1.0 ... 5.904e+03
  * depth_half                           (depth_half) float32 0.0 ... 6.003e+03
  * level_full                           (level_full) int32 1 2 3 4 ... 88 89 90
  * level_half                           (level_half) int32 1 2 3 4 ... 89 90 91
  * soil_depth_energy_level              (soil_depth_energy_level) float32 0....
  * soil_depth_water_level               (soil_depth_water_level) float32 0.0...
  * time                                 (time) datetime64[ns] 2020-01-21 ......
Dimensions without coordinates: cell
Data variables: (12/88)
    a_tracer_v_to                        (time, depth_half, cell) float32 dask.array<chunksize=(1, 33, 196608), meta=np.ndarray>
    atmos_fluxes_frshflux_evaporation    (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    atmos_fluxes_frshflux_precipitation  (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    atmos_fluxes_frshflux_runoff         (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    atmos_fluxes_frshflux_snowfall       (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    atmos_fluxes_heatflux_latent         (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    ...                                   ...
    va                                   (time, level_full, cell) float32 dask.array<chunksize=(1, 30, 196608), meta=np.ndarray>
    vas                                  (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    w                                    (time, depth_half, cell) float32 dask.array<chunksize=(1, 33, 196608), meta=np.ndarray>
    wa_phy                               (time, level_half, cell) float32 dask.array<chunksize=(1, 31, 196608), meta=np.ndarray>
    wind_speed_10m                       (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>
    zos                                  (time, cell) float32 dask.array<chunksize=(1, 196608), meta=np.ndarray>

Visualization

We create two functions:

  1. nnshow to visualize healpix grid on a lat-lon representation as used by matplotlib for the visualization

  2. worldmap to display the selected field on top of a worldmap showing the contours of the continents.

As a first step, we include a number of packages, required for the processing and plotting. In particular cartopy has additional external dependencies, that might not be installed by the package manager (e.g. pip).

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

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)

Plotting

Last but not least we visualize the 2m air temperature (tas) from this ICON dataset.

worldmap(ds.tas.isel(time=0), cmap=cmocean.cm.thermal)
../../_images/ac2c2562692c0e438c7f54a8c54ad1a3eecbb0b9642e36eebe380710236e5873.png