UXarray for Advanced HEALPix Analysis & Visualization¶

In this section, you’ll learn:¶

Using the uxarray package to perform advanced analysis operators over HEALPix data such as non-conservative zonal means, etc.

Prerequisites¶

Concepts	Importance	Notes
UXarray	Necessary
HEALPix overview	Necessary

Time to learn: 30 minutes

import cartopy.crs as ccrs
import intake
import uxarray as ux

Open data catalog¶

Let us open the online catalog from the WCRP’s Digital Earths Global Hackathon 2025 catalog repository using intake and read the output of the ICON run d3hp003, which is stored in the HEALPix format:

cat_url = "https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml"
cat = intake.open_catalog(cat_url)
model_run = cat.online.icon_d3hp003

/home/runner/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/intake/catalog/utils.py:173: UserWarning: Shell command not executed due to getshell=False
  warnings.warn("Shell command not executed due to getshell=False")
/home/runner/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/intake/catalog/utils.py:182: UserWarning: Shell command not executed due to getshell=False
  warnings.warn("Shell command not executed due to getshell=False")

We can look into the highest possible resolution level allowed in this dataset at zoom level = 9 as Xarray.Dataset:

ds = model_run(zoom=9, time="P1D").to_dask()

Create UXarray Datasets from HEALPix¶

We can use UXarray’s from_healpix API as follows to open a HEALPix grid from xarray.Dataset:

uxds = ux.UxDataset.from_healpix(ds)
uxds

Data variable of interest¶

Then let us pick a variable, the surface temperature, from the dataset, which will give us an uxarray.UxDataArray:

uxda = uxds["ts"]
uxda

Global mean and plot¶

Computing the global surface temperature mean (at the first timestep) and also having a quick plot of it would be a good idea to have as references to compare the upcoming analyses & visualizations to them:

print(
    "Global surface temperature average on ", uxda.time[0].values, ": ", uxda.isel(time=0).mean().values, " K"
)

Global surface temperature average on  2020-01-02T00:00:00.000000000 :  286.9924  K

%%time

projection = ccrs.Robinson()

uxda.isel(time=0).plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global surface temperature (Polygon raster)",
    width=700,
);

WARNING:param.Parameterized: Use method 'warning' via param namespace

WARNING:param.main: pandas could not register all extension types imports failed with the following error: cannot import name 'ABCIndexClass' from 'pandas.core.dtypes.generic' (/home/runner/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/pandas/core/dtypes/generic.py)

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

WARNING:param.Dimension: Use method 'get_param_values' via param namespace

CPU times: user 88.4 ms, sys: 7.04 ms, total: 95.5 ms
Wall time: 95.3 ms

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[7], line 1
----> 1 get_ipython().run_cell_magic('time', '', '\nprojection = ccrs.Robinson()\n\nuxda.isel(time=0).plot(\n    projection=projection,\n    cmap="inferno",\n    features=["borders", "coastline"],\n    title="Global surface temperature (Polygon raster)",\n    width=700,\n);\n')

File <timed exec>:3
      1 'Could not get source, probably due dynamically evaluated source code.'

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/uxarray/core/dataarray.py:2236, in UxDataArray.__getattribute__(self, name)
   2233     return method
   2235 # For all other attributes, use the default behavior
-> 2236 return super().__getattribute__(name)

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/xarray/core/utils.py:1178, in UncachedAccessor.__get__(self, obj, cls)
   1175 if obj is None:
   1176     return self._accessor
-> 1178 return self._accessor(obj)

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/uxarray/plot/accessor.py:375, in UxDataArrayPlotAccessor.__init__(self, uxda)
    374 def __init__(self, uxda: UxDataArray) -> None:
--> 375     _ensure_hvplot_imported()
    376     self._uxda = uxda

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/uxarray/plot/accessor.py:40, in _ensure_hvplot_imported()
     31 global _IMPORTED_HVPLOT
     32 if not _IMPORTED_HVPLOT:
     33     # workaround for hvplot issue #1735;
     34     #  import hvplot.pandas and hvplot.xarray always adjust the hvplot.extension().
   (...)     38     #  Store.current_backend defaults to "matplotlib" while Store.registry is empty,
     39     #  and restoring that unloaded backend would break rendering.
---> 40     from holoviews import Store as _store
     42     _backend_orig = _store.current_backend
     43     import hvplot.pandas

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/holoviews/__init__.py:12
      6 import param
      9 __version__ = str(param.version.Version(fpath=__file__, archive_commit="$Format:%h$",
     10                                         reponame="holoviews"))
---> 12 from . import util # noqa (API import)
     13 from .core import archive, config                        # noqa (API import)
     14 from .core.boundingregion import BoundingBox             # noqa (API import)

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/holoviews/util/__init__.py:14
     11 import param
     12 from pyviz_comms import extension as _pyviz_extension
---> 14 from ..core import DynamicMap, HoloMap, Dimensioned, ViewableElement, StoreOptions, Store
     15 from ..core.options import options_policy, Keywords, Options
     16 from ..core.operation import Operation

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/holoviews/core/__init__.py:4
      1 from datetime import date, datetime
      3 from .boundingregion import *  # noqa (API import)
----> 4 from .data import *            # noqa (API import)
      5 from .dimension import *       # noqa (API import)
      6 from .element import *         # noqa (API import)

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/holoviews/core/data/__init__.py:20
     18 from .array import ArrayInterface
     19 from .dictionary import DictInterface
---> 20 from .grid import GridInterface
     21 from .multipath import MultiInterface         # noqa (API import)
     22 from .image import ImageInterface             # noqa (API import)

File ~/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/site-packages/holoviews/core/data/grid.py:5
      3 import sys
      4 import datetime as dt
----> 5 from collections import OrderedDict, defaultdict, Iterable
      7 try:
      8     import itertools.izip as zip

ImportError: cannot import name 'Iterable' from 'collections' (/home/runner/micromamba/envs/healpix-cookbook-dev2/lib/python3.14/collections/__init__.py)

Rasterized point plots¶

When working with a higher-resolution dataset at a global scale, it’s not always practical to render each cell as a polygon. Instead, we can rasterize the center of each pixel.

%%time
projection = ccrs.Robinson()

# Controls the size of each pixel (smaller value leads to larger pixels)
pixel_ratio = 0.5

uxda.isel(time=0).plot.points(
    projection=projection,
    rasterize=True,
    dynamic=False,
    width=1000,
    height=500,
    pixel_ratio=pixel_ratio,
    cmap="inferno",
    title=f"Global surface temperature (Point raster), pixel_ratio={pixel_ratio}",
);

If we decrease the size of each pixel (by setting the pixel ratio to a higher value), we can start to see missing values, which is due to a lower density of points near the poles, leading to some pixels not containing any of our original points.

Because of this, it’s useful to try a few pixel_ratio values and see which one works best for your given resolution.

projection = ccrs.Robinson()

# Controls the size of each pixel (smaller value leads to larger pixels)
pixel_ratio = 2.0

uxda.isel(time=0).plot.points(
    projection=projection,
    rasterize=True,
    dynamic=False,
    width=1000,
    height=500,
    pixel_ratio=pixel_ratio,
    cmap="inferno",
    title=f"Global surface temperature (Point raster) with a bad pixel size selection, pixel_ratio={pixel_ratio}",
);

Cross-sections¶

We can look at constant latitude/longitude cross-sections of an uxarray.UxDataArray:

boulder_lat = 40.0190


# With fine resolutions like zoom level of 9, it is visually hard to observe the cross-sections,
# so we will use a zoom level of 4 for a better visualization
uxda_coarse = ux.UxDataset.from_healpix(model_run(zoom=4, time="P1D").to_dask())["ts"]
uxda_coarse.uxgrid.face_node_connectivity

uxda_lat = uxda_coarse.cross_section.constant_latitude(boulder_lat)
uxda_lat

import geoviews.feature as gf

uxda_lat.isel(time=0).plot(
    rasterize=False,
    projection=projection,
    global_extent=True,
    cmap="inferno",
    clim=(220, 310),
    features=["coastline"],
    title=f"Global surface temperature cross-section at {boulder_lat} degrees latitude",
    width=700,
) * gf.grid(projection=projection);

Let’s also look at the mean of the cross-section:

print(
    f"Mean at {boulder_lat} degrees lat (Boulder, CO, USA): {uxda_lat.mean().values} K"
)

Latitude interval¶

uxda_lat_interval = uxda_coarse.cross_section.constant_latitude_interval(
    [boulder_lat - 15, boulder_lat + 15]
)

uxda_lat_interval.isel(time=0).plot(
    rasterize=False,
    projection=projection,
    global_extent=True,
    cmap="inferno",
    clim=(220, 310),
    features=["coastline"],
    title=f"Global surface temperature cross-section at the latitude interval [{boulder_lat-5},{boulder_lat+5}] degrees",
    width=700,
) * gf.grid(projection=projection);

print(
    f"Mean at the latitude interval of [{boulder_lat-5},{boulder_lat+5}] degrees (-/+15 degrees Boulder, CO, USA): {uxda_lat_interval.mean().values} K"
)

Non-conservative zonal mean¶

Calculating the zonal mean is easy by providing the latitude range between -90 and 90 degrees with a step size in degrees:

%%time
zonal_mean_ts = uxda.isel(time=0).zonal_mean(lat=(-90, 90, 5))

zonal_mean_ts

import hvplot.xarray

(
    uxda.isel(time=0).plot(
        cmap="inferno",
        # periodic_elements="split",
        height=300,
        width=600,
        colorbar=False,
        ylim=(-90, 90),
    )
    + zonal_mean_ts.hvplot.line(
        y="latitudes",
        size=300,
        # width=180,
        # ylabel="",
        ylim=(-90, 90),
        xlim=(220, 310),
        # xticks=[220, 250, 280, 310],
        yticks=[-90, -45, 0, 45, 90],
    )
).opts(title="Temperature and its Zonal means at every 5 degree latitude")

Remapping¶

Now, we will be looking into one of many possible use cases where remapping would be helpful.

The data set we have been using in this section so far belongs to the newer ICON simulation, icon_d3hp003, while there is an older simulation as well, icon_ngc4008, in the same catalog. They are both stored in the HEALPix format in this case, but for most of the model intercomparison workflows in general, they might not be even so. UXarray would still be helpful to remap of those model outputs to other and then make comparisons since it can support several most commonly used unstructrued grid formats.

In this particular case, we still have some use for UXarray’s remapping such that the newer simulation has the zoom = 9 as the maximum available resolution, while the older one has zoom = 10 available. Unfortunately at zoom = 10 though, there is no actual data simulated for ts, the surface temperature. If there was, we could remap the newer simulation’s output to that one, so we could have both of them at zoom = 10, and then we could look into the difference between them for instance. Let’s pretend the highest zoom-level in the newer data is zoom = 8 then, and remap that one into the older simulation’s grid.

Let’s start with opening the older simulation run first:

model_run_older = cat.online.icon_ngc4008

ds_older = model_run_older(zoom=8, time="P1D").to_dask()

uxds_older = ux.UxDataset.from_healpix(ds_older)
uxds_older

Plot the older simulation for reference¶

Let’s have a quick look at how the global surface temperature looks like in the older simulation’s output:

uxds_older["ts"].isel(time=0).plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global surface temperature (zoom = 8) - Older Simulation",
    width=700,
);

Visually there does not seem to be a huge difference between this and the newer simulation’s output we had plotted in the very beginning.

Remap the old simulation output to the newer one¶

Let’s start remapping! For that, we will use the uxgrid of the newer simulation output as the destination grid and use an inverse distance weighted implementation:

%%time
uxda_older_remapped = uxds_older["ts"].isel(time=0).remap.inverse_distance_weighted(
    uxds.uxgrid, k=3, remap_to="face centers", coord_type="cartesian"
)

Plot the difference between the old and newer simulations¶

Now that we have the older and newer model outputs on the same grid, let’s look at the surface temperature differences between the two:

(uxda.isel(time=0) - uxda_older_remapped).plot(
    projection=projection,
    cmap="RdBu_r",
    features=["borders", "coastline"],
    title="Global surface temperature difference between older and newer simulations (zoom = 9)",
    clim=(-25,25),
    width=700,
)