Add HEALPix Support

[philipc2]

Closes #1097 #1028

Overview

• Implements a Grid.from_healpix(zoom, nest, pixels_only) classmethod that constructs a UGRID-compliant representation of a HEALPix Grid, and a ``UxDataset.from_healpix()` classmethod to load in data.
• The HEALPix grid can be loaded as "coordinates-only", meaning that only the pixels are stored, with boundaries (i.e. face_node_connectivity, node_lonandnode_lat` computed when the user requires them (i.e. for plotting)
• User Guide section on working with HEALPix Grids

Examples

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[rajeeja]

how do you get healpix mesh data, what is the format? do you have an example?

[philipc2]

how do you get healpix mesh data, what is the format? do you have an example?

I've put a few examples in #1028 which are a good read

[erogluorhan]

Can you link the equation to a HEALPix resource, I think this one: https://easy.gems.dkrz.de/Processing/healpix/index.html#healpix-spatial-resolution ?

Reply via ReviewNB

[philipc2]

Sure! I need to fix the latex as well.

[erogluorhan]

A few things:

1. The first sentence seems incomplete. Maybe just removing "Since" will fix it.
2. "...requires explicitly defined connectivity...", i.e. "explicitly" instead of "explicit"
3. It is a bit unclear how this paragraph connects to the previous code cell and the use of pixels_only in the following code cell. Also, the message of this paragraph (i.e. The most important one is the face_node_connectivity) and how the current implementation is done (i.e. defaulting to pixels only) seems a bit disconnected. Should it be something like this: "This above code creates a UXarray Grid with only HEALPix pixels' coordinates, i.e. not boundaries. However, much of UXarray's functionality requires explicitly defined connectivity arrays. The most important one is the face_node_connectivity, which contains the indices of the boundaries of each pixel (i.e. corner nodes? in this implementation). Generating the face_node_connectivity also constructs the node_lat and node_lon coordinates, corresponding to the latitude and longitude of the boundaries of each pixel. We can enforce construction of these by using pixels_only as follows:"

Reply via ReviewNB

[erogluorhan]

Slight rephrase: "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."

Reply via ReviewNB

[erogluorhan]

A bit of rephrase maybe: "Once loaded in, UXarray does not care about the underlying grid format as it represents all formats in its unified UGRID-like convention; hence, all existing functionality is supported on HEALPix through the same UXarray interface, such as plotting."

Reply via ReviewNB

[philipc2]

Good idea.

[florianziemen]

Just a quick comment for now:
Looks good at first sight. I'll look deeper next week.
For now this assumes a full global grid. That's a good start, and I guess, we can add regional cut-outs in a later PR.

[philipc2]

Just a quick comment for now: Looks good at first sight. I'll look deeper next week. For now this assumes a full global grid. That's a good start, and I guess, we can add regional cut-outs in a later PR.

Thanks for the initial comment!

Do you have any example datasets that have region cut-outs?

[florianziemen]

This works with the easygems package:

import xarray as xr 
import numpy as np
from easygems import healpix as egh

da = xr.DataArray (np.arange (48,0.,-1), dims = ("cell",))
crs = xr.DataArray(0, coords = dict (crs=np.nan))
crs.attrs = dict (grid_mapping_name = 'healpix', healpix_nside=2**1, healpix_order='nest')
da.attrs['grid_mapping'] = "crs"

ds = xr.Dataset(dict (data= da), coords = dict (cell = np.arange(48),  crs = crs))
subset = ds.isel(cell=slice(3,41))

egh.healpix_show(subset['data'])

Maybe @lkluft has more thoughts on this. He's been working with regional grids.

[philipc2]

I've put together some revisions to the user guide, should be much cleaner now.

[lkluft]

Maybe @lkluft has more thoughts on this. He's been working with regional grids.

To give some more context:

A "regional grid" in our test cases is simply a slice of a full global HEALPix grid, which we store along with its global indices. This way it is still possible to facilitate the fast ang2pix/pix2ang mappings provided by the HEALPix library.

When it comes to plotting, this offers several possibilities. In easygems an algorithm fills the pixel space of the screen with the closest values in the input array. As the values are addressed by their global indices, this works out-of-the-box with our regional slices. Alternatively, one could recreate the full HEALPix grid and fill the regional extent with the available data, but this would increase memory usage.

[rajeeja]

Good job.