Remapping a dataset to HEALPix

In this tutorial we will remap a dataset from the native ICON grid to HEALPix. The dataset used as an example is from the EERIE project and is available online. It consists of five 2d variables and covers 23 years with a 6 hour time step.

import healpix as hp
import numpy as np
import xarray as xr

import easygems.healpix as egh
import easygems.remap as egr


ds = xr.open_dataset(
    "https://eerie.cloud.dkrz.de/datasets/icon-esm-er.eerie-control-1950.v20240618.atmos.native.2d_6h_inst/kerchunk",
    engine="zarr",
    chunks={},
)
ds

<xarray.Dataset> Size: 8TB
Dimensions:             (ncells: 5242880, height: 1, height_2: 1, time: 73051)
Coordinates:
    cell_sea_land_mask  (ncells) float64 42MB dask.array<chunksize=(5242880,), meta=np.ndarray>
  * height              (height) float64 8B 10.0
  * height_2            (height_2) float64 8B 2.0
    lat                 (ncells) float64 42MB dask.array<chunksize=(5242880,), meta=np.ndarray>
    lon                 (ncells) float64 42MB dask.array<chunksize=(5242880,), meta=np.ndarray>
  * time                (time) datetime64[ns] 584kB 1991-01-01T05:59:59 ... 2...
Dimensions without coordinates: ncells
Data variables:
    hus2m               (time, height_2, ncells) float32 2TB dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>
    psl                 (time, ncells) float32 2TB dask.array<chunksize=(1, 5242880), meta=np.ndarray>
    ts                  (time, ncells) float32 2TB dask.array<chunksize=(1, 5242880), meta=np.ndarray>
    uas                 (time, height, ncells) float32 2TB dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>
    vas                 (time, height, ncells) float32 2TB dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>
Attributes: (12/31)
    Conventions:           CF-1.7 CMIP-6.2
    activity_id:           EERIE
    data_specs_version:    01.00.32
    forcing_index:         1
    initialization_index:  1
    license:               EERIE model data produced by MPI-M is licensed und...
    ...                    ...
    parent_activity_id:    EERIE
    sub_experiment_id:     none
    experiment:            coupled control with fixed 1950's forcing (HighRes...
    source:                ICON-ESM-ER (2023): \naerosol: none, prescribed MA...
    institution:           Max Planck Institute for Meteorology, Hamburg 2014...
    sub_experiment:        none

xarray.Dataset

• Dimensions:
  • ncells: 5242880
  • height: 1
  • height_2: 1
  • time: 73051
• Coordinates: (6)
  • cell_sea_land_mask
    (ncells)
    float64
    dask.array<chunksize=(5242880,), meta=np.ndarray>

    grid_type :

    unstructured

    long_name :

    sea (-2 inner, -1 boundary) land (2 inner, 1 boundary) mask for the cell

    number_of_grid_in_reference :

    1

    units :

    2,1,-1,-

    Array Chunk Bytes 40.00 MiB 40.00 MiB Shape (5242880,) (5242880,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • height
    (height)
    float64
    10.0

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    array([10.])

  • height_2
    (height_2)
    float64
    2.0

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    array([2.])

  • lat
    (ncells)
    float64
    dask.array<chunksize=(5242880,), meta=np.ndarray>

    bounds :

    clat_vertices

    long_name :

    center latitude

    standard_name :

    grid_latitude

    units :

    radian

    Array Chunk Bytes 40.00 MiB 40.00 MiB Shape (5242880,) (5242880,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • lon
    (ncells)
    float64
    dask.array<chunksize=(5242880,), meta=np.ndarray>

    bounds :

    clon_vertices

    long_name :

    center longitude

    standard_name :

    grid_longitude

    units :

    radian

    Array Chunk Bytes 40.00 MiB 40.00 MiB Shape (5242880,) (5242880,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • time
    (time)
    datetime64[ns]
    1991-01-01T05:59:59 ... 2040-12-...

    standard_name :

    time

    array(['1991-01-01T05:59:59.000000000', '1991-01-01T11:59:59.000000000',
           '1991-01-01T17:59:59.000000000', ..., '2040-12-31T05:59:59.000000000',
           '2040-12-31T11:59:59.000000000', '2040-12-31T17:59:59.000000000'],
          shape=(73051,), dtype='datetime64[ns]')

• Data variables: (5)
  • hus2m
    (time, height_2, ncells)
    float32
    dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    specific humidity in 2m

    number_of_grid_in_reference :

    1

    param :

    6.0.0

    standard_name :

    qv2m

    units :

    kg kg-1

    Array Chunk Bytes 1.39 TiB 20.00 MiB Shape (73051, 1, 5242880) (1, 1, 5242880) Dask graph 73051 chunks in 2 graph layers Data type float32 numpy.ndarray

  • psl
    (time, ncells)
    float32
    dask.array<chunksize=(1, 5242880), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    mean sea level pressure

    number_of_grid_in_reference :

    1

    param :

    1.3.0

    standard_name :

    mean sea level pressure

    units :

    Pa

    Array Chunk Bytes 1.39 TiB 20.00 MiB Shape (73051, 5242880) (1, 5242880) Dask graph 73051 chunks in 2 graph layers Data type float32 numpy.ndarray

  • ts
    (time, ncells)
    float32
    dask.array<chunksize=(1, 5242880), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    surface temperature

    number_of_grid_in_reference :

    1

    param :

    0.0.0

    standard_name :

    surface_temperature

    units :

    K

    Array Chunk Bytes 1.39 TiB 20.00 MiB Shape (73051, 5242880) (1, 5242880) Dask graph 73051 chunks in 2 graph layers Data type float32 numpy.ndarray

  • uas
    (time, height, ncells)
    float32
    dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    zonal wind in 10m

    number_of_grid_in_reference :

    1

    param :

    2.2.0

    standard_name :

    uas

    units :

    m s-1

    Array Chunk Bytes 1.39 TiB 20.00 MiB Shape (73051, 1, 5242880) (1, 1, 5242880) Dask graph 73051 chunks in 2 graph layers Data type float32 numpy.ndarray

  • vas
    (time, height, ncells)
    float32
    dask.array<chunksize=(1, 1, 5242880), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    meridional wind in 10m

    number_of_grid_in_reference :

    1

    param :

    3.2.0

    standard_name :

    vas

    units :

    m s-1

    Array Chunk Bytes 1.39 TiB 20.00 MiB Shape (73051, 1, 5242880) (1, 1, 5242880) Dask graph 73051 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • height
    PandasIndex

    PandasIndex(Index([10.0], dtype='float64', name='height'))

  • height_2
    PandasIndex

    PandasIndex(Index([2.0], dtype='float64', name='height_2'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1991-01-01 05:59:59', '1991-01-01 11:59:59',
                   '1991-01-01 17:59:59', '1991-01-01 23:59:59',
                   '1991-01-02 05:59:59', '1991-01-02 11:59:59',
                   '1991-01-02 17:59:59', '1991-01-02 23:59:59',
                   '1991-01-03 05:59:59', '1991-01-03 11:59:59',
                   ...
                   '2040-12-29 11:59:59', '2040-12-29 17:59:59',
                   '2040-12-29 23:59:59', '2040-12-30 05:59:59',
                   '2040-12-30 11:59:59', '2040-12-30 17:59:59',
                   '2040-12-30 23:59:59', '2040-12-31 05:59:59',
                   '2040-12-31 11:59:59', '2040-12-31 17:59:59'],
                  dtype='datetime64[ns]', name='time', length=73051, freq=None))

• Attributes: (31)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  EERIE

  data_specs_version :

  01.00.32

  forcing_index :

  1

  initialization_index :

  1

  license :

  EERIE model data produced by MPI-M is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses). The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  EERIE

  project_id :

  EERIE

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  source_type :

  AOGCM

  variant_label :

  r1i1p1f1

  creation_date :

  2024-12-16T16:05:56

  institution_id :

  MPI-M

  source_id :

  ICON-ESM-ER

  experiment_id :

  eerie-control-1950

  contact :

  juergen.kroeger@mpimet.mpg.de

  references :

  Hohenegger et al., ICON-Sapphire: simulating the components of the Earth system and their interactions at kilometer and subkilometer scales. Geosci. Model Dev., 16, 779–811, 2023, https://doi.org/10.5194/gmd-16-779-2023

  nominal_resolution :

  10 km

  version_id :

  v20240618

  realm :

  atmos

  grid :

  native

  grid_label :

  gn

  parent_experiment_id :

  eerie-spinup-1950

  parent_activity_id :

  EERIE

  sub_experiment_id :

  none

  experiment :

  coupled control with fixed 1950's forcing (HighResMIP equivalent of pre-industrial control)

  source :

  ICON-ESM-ER (2023): aerosol: none, prescribed MACv2-SP atmos: ICON-A (icosahedral/triangles; 10 km; 90 levels; top level 80 km) atmosChem: none land: JSBACH4.20 landIce: none/prescribed ocean: ICON-O (icosahedral/triangles; 5 km; 72 levels; top grid cell 0-2 m) ocnBgchem: none seaIce: unnamed (thermodynamic (Semtner zero-layer) dynamic (Hibler 79) sea ice model)

  institution :

  Max Planck Institute for Meteorology, Hamburg 20146, Germany

  sub_experiment :

  none

The first step is to create a HEALPix grid that is close to the resolution of our source grid. Here we will choose an order of 9 (also known as “zoom”) and “nest” ordering.

order = zoom = 9
nside = hp.order2nside(order)
npix = hp.nside2npix(nside)

hp_lon, hp_lat = hp.pix2ang(nside=nside, ipix=np.arange(npix), lonlat=True, nest=True)
hp_lon = (hp_lon + 180) % 360 - 180  # [-180, 180)
hp_lon += 360 / (4 * nside) / 4  # shift quarter-width

Next, we can use our defined source and target grids to compute interpolation weights. The easygems package provides a function to compute these weights using the Delaunay triangulation method.

weights = egr.compute_weights_delaunay((ds.lon, ds.lat), (hp_lon, hp_lat))

# You can also save the calculated weights for future use
# weights.to_netcdf("healpix_weights.nc")

These weights can be applied to single fields directly:

egh.healpix_show(egr.apply_weights(ds.ts.isel(time=0), **weights))

<matplotlib.image.AxesImage at 0x7efbe489b980>

/usr/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

<Figure size 640x480 with 0 Axes>

Alternatively, we can use xarray’s apply_ufunc() function to lift the function onto a full dataset. This requires a coupled of additional information from the user, e.g. the input dimension along which the function should be applied, and the resulting output dimensions name and size.

ds_remap = xr.apply_ufunc(
    egr.apply_weights,
    ds,
    kwargs=weights,
    keep_attrs=True,
    input_core_dims=[["ncells"]],
    output_core_dims=[["cell"]],
    output_dtypes=["f4"],
    vectorize=True,
    dask="parallelized",
    dask_gufunc_kwargs={
        "output_sizes": {"cell": npix},
    },
)
ds_remap

<xarray.Dataset> Size: 5TB
Dimensions:   (time: 73051, height_2: 1, cell: 3145728, height: 1)
Coordinates:
  * height    (height) float64 8B 10.0
  * height_2  (height_2) float64 8B 2.0
  * time      (time) datetime64[ns] 584kB 1991-01-01T05:59:59 ... 2040-12-31T...
Dimensions without coordinates: cell
Data variables:
    hus2m     (time, height_2, cell) float32 919GB dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>
    psl       (time, cell) float32 919GB dask.array<chunksize=(1, 3145728), meta=np.ndarray>
    ts        (time, cell) float32 919GB dask.array<chunksize=(1, 3145728), meta=np.ndarray>
    uas       (time, height, cell) float32 919GB dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>
    vas       (time, height, cell) float32 919GB dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>
Attributes: (12/31)
    Conventions:           CF-1.7 CMIP-6.2
    activity_id:           EERIE
    data_specs_version:    01.00.32
    forcing_index:         1
    initialization_index:  1
    license:               EERIE model data produced by MPI-M is licensed und...
    ...                    ...
    parent_activity_id:    EERIE
    sub_experiment_id:     none
    experiment:            coupled control with fixed 1950's forcing (HighRes...
    source:                ICON-ESM-ER (2023): \naerosol: none, prescribed MA...
    institution:           Max Planck Institute for Meteorology, Hamburg 2014...
    sub_experiment:        none

xarray.Dataset

• Dimensions:
  • time: 73051
  • height_2: 1
  • cell: 3145728
  • height: 1
• Coordinates: (3)
  • height
    (height)
    float64
    10.0

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    array([10.])

  • height_2
    (height_2)
    float64
    2.0

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    array([2.])

  • time
    (time)
    datetime64[ns]
    1991-01-01T05:59:59 ... 2040-12-...

    standard_name :

    time

    array(['1991-01-01T05:59:59.000000000', '1991-01-01T11:59:59.000000000',
           '1991-01-01T17:59:59.000000000', ..., '2040-12-31T05:59:59.000000000',
           '2040-12-31T11:59:59.000000000', '2040-12-31T17:59:59.000000000'],
          shape=(73051,), dtype='datetime64[ns]')

• Data variables: (5)
  • hus2m
    (time, height_2, cell)
    float32
    dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    specific humidity in 2m

    number_of_grid_in_reference :

    1

    param :

    6.0.0

    standard_name :

    qv2m

    units :

    kg kg-1

    Array Chunk Bytes 856.07 GiB 12.00 MiB Shape (73051, 1, 3145728) (1, 1, 3145728) Dask graph 73051 chunks in 5 graph layers Data type float32 numpy.ndarray

  • psl
    (time, cell)
    float32
    dask.array<chunksize=(1, 3145728), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    mean sea level pressure

    number_of_grid_in_reference :

    1

    param :

    1.3.0

    standard_name :

    mean sea level pressure

    units :

    Pa

    Array Chunk Bytes 856.07 GiB 12.00 MiB Shape (73051, 3145728) (1, 3145728) Dask graph 73051 chunks in 5 graph layers Data type float32 numpy.ndarray

  • ts
    (time, cell)
    float32
    dask.array<chunksize=(1, 3145728), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    surface temperature

    number_of_grid_in_reference :

    1

    param :

    0.0.0

    standard_name :

    surface_temperature

    units :

    K

    Array Chunk Bytes 856.07 GiB 12.00 MiB Shape (73051, 3145728) (1, 3145728) Dask graph 73051 chunks in 5 graph layers Data type float32 numpy.ndarray

  • uas
    (time, height, cell)
    float32
    dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    zonal wind in 10m

    number_of_grid_in_reference :

    1

    param :

    2.2.0

    standard_name :

    uas

    units :

    m s-1

    Array Chunk Bytes 856.07 GiB 12.00 MiB Shape (73051, 1, 3145728) (1, 1, 3145728) Dask graph 73051 chunks in 5 graph layers Data type float32 numpy.ndarray

  • vas
    (time, height, cell)
    float32
    dask.array<chunksize=(1, 1, 3145728), meta=np.ndarray>

    CDI_grid_type :

    unstructured

    long_name :

    meridional wind in 10m

    number_of_grid_in_reference :

    1

    param :

    3.2.0

    standard_name :

    vas

    units :

    m s-1

    Array Chunk Bytes 856.07 GiB 12.00 MiB Shape (73051, 1, 3145728) (1, 1, 3145728) Dask graph 73051 chunks in 5 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • height
    PandasIndex

    PandasIndex(Index([10.0], dtype='float64', name='height'))

  • height_2
    PandasIndex

    PandasIndex(Index([2.0], dtype='float64', name='height_2'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1991-01-01 05:59:59', '1991-01-01 11:59:59',
                   '1991-01-01 17:59:59', '1991-01-01 23:59:59',
                   '1991-01-02 05:59:59', '1991-01-02 11:59:59',
                   '1991-01-02 17:59:59', '1991-01-02 23:59:59',
                   '1991-01-03 05:59:59', '1991-01-03 11:59:59',
                   ...
                   '2040-12-29 11:59:59', '2040-12-29 17:59:59',
                   '2040-12-29 23:59:59', '2040-12-30 05:59:59',
                   '2040-12-30 11:59:59', '2040-12-30 17:59:59',
                   '2040-12-30 23:59:59', '2040-12-31 05:59:59',
                   '2040-12-31 11:59:59', '2040-12-31 17:59:59'],
                  dtype='datetime64[ns]', name='time', length=73051, freq=None))

• Attributes: (31)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  EERIE

  data_specs_version :

  01.00.32

  forcing_index :

  1

  initialization_index :

  1

  license :

  EERIE model data produced by MPI-M is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses). The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  EERIE

  project_id :

  EERIE

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  source_type :

  AOGCM

  variant_label :

  r1i1p1f1

  creation_date :

  2024-12-16T16:05:56

  institution_id :

  MPI-M

  source_id :

  ICON-ESM-ER

  experiment_id :

  eerie-control-1950

  contact :

  juergen.kroeger@mpimet.mpg.de

  references :

  Hohenegger et al., ICON-Sapphire: simulating the components of the Earth system and their interactions at kilometer and subkilometer scales. Geosci. Model Dev., 16, 779–811, 2023, https://doi.org/10.5194/gmd-16-779-2023

  nominal_resolution :

  10 km

  version_id :

  v20240618

  realm :

  atmos

  grid :

  native

  grid_label :

  gn

  parent_experiment_id :

  eerie-spinup-1950

  parent_activity_id :

  EERIE

  sub_experiment_id :

  none

  experiment :

  coupled control with fixed 1950's forcing (HighResMIP equivalent of pre-industrial control)

  source :

  ICON-ESM-ER (2023): aerosol: none, prescribed MACv2-SP atmos: ICON-A (icosahedral/triangles; 10 km; 90 levels; top level 80 km) atmosChem: none land: JSBACH4.20 landIce: none/prescribed ocean: ICON-O (icosahedral/triangles; 5 km; 72 levels; top grid cell 0-2 m) ocnBgchem: none seaIce: unnamed (thermodynamic (Semtner zero-layer) dynamic (Hibler 79) sea ice model)

  institution :

  Max Planck Institute for Meteorology, Hamburg 20146, Germany

  sub_experiment :

  none

The resulting dataset can then be used as usual and the remapping is performed on demand.

egh.healpix_show(ds_remap.ts.sel(time="1991-01-01").mean("time"))

<matplotlib.image.AxesImage at 0x7efbedb5dc40>

<Figure size 640x480 with 0 Axes>

Storing the coordinate reference system

It is good practice to store map projection information in the Coordinate Reference Systems (CRS). By making the crs coordinate non-dimensional it will “stick” to the dataset even if individual variables are subselected.

ds_remap["crs"] = xr.DataArray(
    name="crs",
    data=0,
    attrs={
        "grid_mapping_name": "healpix",
        "healpix_nside": 2**zoom,
        "healpix_order": "nest",
    },
)