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 uxarray as ux
import xarray as xr
import cartopy.crs as ccrs
import geoviews as gv
import geoviews.feature as gf
import warnings

import holoviews as hv

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 31 ms, sys: 0 ns, total: 31 ms
Wall time: 30.8 ms

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