Constructing a Grid from Points

In many cases, data can be represented as an unstructured series of points, including data from climate models when not paired with any connectivity information or radar observations. UXarray is written around the UGRID conventions, which requires a minimal set of coordinate and connectivity variables to represent a two-dimensional grid. This notebook demonstrates how grid connectivity can be constructured using point data.

import warnings

import cartopy.crs as ccrs
import geoviews as gv
import geoviews.feature as gf
import holoviews as hv
import xarray as xr

import uxarray as ux

hv.extension("matplotlib")

warnings.filterwarnings("ignore")

Types of Point Data

Different types of point data can be encountered depending on the coverage and structure of the data. The following table categorizes these types, providing examples of common use cases.

Domain	Description
Global	Data that covers the entire globe (e.g., atmospheric or climate simulations)
Global with Holes	Data that spans most of the globe but has gaps or regions without observations (e.g., land-only or ocean-only data).
Regional	Data that is limited to a specific area (e.g. local weather forecasts or satellite observations)
Regional with Holes	Regional data with missing sections or observations within the area of interest, often due to obstacles or coverage limitations.

For this notebook, we will be using the coordinates from three testing grids to represent our point data:

• outCSne30.ug: Global Cube Sphere Grid

• oQU480.23010.nc: Global Ocean Grid

• Subset of outCSne30.ug: 9 points centered about (0, 0)

uxgrid_global = ux.open_grid("../../test/meshfiles/ugrid/outCSne30/outCSne30.ug")
uxgrid_global_ocean = ux.open_grid("../../test/meshfiles/mpas/QU/oQU480.231010.nc")
uxgrid_global_ocean.normalize_cartesian_coordinates()
uxgrid_regional = uxgrid_global.subset.nearest_neighbor((0.0, 0.0), k=50)

(
    uxgrid_global.plot.face_centers(
        global_extent=True,
        features=["grid"],
        title="Global Points",
        height=500,
        width=1000,
        s=20,
    )
    + uxgrid_global_ocean.plot.face_centers(
        global_extent=True,
        features=["grid"],
        title="Global Points with Holes",
        height=500,
        width=1000,
        s=20,
    )
    + uxgrid_regional.plot.face_centers(
        global_extent=True,
        features=["grid"],
        title="Regional Points",
        height=500,
        width=1000,
        s=20,
    )
).cols(1).opts(fig_size=300)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/core/dimension.py:1381, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1374 def _repr_mimebundle_(self, include=None, exclude=None):
   1375     """Resolves the class hierarchy for the class rendering the
   1376     object using any display hooks registered on Store.display
   1377     hooks.  The output of all registered display_hooks is then
   1378     combined and returned.
   1379 
   1380     """
-> 1381     return Store.render(self)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/core/options.py:1433, in Store.render(cls, obj)
   1431 data, metadata = {}, {}
   1432 for hook in hooks:
-> 1433     ret = hook(obj)
   1434     if ret is None:
   1435         continue

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/ipython/display_hooks.py:340, in pprint_display(obj)
    338 if not ip.display_formatter.formatters['text/plain'].pprint:
    339     return None
--> 340 return display(obj, raw_output=True)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/ipython/display_hooks.py:311, in display(obj, raw_output, **kwargs)
    309 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    310     with option_state(obj):
--> 311         output = layout_display(obj)
    312 elif isinstance(obj, (HoloMap, DynamicMap)):
    313     with option_state(obj):

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/ipython/display_hooks.py:205, in display_hook.<locals>.wrapped(element)
    203 try:
    204     max_frames = OutputSettings.options['max_frames']
--> 205     mimebundle = fn(element, max_frames=max_frames)
    206     if mimebundle is None:
    207         return {}, {}

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/ipython/display_hooks.py:274, in layout_display(layout, max_frames)
    271     max_frame_warning(max_frames)
    272     return None
--> 274 return render(layout)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/ipython/display_hooks.py:76, in render(obj, **kwargs)
     73 if renderer.fig == 'pdf':
     74     renderer = renderer.instance(fig='png')
---> 76 return renderer.components(obj, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/renderer.py:377, in Renderer.components(self, obj, fmt, comm, **kwargs)
    375     plot = obj
    376 else:
--> 377     plot, fmt = self._validate(obj, fmt)
    379 if not isinstance(plot, Viewable):
    380     html = self._figure_data(plot, fmt, as_script=True, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/renderer.py:308, in Renderer._validate(self, obj, fmt, **kwargs)
    305     plot = HoloViewsPane(obj, center=self.center, backend=self.backend,
    306                          renderer=self)
    307 else:
--> 308     plot = self.get_plot(obj, renderer=self, **kwargs)
    310 all_formats = set(fig_formats + holomap_formats)
    311 if fmt not in all_formats:

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/renderer.py:239, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    236     defaults = [kd.default for kd in plot.dimensions]
    237     init_key = tuple(v if d is None else d for v, d in
    238                      zip(plot.keys[0], defaults, strict=None))
--> 239     plot.update(init_key)
    240 else:
    241     plot = obj

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/mpl/plot.py:293, in MPLPlot.update(self, key)
    291 def update(self, key):
    292     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
--> 293         return self.initialize_plot()
    294     return self.__getitem__(key)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/mpl/plot.py:62, in mpl_rc_context.<locals>.wrapper(self, *args, **kwargs)
     60 def wrapper(self, *args, **kwargs):
     61     with _rc_context(self.fig_rcparams):
---> 62         return f(self, *args, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/mpl/plot.py:1153, in LayoutPlot.initialize_plot(self)
   1148 traverse_fn = lambda x: x.handles.get('bbox_extra_artists', None)
   1149 extra_artists = list(
   1150     chain.from_iterable(artists for artists in self.traverse(traverse_fn)
   1151                         if artists is not None)
   1152 )
-> 1153 fix_aspect(fig, self.rows, self.cols,
   1154            title_obj, extra_artists,
   1155            vspace=self.vspace*self.fig_scale,
   1156            hspace=self.hspace*self.fig_scale)
   1157 colorbars = self.traverse(specs=[lambda x: hasattr(x, 'colorbar')])
   1158 for cbar_plot in colorbars:

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/holoviews/plotting/mpl/util.py:283, in fix_aspect(fig, nrows, ncols, title, extra_artists, vspace, hspace)
    281 if title and title.get_text():
    282     offset = title.get_window_extent().height/fig.dpi
--> 283 fig.set_size_inches(w, (w*aspect)+offset)
    285 # Redraw and adjust title position if defined
    286 fig.canvas.draw()

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/v2026.06.0/lib/python3.14/site-packages/matplotlib/figure.py:3167, in Figure.set_size_inches(self, w, h, forward)
   3165 size = np.array([w, h])
   3166 if not np.isfinite(size).all() or (size < 0).any():
-> 3167     raise ValueError(f'figure size must be positive finite not {size}')
   3168 self.bbox_inches.p1 = size
   3169 if forward:

ValueError: figure size must be positive finite not [12. nan]

:Layout
   .Overlay.I   :Overlay
      .Points.I :Points   [lon,lat]
      .Grid.I   :Feature   [Longitude,Latitude]
   .Overlay.II  :Overlay
      .Points.I :Points   [lon,lat]
      .Grid.I   :Feature   [Longitude,Latitude]
   .Overlay.III :Overlay
      .Points.I :Points   [lon,lat]
      .Grid.I   :Feature   [Longitude,Latitude]

Preparing Point Data

UXarray’s Grid.from_points() method supports both Spherical (lon and lat) and Cartesian (x, y, z) coordinates. It is important to note that the coordinate arrays must be unique in order to run the following methods.

Below we extract the Cartesian (x, y, z) coordinates which we will use for constructing our grid.

x_global, y_global, z_global = (
    uxgrid_global.face_x.values,
    uxgrid_global.face_y.values,
    uxgrid_global.face_z.values,
)
points_global = (x_global, y_global, z_global)

x_global_ocean, y_global_ocean, z_global_ocean = (
    uxgrid_global_ocean.face_x.values,
    uxgrid_global_ocean.face_y.values,
    uxgrid_global_ocean.face_z.values,
)
points_global_ocean = (x_global_ocean, y_global_ocean, z_global_ocean)
boundary_points_global_ocean = uxgrid_global_ocean.boundary_face_indices.values

x_regional, y_regional, z_regional = (
    uxgrid_regional.face_x.values,
    uxgrid_regional.face_y.values,
    uxgrid_regional.face_z.values,
)
points_regional = (x_regional, y_regional, z_regional)

Global Data

The following algorithms will returns grids with a full coverage of the surface of a sphere, which makes them suitable for constructing connectivity from global point data.

Spherical Delaunay

The spherical_delaunay method in the Grid.from_points() function is designed to perform Delaunay triangulation on points distributed over a spherical surface.

How It Works

1. Input Points on the Sphere:

   • The method accepts input points defined in spherical coordinates (e.g., latitude and longitude) or Cartesian coordinates (x, y, z) that lie on the surface of the sphere. They are internally converted to normalized Cartesian coordinates.

2. Computing the Convex Hull:

   • The algorithm computes the Convex Hull of the set of Cartesian points. The convex hull is the smallest convex shape that encloses all the points. In three dimensions, the convex hull forms a polyhedron where each face is a triangle.

3. Extracting Triangles:

   • Once the convex hull is determined, the triangular faces of the hull are extracted. These triangles represent the Delaunay triangulation on the sphere’s surface, ensuring that no point is inside the circumcircle of any triangle, which is a key property of Delaunay triangulations.

%%time
grid_dt = ux.Grid.from_points(points_global, method="spherical_delaunay")

grid_dt.plot(
    projection=ccrs.Robinson(),
    linewidth=0.5,
    periodic_elements="split",
    title="Spherical Delaunay Triangulation",
    height=500,
    width=1000,
)

The resulting grid will always be strictly triangular and cover the entire sphere.

grid_dt.triangular

grid_dt.plot.face_degree_distribution()

grid_dt.plot.face_area_distribution(bins=10)

grid_dt.global_sphere_coverage

Zooming in, we can observe the Delaunay Triangles in detail. The original point coordinates are now the corners of our faces. This means that any data that was originally mapped to the points will reside on the corner nodes.

(grid_dt.plot() * uxgrid_global.plot.face_centers(color="red", s=1000)).opts(
    xlim=(-10, 10),
    ylim=(-5, 5),
    title="Spherical Delaunay Triangles (Zoomed)",
)

Spherical Voronoi

The spherical_voronoi method in the Grid.from_points() function is designed to generate a Voronoi tessellation on points distributed over a spherical surface. This method leverages SciPy’s Spherical Voronoi functionality internally.

How It Works

1. Input Points on the Sphere:

   • The method accepts input points defined in spherical coordinates (e.g., latitude and longitude) or Cartesian coordinates (x, y, z) that lie on the surface of the sphere. They are internally converted to normalized Cartesian coordinates.

2. Computing the Spherical Voronoi Diagram:

   • Using SciPy’s SphericalVoronoi class, the algorithm computes the Voronoi tessellation on the sphere. This involves determining the regions on the sphere where each region contains all points closer to one generating point than to any other.

3. Constructing Voronoi Regions:

   • The Spherical Voronoi algorithm identifies the vertices and edges that define each Voronoi region. Each region corresponds to one input point and consists of a polygon on the sphere’s surface.

%%time
grid_sv = ux.Grid.from_points(points_global, method="spherical_voronoi")

grid_sv.plot(
    projection=ccrs.Robinson(),
    linewidth=0.5,
    periodic_elements="split",
    height=500,
    width=1000,
    title="Spherical Voronoi Tesselation",
)

The resulting grid consists of mostly 6-sided faces, with small numbers of faces with 4, 5, 7, and 8 sides.

grid_sv.plot.face_degree_distribution()

grid_sv.plot.face_area_distribution(bins=15)

grid_sv.global_sphere_coverage

Zooming in, we can observe the Voronoi Regions in detail. The original point coordinates are now the centers of each the faces in the grid. This means that any data that was originally mapped to the points now resides on the faces.

(grid_sv.plot() * uxgrid_global.plot.face_centers(color="red")).opts(
    xlim=(-10, 10),
    ylim=(-5, 5),
    title="Spherical Voronoi Cells (Zoomed)",
)

(grid_sv.plot() * uxgrid_global.plot.face_centers(color="red")).opts(
    xlim=(14.5, 18.5),
    ylim=(5.5, 9.0),
    title="Single Spherical Voronoi Cell (Zoomed)",
)

Limitations of Spherical Methods

The spherical methods discussed above are not appropriate for regional data, as the exterior boundaries of the region will wrap around and connect together, forming extremely large faces.

%%time
grid_dt_regional = ux.Grid.from_points(points_regional, method="spherical_delaunay")

grid_dt_regional.plot.face_area_distribution(
    bins=15, title="Delaunay: Face Area Distributon (Regional)"
)

grid_sv_regional = ux.Grid.from_points(points_regional, method="spherical_voronoi")

grid_sv_regional.plot.face_area_distribution(
    bins=15, title="Voronoi: Face Area Distributon (Regional)"
)

Global Data with Holes

For global point data with holes, the spherical methods can be used, but there are certain considerations that need to be made, since by default, each spherical method returns a grid with full sphere coverage.

Spherical Delaunay

%%time
grid_dt = ux.Grid.from_points(points_global_ocean, method="spherical_delaunay")

grid_dt.global_sphere_coverage

grid_dt.plot(
    projection=ccrs.Robinson(),
    linewidth=0.5,
    periodic_elements="exclude",
    title="Spherical Delaunay Triangulation",
    height=500,
    width=1000,
)

(
    grid_dt.plot(
        features=["coastline"],
    )
    * uxgrid_global_ocean.plot.face_centers(color="red")
).opts(
    xlim=(-20, 20),
    ylim=(-10, 10),
    title="Spherical Delaunay Triangles (Zoomed)",
)

This behavior is not always desired, especially if you do not want elements over previously empty regions. The Grid.from_points() method accepts an optional argument boundary_points, which is an array of indices corresponding to which points lie on a defined boundary.

%%time
grid_dt_no_boundary = ux.Grid.from_points(
    points_global_ocean,
    boundary_points=boundary_points_global_ocean,
    method="spherical_delaunay",
)

When appropriate boundary points are provided, the resulting grid has a partial sphere coverage.

grid_dt_no_boundary.global_sphere_coverage

grid_dt_no_boundary.plot(
    projection=ccrs.Robinson(),
    linewidth=0.5,
    periodic_elements="exclude",
    height=500,
    width=1000,
    title="Spherical Delaunay Triangulation without Boundary Points",
)

(
    grid_dt_no_boundary.plot(
        features=["coastline"],
    )
    * uxgrid_global_ocean.plot.face_centers(color="red")
).opts(
    xlim=(-20, 20),
    ylim=(-10, 10),
    title="Spherical Delaunay Triangles without Boundary Points (Zoomed)",
)

Spherical Voronoi

The Spherical Voronoi method can be used for global poitns with holes, however it does not support a boundary_points parameter, meaning that the resulting Grid will always have a global sphere coverage.

%%time
grid_sv = ux.Grid.from_points(points_global_ocean, method="spherical_voronoi")

grid_sv.global_sphere_coverage

grid_sv.plot(
    projection=ccrs.Robinson(),
    linewidth=0.5,
    periodic_elements="exclude",
    height=500,
    width=1000,
    title="Spherical Voronoi Regions",
)

(
    grid_sv.plot(
        features=["coastline"],
    )
    * uxgrid_global_ocean.plot.face_centers(color="red")
).opts(
    xlim=(-20, 20),
    ylim=(-10, 10),
    title="Spherical Voronoi Regions (Zoomed)",
)

Regional Data

The regional delaunay method can be used to construct a grid from points in a regional area.

How It Works

1. Input Points on the Sphere:

   • The method accepts input points defined in spherical coordinates (e.g., latitude and longitude) or Cartesian coordinates (x, y, z) that lie on the surface of the sphere. They are internally converted to normalized Cartesian coordinates.

2. Computing the Regional Delaunay Diagram:

   • The method projects the points to a 2D plane using stereographic projection, followed by SciPy’s Delaunay triangulation method to construct the grid.

3. Extracting Triangles:

   • Once the triangles of 2D points are determined, the connectivity of the triangular faces are extracted. These triangles represent the Delaunay triangulation on the sphere’s surface, ensuring that no point is inside the circumcircle of any triangle, which is a key property of Delaunay triangulations.

grid_r = ux.Grid.from_points(points_regional, method="regional_delaunay")

grid_r.plot(
    linewidth=0.5,
    periodic_elements="exclude",
    height=500,
    width=1000,
    title="Regional Delaunay Regions",
)

Antimerdian

This also works on regions wrapping the antimerdian

antimerdian_region = uxgrid_global.subset.bounding_circle((-180, 0), 10)

x_antimerdian_region, y_antimerdian_region, z_antimerdian_region = (
    antimerdian_region.face_x.values,
    antimerdian_region.face_y.values,
    antimerdian_region.face_z.values,
)
antimerdian_region = (x_antimerdian_region, y_antimerdian_region, z_antimerdian_region)

grid_r = ux.Grid.from_points(antimerdian_region, method="regional_delaunay")

grid_r.plot(
    linewidth=0.5,
    periodic_elements="exclude",
    height=500,
    width=1000,
    title="Regional Delaunay Regions",
)

Creating a Regional Unstructured Grid from Points

UXarray allows users to create unstructured grids from scattered (lon, lat) point coordinates using Delaunay triangulation. When constructing regional unstructured grids with the method=”regional_delaunay” option, it is critical to explicitly specify boundary points to avoid mesh artifacts and ensure accurate face bounds.

import numpy as np

import uxarray as ux

# Define a regular grid of longitude and latitude points over [0, 60] x [0, 60]
lon, lat = np.meshgrid(
    np.linspace(0, 60.0, 10, dtype=np.float32),
    np.linspace(0, 60.0, 10, dtype=np.float32),
)
lon_flat = lon.ravel()
lat_flat = lat.ravel()

# Identify points along the domain boundary
mask = (
    np.isclose(lon_flat, 0.0)
    | np.isclose(lon_flat, 60.0)
    | np.isclose(lat_flat, 0.0)
    | np.isclose(lat_flat, 60.0)
)
boundary_points = np.flatnonzero(mask)

# Create the unstructured grid using the regional Delaunay method
uxgrid = ux.Grid.from_points(
    (lon_flat, lat_flat), method="regional_delaunay", boundary_points=boundary_points
)

Why Specify boundary_points?

The internal triangulation algorithm assumes a spherical domain. Without user-defined constraints, it may attempt to “wrap” the grid edges to form a closed surface, which is inappropriate for bounded regional domains. This can cause:

• Element bounds that exceed the expected coordinate extents.

• Spatial hash errors in downstream workflows.

• Unexpected overlaps or distortions in edge geometry.

By supplying boundary_points, you ensure:

• Proper triangulation only within the region of interest.

• Accurate face_bounds_lat and face_bounds_lon.

• Stable behavior for spatial indexing and remapping.