UGRID Unstructured Mesh Exploration¶

This notebook demonstrates pyramids' full UGRID support using the UGRID convention NC estuary mesh — a D-Flow FM (Delft3D Flexible Mesh) output covering the UGRID convention NC river estuary in the Netherlands.

We'll cover:

1. Reading a UGRID NetCDF file
2. Exploring mesh topology (nodes, faces, edges)
3. Inspecting metadata and global attributes
4. Computing geometric properties (centroids, areas)
5. Creating synthetic bathymetry & velocity data
6. Plotting mesh data (face & node coloring, wireframe)
7. Spatial operations (clipping, bounding box subsetting)
8. Mesh-to-raster conversion (bridging to regular grids)
9. GeoDataFrame / FeatureCollection interop
10. CRS reprojection
11. Writing UGRID-compliant NetCDF
12. Time dimension support & animation

In [1]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

from pyramids.netcdf.ugrid import (
    Connectivity,
    Mesh2d,
    MeshSpatialIndex,
    MeshTopologyInfo,
    MeshVariable,
    UgridDataset,
)

from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.netcdf.ugrid import ( Connectivity, Mesh2d, MeshSpatialIndex, MeshTopologyInfo, MeshVariable, UgridDataset, )

2026-04-11 21:16:01 | INFO | pyramids.base.config | Logging is configured.

1. Reading a UGRID NetCDF File¶

UgridDataset.read_file() automatically detects the mesh topology, separates data variables from coordinate/connectivity arrays, and builds the Mesh2d object.

In [2]:

ugrid_path = Path("../../../../tests/data/netcdf/ugrid/ugrid.nc")
ds = UgridDataset.read_file(ugrid_path)
print(ds)

ugrid_path = Path("../../../../tests/data/netcdf/ugrid/ugrid.nc") ds = UgridDataset.read_file(ugrid_path) print(ds)

UgridDataset: ..\..\tests\mo\netcdf\ugrid_convention_nc01_waqgeom.nc
  Mesh: mesh2d
  Nodes: 8916, Faces: 8355, Edges: 17270
  Bounds: (17325.978515625, 359012.23738230974, 86177.27534049794, 392033.78125)
  CRS: unknown
  Data variables (2):
    mesh2d_node_z: location=node, shape=(8916,)
    mesh2d_edge_type: location=edge, shape=(17270,)

2. Exploring Mesh Topology¶

The mesh consists of nodes (vertices), faces (polygons), and edges (connections between nodes).

In [3]:

mesh = ds.mesh

print(f"Nodes:  {mesh.n_node:,}")
print(f"Faces:  {mesh.n_face:,}")
print(f"Edges:  {mesh.n_edge:,}")
print(f"Bounds: {mesh.bounds}")
print()

# Face-node connectivity: which nodes define each face
fnc = mesh.face_node_connectivity
print(f"Face-node connectivity shape: {fnc.data.shape}")
print(f"Max nodes per face: {fnc.max_nodes_per_element}")
print(f"All triangular? {fnc.is_triangular()}")
print(f"Nodes per face (first 10): {fnc.nodes_per_element()[:10]}")
print()

# Inspect a single face
face_idx = 100
nodes = mesh.get_face_nodes(face_idx)
polygon = mesh.get_face_polygon(face_idx)
print(f"Face {face_idx} nodes: {nodes}")
print(f"Face {face_idx} polygon coordinates:\n{polygon}")

mesh = ds.mesh print(f"Nodes: {mesh.n_node:,}") print(f"Faces: {mesh.n_face:,}") print(f"Edges: {mesh.n_edge:,}") print(f"Bounds: {mesh.bounds}") print() # Face-node connectivity: which nodes define each face fnc = mesh.face_node_connectivity print(f"Face-node connectivity shape: {fnc.data.shape}") print(f"Max nodes per face: {fnc.max_nodes_per_element}") print(f"All triangular? {fnc.is_triangular()}") print(f"Nodes per face (first 10): {fnc.nodes_per_element()[:10]}") print() # Inspect a single face face_idx = 100 nodes = mesh.get_face_nodes(face_idx) polygon = mesh.get_face_polygon(face_idx) print(f"Face {face_idx} nodes: {nodes}") print(f"Face {face_idx} polygon coordinates:\n{polygon}")

Nodes:  8,916
Faces:  8,355
Edges:  17,270
Bounds: (17325.978515625, 359012.23738230974, 86177.27534049794, 392033.78125)

Face-node connectivity shape: (8355, 4)
Max nodes per face: 4
All triangular? False
Nodes per face (first 10): [3 3 3 3 3 3 3 3 3 3]

Face 100 nodes: [1175 1184 8915]
Face 100 polygon coordinates:
[[ 73692.03640027 365995.09281191]
 [ 73573.99345526 366073.37504154]
 [ 73570.51812025 365981.23808658]]

In [4]:

# Edge connectivity
enc = mesh.edge_node_connectivity
print(f"Edge-node connectivity shape: {enc.data.shape}")
print(f"Edge 0 connects nodes: {enc.get_element(0)}")

start, end = mesh.get_edge_coords(0)
print(f"Edge 0: ({start[0]:.1f}, {start[1]:.1f}) -> ({end[0]:.1f}, {end[1]:.1f})")

# Edge connectivity enc = mesh.edge_node_connectivity print(f"Edge-node connectivity shape: {enc.data.shape}") print(f"Edge 0 connects nodes: {enc.get_element(0)}") start, end = mesh.get_edge_coords(0) print(f"Edge 0: ({start[0]:.1f}, {start[1]:.1f}) -> ({end[0]:.1f}, {end[1]:.1f})")

Edge-node connectivity shape: (17270, 2)
Edge 0 connects nodes: [   1 8907]
Edge 0: (77140.1, 369348.9) -> (77041.0, 369670.8)

3. Metadata and Global Attributes¶

In [5]:

# Global attributes
for key, value in ds.global_attributes.items():
    print(f"  {key}: {value}")

print(f"\nMesh name: {ds.mesh_name}")
print(f"Data variables: {ds.data_variable_names}")
print(f"CRS/EPSG: {ds.epsg}")

# Global attributes for key, value in ds.global_attributes.items(): print(f" {key}: {value}") print(f"\nMesh name: {ds.mesh_name}") print(f"Data variables: {ds.data_variable_names}") print(f"CRS/EPSG: {ds.epsg}")

  institution: Deltares
  references: http://www.deltares.nl
  source: D-Flow FM 1.2.63.64757M. Model:
  history: Created on 2019-09-18T15:20:03+0200, D-Flow FM
  date_created: 2019-09-18T15:20:03+0200
  date_modified: 2019-09-18T15:20:03+0200
  Conventions: CF-1.8 UGRID-1.0 Deltares-0.10

Mesh name: mesh2d
Data variables: ['mesh2d_node_z', 'mesh2d_edge_type']
CRS/EPSG: None

In [6]:

# Inspect data variables
for name in ds.data_variable_names:
    var = ds[name]
    print(f"{name}:")
    print(f"  location:      {var.location}")
    print(f"  shape:         {var.shape}")
    print(f"  dtype:         {var.dtype}")
    print(f"  standard_name: {var.standard_name}")
    print(f"  has_time:      {var.has_time}")
    print(f"  n_elements:    {var.n_elements}")
    print()

# Inspect data variables for name in ds.data_variable_names: var = ds[name] print(f"{name}:") print(f" location: {var.location}") print(f" shape: {var.shape}") print(f" dtype: {var.dtype}") print(f" standard_name: {var.standard_name}") print(f" has_time: {var.has_time}") print(f" n_elements: {var.n_elements}") print()

mesh2d_node_z:
  location:      node
  shape:         (8916,)
  dtype:         float64
  standard_name: altitude
  has_time:      False
  n_elements:    8916

mesh2d_edge_type:
  location:      edge
  shape:         (17270,)
  dtype:         int32
  standard_name: None
  has_time:      False
  n_elements:    17270

In [7]:

# Full metadata summary
meta = ds.metadata
print(f"Conventions: {meta.conventions}")
print(f"Topologies:  {len(meta.mesh_topologies)}")
print(f"Variables:   {meta.data_variables}")

# Full metadata summary meta = ds.metadata print(f"Conventions: {meta.conventions}") print(f"Topologies: {len(meta.mesh_topologies)}") print(f"Variables: {meta.data_variables}")

Conventions: CF-1.8 UGRID-1.0 Deltares-0.10
Topologies:  1
Variables:   {'mesh2d_node_z': 'node', 'mesh2d_edge_type': 'edge'}

4. Geometric Properties¶

Face centroids, areas, and fan triangulation are computed lazily and cached on first access.

In [8]:

# Face centroids
cx, cy = mesh.face_centroids
print(
    f"Centroid range: x=[{cx.min():.0f}, {cx.max():.0f}], y=[{cy.min():.0f}, {cy.max():.0f}]"
)

# Face areas (shoelace formula, vectorized)
areas = mesh.face_areas
print(
    f"Face areas: min={areas.min():.0f} m², max={areas.max():.0f} m², mean={areas.mean():.0f} m²"
)

# Fan triangulation (pure numpy, no matplotlib needed)
tri = mesh.fan_triangles
print(f"Fan triangles: {tri.shape[0]} triangles from {mesh.n_face} faces")

# Face centroids cx, cy = mesh.face_centroids print( f"Centroid range: x=[{cx.min():.0f}, {cx.max():.0f}], y=[{cy.min():.0f}, {cy.max():.0f}]" ) # Face areas (shoelace formula, vectorized) areas = mesh.face_areas print( f"Face areas: min={areas.min():.0f} m², max={areas.max():.0f} m², mean={areas.mean():.0f} m²" ) # Fan triangulation (pure numpy, no matplotlib needed) tri = mesh.fan_triangles print(f"Fan triangles: {tri.shape[0]} triangles from {mesh.n_face} faces")

Centroid range: x=[17519, 86141], y=[359201, 391746]
Face areas: min=220 m², max=262424 m², mean=52640 m²
Fan triangles: 16583 triangles from 8355 faces

In [9]:

# Plot face area distribution
fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(areas, bins=50, edgecolor="black", alpha=0.7)
ax.set_xlabel("Face Area (m²)")
ax.set_ylabel("Count")
ax.set_title("Distribution of Mesh Face Areas")
plt.tight_layout()
plt.show()

# Plot face area distribution fig, ax = plt.subplots(figsize=(8, 4)) ax.hist(areas, bins=50, edgecolor="black", alpha=0.7) ax.set_xlabel("Face Area (m²)") ax.set_ylabel("Count") ax.set_title("Distribution of Mesh Face Areas") plt.tight_layout() plt.show()

5. Creating Synthetic Bathymetry & Velocity Data¶

The UGRID convention NC file is geometry-only (all node_z = -999). Let's create realistic synthetic data to demonstrate plotting.

In [10]:

# Create synthetic bathymetry (depth) on faces
# Deeper in the center of the estuary, shallower near edges
xmin, ymin, xmax, ymax = mesh.bounds
x_norm = (cx - xmin) / (xmax - xmin)
y_norm = (cy - ymin) / (ymax - ymin)

# Channel center line roughly follows the estuary
channel_dist = np.abs(y_norm - 0.5)  # distance from center
bathymetry = -15.0 + 12.0 * channel_dist + 3.0 * np.sin(x_norm * 4 * np.pi)
bathymetry += np.random.RandomState(42).normal(0, 0.5, mesh.n_face)

print(f"Bathymetry range: {bathymetry.min():.1f} to {bathymetry.max():.1f} m")

# Create synthetic velocity on nodes
node_x_norm = (mesh.node_x - xmin) / (xmax - xmin)
node_y_norm = (mesh.node_y - ymin) / (ymax - ymin)
velocity = (
    1.5
    * (1 - 2 * np.abs(node_y_norm - 0.5))
    * (0.5 + 0.5 * np.sin(node_x_norm * 2 * np.pi))
)
velocity = np.clip(velocity, 0, None)

print(f"Velocity range: {velocity.min():.2f} to {velocity.max():.2f} m/s")

# Create synthetic bathymetry (depth) on faces # Deeper in the center of the estuary, shallower near edges xmin, ymin, xmax, ymax = mesh.bounds x_norm = (cx - xmin) / (xmax - xmin) y_norm = (cy - ymin) / (ymax - ymin) # Channel center line roughly follows the estuary channel_dist = np.abs(y_norm - 0.5) # distance from center bathymetry = -15.0 + 12.0 * channel_dist + 3.0 * np.sin(x_norm * 4 * np.pi) bathymetry += np.random.RandomState(42).normal(0, 0.5, mesh.n_face) print(f"Bathymetry range: {bathymetry.min():.1f} to {bathymetry.max():.1f} m") # Create synthetic velocity on nodes node_x_norm = (mesh.node_x - xmin) / (xmax - xmin) node_y_norm = (mesh.node_y - ymin) / (ymax - ymin) velocity = ( 1.5 * (1 - 2 * np.abs(node_y_norm - 0.5)) * (0.5 + 0.5 * np.sin(node_x_norm * 2 * np.pi)) ) velocity = np.clip(velocity, 0, None) print(f"Velocity range: {velocity.min():.2f} to {velocity.max():.2f} m/s")

Bathymetry range: -19.3 to -5.8 m
Velocity range: 0.00 to 1.45 m/s

In [11]:

# Build a UgridDataset with the synthetic data
synth_ds = UgridDataset.create_from_arrays(
    node_x=mesh.node_x,
    node_y=mesh.node_y,
    face_node_connectivity=mesh.face_node_connectivity.data,
    data={
        "bathymetry": bathymetry,
        "velocity": velocity,
    },
    data_locations={
        "bathymetry": "face",
        "velocity": "node",
    },
    epsg=28992,  # Rijksdriehoekstelsel (Dutch national grid)
    mesh_name="mesh2d",
)
print(synth_ds)

# Build a UgridDataset with the synthetic data synth_ds = UgridDataset.create_from_arrays( node_x=mesh.node_x, node_y=mesh.node_y, face_node_connectivity=mesh.face_node_connectivity.data, data={ "bathymetry": bathymetry, "velocity": velocity, }, data_locations={ "bathymetry": "face", "velocity": "node", }, epsg=28992, # Rijksdriehoekstelsel (Dutch national grid) mesh_name="mesh2d", ) print(synth_ds)

UgridDataset: (in-memory)
  Mesh: mesh2d
  Nodes: 8916, Faces: 8355, Edges: 0
  Bounds: (17325.978515625, 359012.23738230974, 86177.27534049794, 392033.78125)
  CRS: 28992
  Data variables (2):
    bathymetry: location=face, shape=(8355,)
    velocity: location=node, shape=(8916,)

6. Plotting Mesh Data¶

Plotting uses cleopatra.MeshGlyph under the hood. Face data uses tripcolor, node data uses tricontourf.

In [12]:

# Plot face-centered bathymetry
glyph = synth_ds.plot(
    "bathymetry",
    cmap="terrain",
    title="UGRID convention NC Bathymetry (synthetic)",
)
plt.show()

# Plot face-centered bathymetry glyph = synth_ds.plot( "bathymetry", cmap="terrain", title="UGRID convention NC Bathymetry (synthetic)", ) plt.show()

Matplotlib backend set to inline for static plots in Jupyter notebook

In [13]:

# Plot node-centered velocity with smooth contours
glyph = synth_ds.plot(
    "velocity",
    cmap="YlOrRd",
    title="Flow Velocity (synthetic)",
)
plt.show()

# Plot node-centered velocity with smooth contours glyph = synth_ds.plot( "velocity", cmap="YlOrRd", title="Flow Velocity (synthetic)", ) plt.show()

In [14]:

# Wireframe overlay on bathymetry
glyph = synth_ds.plot(
    "bathymetry",
    cmap="ocean_r",
    title="Bathymetry with Mesh Wireframe",
    edgecolor="gray",
)
plt.show()

# Wireframe overlay on bathymetry glyph = synth_ds.plot( "bathymetry", cmap="ocean_r", title="Bathymetry with Mesh Wireframe", edgecolor="gray", ) plt.show()

In [15]:

# Wireframe only
glyph = synth_ds.plot_outline(color="steelblue", linewidth=0.2)
plt.title("UGRID convention NC Mesh Structure")
plt.show()

# Wireframe only glyph = synth_ds.plot_outline(color="steelblue", linewidth=0.2) plt.title("UGRID convention NC Mesh Structure") plt.show()

7. Spatial Operations¶

Spatial Indexing¶

In [16]:

# Create a spatial index for fast queries
idx = MeshSpatialIndex(synth_ds.mesh)

# Find nearest face to a point (mid-estuary)
mid_x = (xmin + xmax) / 2
mid_y = (ymin + ymax) / 2
nearest = idx.locate_nearest_face(mid_x, mid_y, k=3)
print(f"3 nearest faces to ({mid_x:.0f}, {mid_y:.0f}): {nearest.flatten()}")

# Point-in-face query
face_ids = idx.locate_faces(
    np.array([mid_x, xmin + 1000]),
    np.array([mid_y, ymin + 1000]),
)
print(f"Point-in-face results: {face_ids}")

# Bounding box query
faces_in_box = idx.locate_faces_in_bounds(40000, 375000, 60000, 385000)
print(f"Faces in box [40k-60k, 375k-385k]: {len(faces_in_box)}")

# Create a spatial index for fast queries idx = MeshSpatialIndex(synth_ds.mesh) # Find nearest face to a point (mid-estuary) mid_x = (xmin + xmax) / 2 mid_y = (ymin + ymax) / 2 nearest = idx.locate_nearest_face(mid_x, mid_y, k=3) print(f"3 nearest faces to ({mid_x:.0f}, {mid_y:.0f}): {nearest.flatten()}") # Point-in-face query face_ids = idx.locate_faces( np.array([mid_x, xmin + 1000]), np.array([mid_y, ymin + 1000]), ) print(f"Point-in-face results: {face_ids}") # Bounding box query faces_in_box = idx.locate_faces_in_bounds(40000, 375000, 60000, 385000) print(f"Faces in box [40k-60k, 375k-385k]: {len(faces_in_box)}")

3 nearest faces to (51752, 375523): [6015 5852 6178]
Point-in-face results: [6015   -1]
Faces in box [40k-60k, 375k-385k]: 2273

Bounding Box Subsetting¶

In [17]:

# Subset to the central part of the estuary
subset = synth_ds.subset_by_bounds(35000, 370000, 65000, 390000)
print(f"Original: {synth_ds.n_face} faces, {synth_ds.n_node} nodes")
print(f"Subset:   {subset.n_face} faces, {subset.n_node} nodes")

glyph = subset.plot(
    "bathymetry",
    cmap="terrain",
    title="Central Estuary Subset",
)
plt.show()

# Subset to the central part of the estuary subset = synth_ds.subset_by_bounds(35000, 370000, 65000, 390000) print(f"Original: {synth_ds.n_face} faces, {synth_ds.n_node} nodes") print(f"Subset: {subset.n_face} faces, {subset.n_node} nodes") glyph = subset.plot( "bathymetry", cmap="terrain", title="Central Estuary Subset", ) plt.show()

Original: 8355 faces, 8916 nodes
Subset:   3729 faces, 3987 nodes

Polygon Clipping¶

In [18]:

from shapely.geometry import Polygon, box

# Clip to a custom polygon (diamond shape in the estuary)
clip_poly = Polygon(
    [
        (45000, 377000),
        (55000, 383000),
        (65000, 377000),
        (55000, 371000),
    ]
)

clipped = synth_ds.clip(clip_poly, touch=True)
print(f"Clipped: {clipped.n_face} faces, {clipped.n_node} nodes")

glyph = clipped.plot(
    "bathymetry",
    cmap="terrain",
    title="Diamond-Clipped Region",
)
plt.show()

from shapely.geometry import Polygon, box # Clip to a custom polygon (diamond shape in the estuary) clip_poly = Polygon( [ (45000, 377000), (55000, 383000), (65000, 377000), (55000, 371000), ] ) clipped = synth_ds.clip(clip_poly, touch=True) print(f"Clipped: {clipped.n_face} faces, {clipped.n_node} nodes") glyph = clipped.plot( "bathymetry", cmap="terrain", title="Diamond-Clipped Region", ) plt.show()

Clipped: 1445 faces, 1602 nodes

8. Mesh-to-Raster Conversion¶

to_dataset() converts mesh data to a regular-grid pyramids Dataset, bridging unstructured and structured worlds.

In [19]:

# Convert bathymetry to a 200m resolution raster
raster = synth_ds.to_dataset(
    "bathymetry",
    cell_size=200.0,
    method="nearest",
)
print(f"Raster: {raster.rows} x {raster.columns}, cell_size={raster.cell_size}")
print(f"Raster bounds: {raster.bounds}")

# Plot the raster
raster.plot(cmap="terrain", title="Bathymetry Raster (200m, nearest)")
plt.show()

# Convert bathymetry to a 200m resolution raster raster = synth_ds.to_dataset( "bathymetry", cell_size=200.0, method="nearest", ) print(f"Raster: {raster.rows} x {raster.columns}, cell_size={raster.cell_size}") print(f"Raster bounds: {raster.bounds}") # Plot the raster raster.plot(cmap="terrain", title="Bathymetry Raster (200m, nearest)") plt.show()

Raster: 166 x 345, cell_size=200.0
Raster bounds:                                             geometry
0  POLYGON ((17325.979 392033.781, 17325.979 3588...

In [20]:

# Linear interpolation for smoother result
raster_lin = synth_ds.to_dataset(
    "velocity",
    cell_size=200.0,
    method="linear",
)
raster_lin.plot(cmap="YlOrRd", title="Velocity Raster (200m, linear)")
plt.show()

# Linear interpolation for smoother result raster_lin = synth_ds.to_dataset( "velocity", cell_size=200.0, method="linear", ) raster_lin.plot(cmap="YlOrRd", title="Velocity Raster (200m, linear)") plt.show()

In [21]:

# Save raster to GeoTIFF
raster.to_file("bathymetry_raster.tif")
print("Saved bathymetry_raster.tif")

# Save raster to GeoTIFF raster.to_file("bathymetry_raster.tif") print("Saved bathymetry_raster.tif")

Saved bathymetry_raster.tif

9. GeoDataFrame / FeatureCollection Interop¶

Convert mesh elements to vector geometries for GIS analysis.

In [22]:

# Faces as Polygons with bathymetry column
gdf_faces = subset.to_geodataframe("bathymetry", location="face")
print(f"GeoDataFrame: {len(gdf_faces)} rows")
print(f"Columns: {list(gdf_faces.columns)}")
print(f"Geometry type: {gdf_faces.geometry.iloc[0].geom_type}")
gdf_faces.head()

# Faces as Polygons with bathymetry column gdf_faces = subset.to_geodataframe("bathymetry", location="face") print(f"GeoDataFrame: {len(gdf_faces)} rows") print(f"Columns: {list(gdf_faces.columns)}") print(f"Geometry type: {gdf_faces.geometry.iloc[0].geom_type}") gdf_faces.head()

GeoDataFrame: 3729 rows
Columns: ['bathymetry', 'geometry']
Geometry type: Polygon

Out[22]:

	bathymetry	geometry
0	-12.432332	POLYGON ((38092.773 389066.469, 37824.484 3889...
1	-12.481534	POLYGON ((38331.969 388645.5, 38058.547 388489...
2	-12.153503	POLYGON ((38622.91 388197, 38318.82 388005.312...
3	-12.913499	POLYGON ((38863.301 387844.906, 38520.699 3876...
4	-12.451985	POLYGON ((39549.359 386816.688, 39243.75 38663...

In [23]:

# Plot the GeoDataFrame
fig, ax = plt.subplots(figsize=(10, 6))
gdf_faces.plot(
    column="bathymetry",
    cmap="terrain",
    ax=ax,
    legend=True,
    edgecolor="gray",
    linewidth=0.1,
)
ax.set_title("Face Polygons as GeoDataFrame")
ax.set_aspect("equal")
plt.tight_layout()
plt.show()

# Plot the GeoDataFrame fig, ax = plt.subplots(figsize=(10, 6)) gdf_faces.plot( column="bathymetry", cmap="terrain", ax=ax, legend=True, edgecolor="gray", linewidth=0.1, ) ax.set_title("Face Polygons as GeoDataFrame") ax.set_aspect("equal") plt.tight_layout() plt.show()

In [24]:

# Nodes as Points
gdf_nodes = subset.to_geodataframe("velocity", location="node")
print(f"Node GeoDataFrame: {len(gdf_nodes)} points")
gdf_nodes.head()

# Nodes as Points gdf_nodes = subset.to_geodataframe("velocity", location="node") print(f"Node GeoDataFrame: {len(gdf_nodes)} points") gdf_nodes.head()

Node GeoDataFrame: 3987 points

Out[24]:

	velocity	geometry
0	0.262553	POINT (38092.773 389066.469)
1	0.298695	POINT (38331.969 388645.5)
2	0.336605	POINT (38622.91 388197)
3	0.365932	POINT (38863.301 387844.906)
4	0.395183	POINT (39074.945 387492.125)

In [25]:

# FeatureCollection (pyramids native vector type)
fc = subset.to_feature_collection("bathymetry", location="face")
fc
print(f"FeatureCollection: {len(fc)} features")
print(type(fc))

# FeatureCollection (pyramids native vector type) fc = subset.to_feature_collection("bathymetry", location="face") fc print(f"FeatureCollection: {len(fc)} features") print(type(fc))

FeatureCollection: 3729 features
<class 'pyramids.feature.feature.FeatureCollection'>

10. CRS Reprojection¶

In [26]:

# Reproject from Dutch national grid (EPSG:28992) to WGS84 (EPSG:4326)
ds_wgs84 = synth_ds.to_crs(4326)
print(f"Original EPSG: {synth_ds.epsg}")
print(f"Reprojected EPSG: {ds_wgs84.epsg}")
print(f"Original bounds: {synth_ds.bounds}")
print(f"WGS84 bounds:    {ds_wgs84.bounds}")
print()
print(
    f"Data preserved: {np.array_equal(synth_ds['bathymetry'].data, ds_wgs84['bathymetry'].data)}"
)

# Reproject from Dutch national grid (EPSG:28992) to WGS84 (EPSG:4326) ds_wgs84 = synth_ds.to_crs(4326) print(f"Original EPSG: {synth_ds.epsg}") print(f"Reprojected EPSG: {ds_wgs84.epsg}") print(f"Original bounds: {synth_ds.bounds}") print(f"WGS84 bounds: {ds_wgs84.bounds}") print() print( f"Data preserved: {np.array_equal(synth_ds['bathymetry'].data, ds_wgs84['bathymetry'].data)}" )

Original EPSG: 28992
Reprojected EPSG: 4326
Original bounds: (17325.978515625, 359012.23738230974, 86177.27534049794, 392033.78125)
WGS84 bounds:    (3.4077401318623934, 51.216226709146426, 4.401786735741385, 51.50166760826841)

Data preserved: True

11. Writing UGRID-Compliant NetCDF¶

In [27]:

# Write the synthetic dataset to a new UGRID NetCDF file
out_path = Path("ugrid_convention_nc_synthetic.nc")
synth_ds.to_file(out_path)
print(f"Written to {out_path} ({out_path.stat().st_size / 1024:.0f} KB)")

# Read it back and verify
ds_roundtrip = UgridDataset.read_file(out_path)
print(f"\nRound-trip verification:")
print(f"  Nodes: {ds_roundtrip.n_node} (expected {synth_ds.n_node})")
print(f"  Faces: {ds_roundtrip.n_face} (expected {synth_ds.n_face})")
print(f"  Variables: {ds_roundtrip.data_variable_names}")

# Write the synthetic dataset to a new UGRID NetCDF file out_path = Path("ugrid_convention_nc_synthetic.nc") synth_ds.to_file(out_path) print(f"Written to {out_path} ({out_path.stat().st_size / 1024:.0f} KB)") # Read it back and verify ds_roundtrip = UgridDataset.read_file(out_path) print(f"\nRound-trip verification:") print(f" Nodes: {ds_roundtrip.n_node} (expected {synth_ds.n_node})") print(f" Faces: {ds_roundtrip.n_face} (expected {synth_ds.n_face})") print(f" Variables: {ds_roundtrip.data_variable_names}")

Written to ugrid_convention_nc_synthetic.nc (420 KB)

Round-trip verification:
  Nodes: 8916 (expected 8916)
  Faces: 8355 (expected 8355)
  Variables: ['bathymetry', 'velocity']

12. Time Dimension Support & Animation¶

Create a temporal dataset simulating tidal water levels.

In [28]:

# Simulate 12 hours of tidal water levels (hourly)
n_hours = 12
time_hours = np.arange(n_hours)

# Tidal signal: water level varies with time and position
water_levels = np.zeros((n_hours, subset.n_face))
sub_cx, sub_cy = subset.mesh.face_centroids
sub_xmin, sub_ymin, sub_xmax, sub_ymax = subset.bounds
sub_x_norm = (sub_cx - sub_xmin) / (sub_xmax - sub_xmin)

for t in range(n_hours):
    phase = 2 * np.pi * t / 12  # 12-hour tidal cycle
    water_levels[t] = 1.5 * np.sin(phase + sub_x_norm * np.pi) + 0.1 * np.random.randn(
        subset.n_face
    )

print(f"Water level array shape: {water_levels.shape}")
print(f"Range: [{water_levels.min():.2f}, {water_levels.max():.2f}] m")

# Simulate 12 hours of tidal water levels (hourly) n_hours = 12 time_hours = np.arange(n_hours) # Tidal signal: water level varies with time and position water_levels = np.zeros((n_hours, subset.n_face)) sub_cx, sub_cy = subset.mesh.face_centroids sub_xmin, sub_ymin, sub_xmax, sub_ymax = subset.bounds sub_x_norm = (sub_cx - sub_xmin) / (sub_xmax - sub_xmin) for t in range(n_hours): phase = 2 * np.pi * t / 12 # 12-hour tidal cycle water_levels[t] = 1.5 * np.sin(phase + sub_x_norm * np.pi) + 0.1 * np.random.randn( subset.n_face ) print(f"Water level array shape: {water_levels.shape}") print(f"Range: [{water_levels.min():.2f}, {water_levels.max():.2f}] m")

Water level array shape: (12, 3729)
Range: [-1.83, 1.82] m

In [29]:

# Create a temporal UgridDataset
temporal_ds = UgridDataset.create_from_arrays(
    node_x=subset.mesh.node_x,
    node_y=subset.mesh.node_y,
    face_node_connectivity=subset.mesh.face_node_connectivity.data,
    data={"water_level": water_levels},
    data_locations={"water_level": "face"},
    epsg=28992,
)

wl_var = temporal_ds["water_level"]
print(f"Has time dimension: {wl_var.has_time}")
print(f"Time steps: {wl_var.n_time_steps}")
print(f"Shape: {wl_var.shape}")

# Create a temporal UgridDataset temporal_ds = UgridDataset.create_from_arrays( node_x=subset.mesh.node_x, node_y=subset.mesh.node_y, face_node_connectivity=subset.mesh.face_node_connectivity.data, data={"water_level": water_levels}, data_locations={"water_level": "face"}, epsg=28992, ) wl_var = temporal_ds["water_level"] print(f"Has time dimension: {wl_var.has_time}") print(f"Time steps: {wl_var.n_time_steps}") print(f"Shape: {wl_var.shape}")

Has time dimension: True
Time steps: 12
Shape: (12, 3729)

In [30]:

# Select a single time step
ds_t3 = temporal_ds.sel_time(3)
print(f"Time step 3: {ds_t3['water_level'].shape}")
print(f"Has time: {ds_t3['water_level'].has_time}")

# Select a single time step ds_t3 = temporal_ds.sel_time(3) print(f"Time step 3: {ds_t3['water_level'].shape}") print(f"Has time: {ds_t3['water_level'].has_time}")

Time step 3: (3729,)
Has time: False

In [31]:

# Plot multiple time steps side by side
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
for i, ax in enumerate(axes.flat):
    t = i * 2  # every 2 hours
    from pyramids.netcdf.ugrid.plot import plot_mesh_data

    glyph = plot_mesh_data(
        temporal_ds.mesh,
        water_levels[t],
        location="face",
        ax=ax,
        cmap="RdBu_r",
        vmin=-2.0,
        vmax=2.0,
        colorbar=False,
        title=f"Hour {t}",
    )

fig.suptitle("Tidal Water Levels (every 2 hours)", fontsize=14)
plt.tight_layout()
plt.show()

# Plot multiple time steps side by side fig, axes = plt.subplots(2, 3, figsize=(16, 10)) for i, ax in enumerate(axes.flat): t = i * 2 # every 2 hours from pyramids.netcdf.ugrid.plot import plot_mesh_data glyph = plot_mesh_data( temporal_ds.mesh, water_levels[t], location="face", ax=ax, cmap="RdBu_r", vmin=-2.0, vmax=2.0, colorbar=False, title=f"Hour {t}", ) fig.suptitle("Tidal Water Levels (every 2 hours)", fontsize=14) plt.tight_layout() plt.show()

In [32]:

# Animate the tidal cycle using cleopatra's MeshGlyph.animate()
from pyramids.netcdf.ugrid.plot import _mesh_to_glyph

glyph = _mesh_to_glyph(temporal_ds.mesh)
anim = glyph.animate(
    water_levels,
    time=[f"Hour {t}" for t in range(n_hours)],
    location="face",
    cmap="RdBu_r",
    vmin=-2.0,
    vmax=2.0,
    interval=500,
    title="Tidal Water Level",
    cbar_label="Water Level (m)",
)

# In a Jupyter notebook, display with:
from IPython.display import HTML

HTML(anim.to_jshtml())

# Animate the tidal cycle using cleopatra's MeshGlyph.animate() from pyramids.netcdf.ugrid.plot import _mesh_to_glyph glyph = _mesh_to_glyph(temporal_ds.mesh) anim = glyph.animate( water_levels, time=[f"Hour {t}" for t in range(n_hours)], location="face", cmap="RdBu_r", vmin=-2.0, vmax=2.0, interval=500, title="Tidal Water Level", cbar_label="Water Level (m)", ) # In a Jupyter notebook, display with: from IPython.display import HTML HTML(anim.to_jshtml())

Out[32]:

Once Loop Reflect

In [33]:

# Save animation to GIF
# glyph.save_animation("tidal_animation.gif", fps=2)

# Save animation to GIF # glyph.save_animation("tidal_animation.gif", fps=2)

Summary¶

This notebook demonstrated the complete UGRID workflow in pyramids:

Operation	Method	Notes
Read	UgridDataset.read_file()	Auto-detects UGRID topology
Inspect	.mesh, .metadata, ds[var_name]	Lazy-computed properties
Plot	.plot(), .plot_outline()	Via cleopatra MeshGlyph
Clip	.clip(polygon)	STRtree-accelerated
Subset	.subset_by_bounds()	Centroid-based, vectorized
Rasterize	.to_dataset(cell_size=...)	Nearest or linear interpolation
Vector	.to_geodataframe()	Faces→Polygons, Nodes→Points
Reproject	.to_crs(epsg)	Via pyproj Transformer
Write	.to_file()	UGRID-compliant NetCDF
Time	.sel_time(), .sel_time_range()	Temporal subsetting
Animate	MeshGlyph.animate()	Via cleopatra
Create	.create_from_arrays()	Programmatic construction