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.

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.

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([[674, 500, 709, 345],
       [663, 615, 500, 674],
       [500, 197, 682, 709],
       ...,
       [170,   7,  41, 235],
       [ 41, 267, 381, 696],
       [  7, 303, 267,  41]], 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.

Our original data variable now resides on our HEALPix grid.

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

We can now chain together a few functions to convert our output from a UXarray Dataset to a NetCDF file, with the grid represented in the HEALPix format.

to_dataset(): Converts the ux.UxDataArray to a ux.UxDataset
to_xarray(grid_format="HEALPix"): Converts our ux.UxDataset to a xr.Dataset encoded in the HEALPix format.
to_netcdf("psi_healpix.nc"): Saves our dataset to disk as a NetCDF file.

psi_hp.to_dataset().to_xarray(grid_format="HEALPix").to_netcdf("psi_healpix.nc")

HEALPix Equal Area Property#

HEALPix’s fundamental property is that all pixels at a given resolution have exactly the same spherical area. This “equal area” property is crucial for scientific applications such as global averaging, zonal means, and conservation properties in regridding operations.

The theoretical area of each HEALPix pixel is given by:

\[ A_{\text{pixel}} = \frac{4\pi}{N_{\text{pix}}} = \frac{4\pi}{12 \cdot 4^{\text{resolution}}} \]

where
\(N_{\text{pix}} = 12 \cdot 4^{\text{resolution}}\) is the total number of pixels at a given resolution level, and
\(4\pi\) is the total surface area of the unit sphere in steradians.

Geometric Representation Considerations#

UXarray represents HEALPix grids using Great Circle Arcs (GCAs) to define pixel boundaries, following UGRID conventions. However, this geometric representation can introduce systematic errors when computing areas numerically, potentially violating HEALPix’s equal-area property.

Note

To alleviate the impacts of this systematic differences between UXarray and HEALPix, we adjust our Grid.face_areas property to fulfill the HEALPix equal area property, making sure that all the faces in a HEALPix mesh have the same theoretical HEALPix area.

Let’s demonstrate this with HEALPix’s 12 base pixels:

import numpy as np

# Create HEALPix grid at resolution 0 (12 base pixels)
grid = ux.Grid.from_healpix(0, pixels_only=False)

# Using face_areas property ensures equal areas for HEALPix
hp_areas = grid.face_areas.values
print(f"Standard deviation: {np.std(hp_areas):.2e}")
print(f"All pixels have area: {hp_areas[0]:.6f} steradians")

Standard deviation: 2.22e-16
All pixels have area: 1.047198 steradians

Still need geometric face area calculations for a HEALPix mesh?#

For most use cases, the Grid.face_areas property provides the recommended approach for accessing face areas. However, if you specifically need geometric calculations of individual HEALPix faces as they are represented in UXarray (rather than the theoretical equal areas), you may want to access the internal computation method, Grid._compute_face_areas(). Note that this approach may not preserve HEALPix’s equal-area property due to geometric representation differences.

Look at the following example:

# For advanced use cases: access geometric calculations (may not preserve equal-area property)
hp_geometric_areas, face_jacobians = grid._compute_face_areas(
    quadrature_rule="triangular", order=4
)
print(
    f"Geometric std deviation: {hp_geometric_areas.std():.2e} (vs theoretical equal areas)"
)

Geometric std deviation: 1.50e-01 (vs theoretical equal areas)

If you are interested in futher details of the systematic errors due to geometric calculation with a clear spatial pattern, our analysis across resolution levels 0-7 shows:

Equatorial pixels (lat≈0°): +20.5% area error
Mid-latitude pixels (lat≈±42°): -9.9% area error
Maximum errors: Persist at ~10% even at fine resolutions

Future Work#

Support for the HEALPix format is still in its early stages. We welcome your suggestions and contributions! To get started, please check out our Contributors Guide.