Working with HEALPix Data

Working with HEALPix Data#

In this notebook, we showcase how UXarray can be used with HEALPix data.

  • A brief overview into HEALPix

  • Representing a HEALPix grid in the UGRID conventions

  • Loading data that resides on HEALPix Grids

import cartopy.crs as ccrs

import uxarray as ux

What is HEALPix?#

The Hierarchical Equal Area isoLatitude Pixelisation (HEALPix) algorithm is a method for the pixelisation of the 2-sphere. It has three defining qualities.

  • The sphere is hierarchically tessellated into curvilinear quadrilaterals

  • Areas of all pixels at a given resolution are identical

  • Pixels are distributed on lines of constant latitude

Note

For more information on HEALPix, you can refer to this excellent overview in DKRZ’s Easy Gems documentation.

Representing a HEALPix Grid in the UGRID Conventions#

Most datasets loaded into UXarray come with a dedicated grid file. However, for HEALPix datasets, the underlying grid structure is inferred from the number of cells.

For a HEALPix grid in a nested layout, the number of cells is defined by:

\[ N_{\text{cells}} = 12 \cdot z^4, \]

where z is the zoom level. A higher zoom level increases the effective resolution of the grid.

The ux.Grid.from_healpix(zoom) classmethod can be used to construct a ux.Grid instance representing a HEALPix grid in the UGRID conventions.

ux.Grid.from_healpix(zoom=3)

<uxarray.Grid>
Original Grid Type: HEALPix
Grid Dimensions:
  * n_face: 768
Grid Coordinates (Spherical):
  * face_lon: (768,)
  * face_lat: (768,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
Grid Descriptor Variables:

By default, our Grid will only contain the face coordinate values, which represent the HEALPix pixel locations. We can see this above with only the face_lon and face_lat variables populated.

However, much of UXarray functionality requires the boundaries of each cell to be defined. This is represented using the following variables:

  • node_lon: Longitude of the corners of each cell

  • node_lat: Latitude of the corners of each cell

  • face_node_connectivity: The indices of the nodes that surround each face

An optional parameter, pixels_only, is provided which can be set to True if you require the boundaries to be computed upfront.

ux.Grid.from_healpix(zoom=3, pixels_only=False)

<uxarray.Grid>
Original Grid Type: HEALPix
Grid Dimensions:
  * n_node: 770
  * n_face: 768
  * n_max_face_nodes: 4
Grid Coordinates (Spherical):
  * node_lon: (770,)
  * node_lat: (770,)
  * face_lon: (768,)
  * face_lat: (768,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (768, 4)
Grid Descriptor Variables:

However, even if you load in a HEALPix grid specifying that you do not want the connectivity upfront, they can still be constructed when desired because of UXarray’s internal design.

ux.Grid.from_healpix(zoom=3, pixels_only=True).face_node_connectivity

<xarray.DataArray 'face_node_connectivity' (n_face: 768, n_max_face_nodes: 4)> Size: 25kB
array([[255, 193, 662, 219],
       [267, 594, 193, 255],
       [193, 427, 393, 662],
       ...,
       [338, 112, 541,  86],
       [541, 448, 319,  90],
       [112, 406, 448, 541]], shape=(768, 4))
Dimensions without coordinates: n_face, n_max_face_nodes
Attributes:
    cf_role:      face_node_connectivity
    long name:    Maps every face to its corner nodes.
    start_index:  0
    _FillValue:   -9223372036854775808

Grid-Agnostic Functionality#

Once loaded in, UXarray does not care about the underlying grid format, as it represents all formats in its unified UGRID-like convention. This means that all existing functionality is supported through the same UXarray inferface. Below are basic examples of plotting & remapping.

Plotting#

By using UXarray’s plotting functionality, we can observe how increasing the zoom level leads to a higher-resolution grid.

(
    ux.Grid.from_healpix(zoom=2).plot(
        periodic_elements="split", projection=ccrs.Orthographic(), title="zoom=2"
    )
    + ux.Grid.from_healpix(zoom=3).plot(
        periodic_elements="split", projection=ccrs.Orthographic(), title="zoom=3"
    )
    + ux.Grid.from_healpix(zoom=4).plot(
        periodic_elements="split", projection=ccrs.Orthographic(), title="zoom=4"
    )
).cols(1)

Remapping#

Below is a simple example of remapping UGRID data residing on a CSne30 grid to a HEALPix grid with a zoom level of 5.

source_uxds = ux.open_dataset(
    "../../test/meshfiles/ugrid/outCSne30/outCSne30.ug",
    "../../test/meshfiles/ugrid/outCSne30/outCSne30_vortex.nc",
)

# source data & grid
source_uxds["psi"].plot(
    cmap="inferno", projection=ccrs.Orthographic(), title="Source Data"
)

# destination grid
hp_grid = ux.Grid.from_healpix(zoom=5)

Before remapping, we can plot our Source and Destination grids.

source_uxds.uxgrid.plot(
    projection=ccrs.Orthographic(), title="Source Grid (CSne30)"
) + hp_grid.plot(projection=ccrs.Orthographic(), title="Destination Grid (HEALPix)")

We can now perform our remapping. In this case, we apply a simple nearest neighbor remapping.

psi_hp = source_uxds["psi"].remap.nearest_neighbor(destination_grid=hp_grid)
psi_hp

/home/docs/checkouts/readthedocs.org/user_builds/uxarray/conda/latest/lib/python3.13/site-packages/uxarray/grid/coordinates.py:254: UserWarning: This cannot be guaranteed to work correctly on concave polygons
  warnings.warn("This cannot be guaranteed to work correctly on concave polygons")

<xarray.UxDataArray 'psi' (n_face: 12288)> Size: 98kB
array([0.63179368, 0.63440432, 0.6660493 , ..., 1.42462846, 1.43725128,
       1.43774618], shape=(12288,))
Dimensions without coordinates: n_face

Our original data variable now resides on our HEALPix grid.

psi_hp.plot(cmap="inferno", projection=ccrs.Orthographic(), title="Remapped Data")