Constructing a Grid from Points

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)

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")

CPU times: user 22.2 ms, sys: 26 μs, total: 22.2 ms
Wall time: 22.1 ms

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

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 1
----> 1 grid_dt.plot(
      2     projection=ccrs.Robinson(),
      3     linewidth=0.5,
      4     periodic_elements="split",
      5     title="Spherical Delaunay Triangulation",
      6     height=500,
      7     width=1000,
      8 )

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/plot/accessor.py:42, in GridPlotAccessor.__call__(self, **kwargs)
     41 def __call__(self, **kwargs) -> Any:
---> 42     return self.edges(**kwargs)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/plot/accessor.py:232, in GridPlotAccessor.edges(self, periodic_elements, backend, engine, **kwargs)
    229         central_longitude = 0.0
    230     kwargs["crs"] = ccrs.PlateCarree(central_longitude=central_longitude)
--> 232 gdf = self._uxgrid.to_geodataframe(
    233     periodic_elements=periodic_elements,
    234     projection=kwargs.get("projection"),
    235     engine=engine,
    236     project=False,
    237 )
    239 return gdf.hvplot.paths(geo=True, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/grid/grid.py:2278, in Grid.to_geodataframe(self, periodic_elements, projection, cache, override, engine, exclude_antimeridian, return_non_nan_polygon_indices, exclude_nan_polygons, **kwargs)
   2275         return self._gdf_cached_parameters["gdf"]
   2277 # construct a GeoDataFrame with the faces stored as polygons as the geometry
-> 2278 gdf, non_nan_polygon_indices = _grid_to_polygon_geodataframe(
   2279     self, periodic_elements, projection, project, engine
   2280 )
   2282 if exclude_nan_polygons and non_nan_polygon_indices is not None:
   2283     # exclude any polygons that contain NaN values
   2284     gdf = GeoDataFrame({"geometry": gdf["geometry"][non_nan_polygon_indices]})

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/grid/geometry.py:224, in _grid_to_polygon_geodataframe(grid, periodic_elements, projection, project, engine)
    221 grid._gdf_cached_parameters["antimeridian_face_indices"] = antimeridian_face_indices
    223 if periodic_elements == "split":
--> 224     gdf = _build_geodataframe_with_antimeridian(
    225         polygon_shells,
    226         projected_polygon_shells,
    227         antimeridian_face_indices,
    228         engine=engine,
    229     )
    230 elif periodic_elements == "ignore":
    231     if engine == "geopandas":
    232         # create a geopandas.GeoDataFrame

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/grid/geometry.py:303, in _build_geodataframe_with_antimeridian(polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine)
    300 import spatialpandas
    301 from spatialpandas.geometry import MultiPolygonArray
--> 303 polygons = _build_corrected_shapely_polygons(
    304     polygon_shells, projected_polygon_shells, antimeridian_face_indices
    305 )
    306 if engine == "geopandas":
    307     # Create a geopandas.GeoDataFrame
    308     gdf = geopandas.GeoDataFrame({"geometry": polygons})

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/uxarray/grid/geometry.py:338, in _build_corrected_shapely_polygons(polygon_shells, projected_polygon_shells, antimeridian_face_indices)
    335 antimeridian_polygons = Polygons(polygon_shells[antimeridian_face_indices])
    337 # correct each antimeridian polygon
--> 338 corrected_polygons = [antimeridian.fix_polygon(P) for P in antimeridian_polygons]
    340 # insert correct polygon back into original array
    341 for i in reversed(antimeridian_face_indices):

File ~/checkouts/readthedocs.org/user_builds/uxarray/conda/1454/lib/python3.14/site-packages/antimeridian/_implementation.py:368, in fix_polygon(polygon, force_north_pole, force_south_pole, fix_winding, great_circle)
    366             return pole_covering_polygon
    367         else:
--> 368             raise ValueError(
    369                 "Fixed polygon is invalid, check your input polygon for validity. "
    370                 "Reason your polygon is invalid: "
    371                 + shapely.is_valid_reason(polygon)
    372             )
    374 else:
    375     return MultiPolygon(polygons)

ValueError: Fixed polygon is invalid, check your input polygon for validity. Reason your polygon is invalid: Self-intersection[-45 87.8791656494141]

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.

Contents