Investigate HEALPix Support

[philipc2]

References

https://healpix.sourceforge.io/

https://pypi.org/project/healpix/

https://healpix.jpl.nasa.gov/pdf/intro.pdf

https://easy.gems.dkrz.de/Processing/healpix/index.html

[philipc2]

Following a similar track to what I've done with other grid formats, I've initially attempted to represent a HEALPix grid in the UGRID conventions.

The current expectation from UXarray (and the UGRID conventions) is that at least the following are required to represent a 2D unstructured grid:

• node_lon: Longitude of the corner nodes that make up each face/cell
• node_lat: Latitude of the corner nodes that make up each face/cell
• face_node_connectivity: Indices of the corner nodes that make up each face/cell

If we posses these three variables, a Grid can be constructed even if we do not support a format directly:

uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)

I used the healpy library to iterate over each pixel to obtain its boundaries using hp.boundaries and used the outputs to generate the face_node_connectivity

import uxarray as ux
import healpy as hp
import numpy as np

n_side = 4
n_pix = hp.nside2npix(n_side)  


unique_lonlat_pairs = dict()
face_node_connectivity = np.empty((n_pix, 4), dtype=ux.INT_DTYPE)

for pix in range(n_pix):
    vec = hp.boundaries(n_side, pix, step=1, nest=False).T  

    # Convert Cartesian coordinates to spherical coordinates
    lon, lat = hp.vec2ang(vec, lonlat=True)
    lon = (lon + 180) % 360 - 180
    
    cur_face_nodes = []
    
    for i in range(4):
        # 4 corner nodes for each pixel
        if lon[i] == -180.0:
            # avoid duplicates
            lon[i] = 180
        
        lonlat = (lon[i], lat[i])
        if lonlat not in unique_lonlat_pairs:
            unique_lonlat_pairs[lonlat] = len(unique_lonlat_pairs)
            
        face_node_connectivity[pix, i] = unique_lonlat_pairs[lonlat]
        
node_lon = np.empty(len(unique_lonlat_pairs), dtype=np.float32)
node_lat = np.empty(len(unique_lonlat_pairs), dtype=np.float32)

for i, lonlat in enumerate(unique_lonlat_pairs):
    node_lon[i] = lonlat[0]
    node_lat[i] = lonlat[1]

uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)

This was a very initial exploration, but the results appear to be on the right track.

What I'm the most curious about is whether a UGRID to HEALPix conversion is something that is desired by the community, and/or if it is even practical.