UGRID GIS Operations¶

This notebook demonstrates the GIS analysis methods available on UgridDataset. These methods follow the kernel + facade pattern: standalone numpy/scipy kernels in focused modules, wired through the _UgridFacades mixin.

We build a small synthetic mesh so every cell runs fast and the examples are self-contained (no external data files needed).

In [1]:

%matplotlib inline

import numpy as np

from pyramids.netcdf.ugrid import UgridDataset

%matplotlib inline import numpy as np from pyramids.netcdf.ugrid import UgridDataset

1. Build a Synthetic Mesh¶

We create a simple 4×4 quad mesh (16 faces) with synthetic bathymetry data so every example below has something to work with.

In [2]:

nx, ny = 5, 5
x = np.linspace(0, 400, nx)
y = np.linspace(0, 400, ny)
xx, yy = np.meshgrid(x, y)
node_x = xx.ravel().astype(np.float64)
node_y = yy.ravel().astype(np.float64)

faces = []
for j in range(ny - 1):
    for i in range(nx - 1):
        bl = j * nx + i
        br = bl + 1
        tl = bl + nx
        tr = tl + 1
        faces.append([bl, br, tr, tl])
fnc = np.array(faces, dtype=np.intp)

depth = np.array(
    [1.0, 2.0, 3.0, 4.0, 2.0, np.nan, 5.0, 6.0, 3.0, 4.0, 7.0, 8.0, 4.0, 5.0, 6.0, 9.0],
    dtype=np.float64,
)

ds = UgridDataset.from_arrays(
    node_x=node_x,
    node_y=node_y,
    face_node_connectivity=fnc,
    data={"depth": depth},
    data_locations={"depth": "face"},
    epsg=28992,
    mesh_name="mesh2d",
)
print(ds)
print(f"Faces: {ds.n_face}, Nodes: {ds.n_node}, CRS: {ds.epsg}")

nx, ny = 5, 5 x = np.linspace(0, 400, nx) y = np.linspace(0, 400, ny) xx, yy = np.meshgrid(x, y) node_x = xx.ravel().astype(np.float64) node_y = yy.ravel().astype(np.float64) faces = [] for j in range(ny - 1): for i in range(nx - 1): bl = j * nx + i br = bl + 1 tl = bl + nx tr = tl + 1 faces.append([bl, br, tr, tl]) fnc = np.array(faces, dtype=np.intp) depth = np.array( [1.0, 2.0, 3.0, 4.0, 2.0, np.nan, 5.0, 6.0, 3.0, 4.0, 7.0, 8.0, 4.0, 5.0, 6.0, 9.0], dtype=np.float64, ) ds = UgridDataset.from_arrays( node_x=node_x, node_y=node_y, face_node_connectivity=fnc, data={"depth": depth}, data_locations={"depth": "face"}, epsg=28992, mesh_name="mesh2d", ) print(ds) print(f"Faces: {ds.n_face}, Nodes: {ds.n_node}, CRS: {ds.epsg}")

UgridDataset: (in-memory)
  Mesh: mesh2d
  Nodes: 25, Faces: 16, Edges: 0
  Bounds: (0.0, 0.0, 400.0, 400.0)
  CRS: 28992
  Data variables (1):
    depth: location=face, shape=(16,)
Faces: 16, Nodes: 25, CRS: 28992

2. CRS Management¶

set_crs — assign / override a CRS without reprojecting¶

from_arrays defaults to EPSG:4326 (WGS 84). If the coordinates are actually in a local grid, use set_crs to correct the CRS declaration without touching the coordinates.

In [3]:

ds_default_crs = UgridDataset.from_arrays(
    node_x=node_x,
    node_y=node_y,
    face_node_connectivity=fnc,
    data={"depth": depth},
    data_locations={"depth": "face"},
    mesh_name="mesh2d",
)
print(f"Before (default): CRS = EPSG:{ds_default_crs.epsg}")

ds_default_crs.set_crs(epsg=28992)
print(f"After set_crs:    CRS = EPSG:{ds_default_crs.epsg}")

ds_default_crs = UgridDataset.from_arrays( node_x=node_x, node_y=node_y, face_node_connectivity=fnc, data={"depth": depth}, data_locations={"depth": "face"}, mesh_name="mesh2d", ) print(f"Before (default): CRS = EPSG:{ds_default_crs.epsg}") ds_default_crs.set_crs(epsg=28992) print(f"After set_crs: CRS = EPSG:{ds_default_crs.epsg}")

Before (default): CRS = EPSG:4326
After set_crs:    CRS = EPSG:28992

to_crs / to_epsg — reproject mesh coordinates¶

In [4]:

ds_wgs84 = ds.to_crs(4326)
print(f"Original bounds (EPSG:28992): {ds.bounds}")
print(f"Reprojected bounds (WGS 84): {ds_wgs84.bounds}")
print(f"New CRS: {ds_wgs84.epsg}")

ds_wgs84 = ds.to_crs(4326) print(f"Original bounds (EPSG:28992): {ds.bounds}") print(f"Reprojected bounds (WGS 84): {ds_wgs84.bounds}") print(f"New CRS: {ds_wgs84.epsg}")

Original bounds (EPSG:28992): (0.0, 0.0, 400.0, 400.0)
Reprojected bounds (WGS 84): (3.3134090710064177, 47.97476588466149, 3.3189057210154167, 47.9784566419206)
New CRS: 4326

3. No-Data Handling¶

fillna — replace missing values¶

In [5]:

filled = ds.fillna(-9999.0)
original = ds.get_data("depth").data
filled_arr = filled.get_data("depth").data
print(f"Original (has NaN): {original}")
print(f"Filled:             {filled_arr}")

filled = ds.fillna(-9999.0) original = ds.get_data("depth").data filled_arr = filled.get_data("depth").data print(f"Original (has NaN): {original}") print(f"Filled: {filled_arr}")

Original (has NaN): [ 1.  2.  3.  4.  2. nan  5.  6.  3.  4.  7.  8.  4.  5.  6.  9.]
Filled:             [ 1.000e+00  2.000e+00  3.000e+00  4.000e+00  2.000e+00 -9.999e+03
  5.000e+00  6.000e+00  3.000e+00  4.000e+00  7.000e+00  8.000e+00
  4.000e+00  5.000e+00  6.000e+00  9.000e+00]

isnull / notnull — boolean masks¶

In [6]:

null_mask = ds.isnull()
valid_mask = ds.notnull()
print(f"isnull:  {null_mask.get_data('depth').data}")
print(f"notnull: {valid_mask.get_data('depth').data}")
print(f"NaN count: {null_mask.get_data('depth').data.sum()}")

null_mask = ds.isnull() valid_mask = ds.notnull() print(f"isnull: {null_mask.get_data('depth').data}") print(f"notnull: {valid_mask.get_data('depth').data}") print(f"NaN count: {null_mask.get_data('depth').data.sum()}")

isnull:  [False False False False False  True False False False False False False
 False False False False]
notnull: [ True  True  True  True  True False  True  True  True  True  True  True
  True  True  True  True]
NaN count: 1

4. Value Transforms¶

clip_values — clamp to a range¶

In [7]:

clamped = ds.clip_values(2.0, 7.0)
print(f"Original: {ds.get_data('depth').data}")
print(f"Clamped:  {clamped.get_data('depth').data}")

clamped = ds.clip_values(2.0, 7.0) print(f"Original: {ds.get_data('depth').data}") print(f"Clamped: {clamped.get_data('depth').data}")

Original: [ 1.  2.  3.  4.  2. nan  5.  6.  3.  4.  7.  8.  4.  5.  6.  9.]
Clamped:  [ 2.  2.  3.  4.  2. nan  5.  6.  3.  4.  7.  7.  4.  5.  6.  7.]

round — round to N decimals¶

In [8]:

noisy = ds.fillna(0.0)
noisy = noisy.with_variable(
    "depth",
    noisy.get_data("depth").data + np.random.default_rng(42).normal(0, 0.1, 16),
    location="face",
)
rounded = noisy.round(1)
print(f"Noisy:   {noisy.get_data('depth').data[:6]}")
print(f"Rounded: {rounded.get_data('depth').data[:6]}")

noisy = ds.fillna(0.0) noisy = noisy.with_variable( "depth", noisy.get_data("depth").data + np.random.default_rng(42).normal(0, 0.1, 16), location="face", ) rounded = noisy.round(1) print(f"Noisy: {noisy.get_data('depth').data[:6]}") print(f"Rounded: {rounded.get_data('depth').data[:6]}")

Noisy:   [ 1.03047171  1.89600159  3.07504512  4.09405647  1.80489648 -0.13021795]
Rounded: [ 1.   1.9  3.1  4.1  1.8 -0.1]

astype — cast to a different dtype¶

In [9]:

filled_int = filled.astype(np.int32)
print(f"Original dtype: {filled.get_data('depth').data.dtype}")
print(f"Cast dtype:     {filled_int.get_data('depth').data.dtype}")
print(f"Values:         {filled_int.get_data('depth').data}")

filled_int = filled.astype(np.int32) print(f"Original dtype: {filled.get_data('depth').data.dtype}") print(f"Cast dtype: {filled_int.get_data('depth').data.dtype}") print(f"Values: {filled_int.get_data('depth').data}")

Original dtype: float64
Cast dtype:     int32
Values:         [    1     2     3     4     2 -9999     5     6     3     4     7     8
     4     5     6     9]

isin — membership test¶

In [10]:

membership = ds.isin([1.0, 3.0, 5.0, 7.0, 9.0])
print(f"Values:    {ds.get_data('depth').data}")
print(f"Is in set: {membership.get_data('depth').data}")

membership = ds.isin([1.0, 3.0, 5.0, 7.0, 9.0]) print(f"Values: {ds.get_data('depth').data}") print(f"Is in set: {membership.get_data('depth').data}")

Values:    [ 1.  2.  3.  4.  2. nan  5.  6.  3.  4.  7.  8.  4.  5.  6.  9.]
Is in set: [ True False  True False False False  True False  True False  True False
 False  True False  True]

5. Domain Metrics¶

domain_area — total mesh area¶

In [11]:

area = ds.domain_area()
print(f"Total mesh area: {area:.0f} square units")
print(f"Expected (400×400): {400 * 400}")

area = ds.domain_area() print(f"Total mesh area: {area:.0f} square units") print(f"Expected (400×400): {400 * 400}")

Total mesh area: 160000 square units
Expected (400×400): 160000

count_domain_cells — number of active faces¶

In [12]:

n_all = ds.count_domain_cells()
print(f"Total faces: {n_all}")
print(f"Expected: {(nx - 1) * (ny - 1)}")

n_all = ds.count_domain_cells() print(f"Total faces: {n_all}") print(f"Expected: {(nx - 1) * (ny - 1)}")

Total faces: 16
Expected: 16

6. Dataset Comparison¶

equals — value equality (ignores dtype)¶

In [13]:

ds_copy = ds.copy()
print(f"ds.equals(ds_copy): {ds.equals(ds_copy)}")

ds_modified = ds.fillna(0.0)
print(f"ds.equals(ds_modified): {ds.equals(ds_modified)}")

ds_copy = ds.copy() print(f"ds.equals(ds_copy): {ds.equals(ds_copy)}") ds_modified = ds.fillna(0.0) print(f"ds.equals(ds_modified): {ds.equals(ds_modified)}")

ds.equals(ds_copy): True
ds.equals(ds_modified): False

identical — strict equality (dtype, CRS, nodata, topology)¶

In [14]:

print(f"ds.identical(ds_copy): {ds.identical(ds_copy)}")

ds_int = ds.fillna(0.0).astype(np.int32)
ds_float = ds.fillna(0.0)
print(f"Same values, different dtype: {ds_int.identical(ds_float)}")

print(f"ds.identical(ds_copy): {ds.identical(ds_copy)}") ds_int = ds.fillna(0.0).astype(np.int32) ds_float = ds.fillna(0.0) print(f"Same values, different dtype: {ds_int.identical(ds_float)}")

ds.identical(ds_copy): True
Same values, different dtype: False

7. Point Extraction¶

extract — sample values at specific coordinates¶

In [15]:

query_x = np.array([50.0, 150.0, 250.0, 350.0])
query_y = np.array([50.0, 150.0, 250.0, 350.0])
values = ds.extract(query_x, query_y)
print("Extracted values at diagonal points:")
for xi, yi, vi in zip(query_x, query_y, values):
    print(f"  ({xi:.0f}, {yi:.0f}) → {vi}")

query_x = np.array([50.0, 150.0, 250.0, 350.0]) query_y = np.array([50.0, 150.0, 250.0, 350.0]) values = ds.extract(query_x, query_y) print("Extracted values at diagonal points:") for xi, yi, vi in zip(query_x, query_y, values): print(f" ({xi:.0f}, {yi:.0f}) → {vi}")

Extracted values at diagonal points:
  (50, 50) → 1.0
  (150, 150) → nan
  (250, 250) → 7.0
  (350, 350) → 9.0

sample_random — random spatial samples¶

In [16]:

rx, ry, rvals = ds.sample_random(5, seed=42)
print("Random samples (n=5, seed=42):")
for xi, yi, vi in zip(rx, ry, rvals):
    print(f"  ({xi:.1f}, {yi:.1f}) → {vi}")

rx, ry, rvals = ds.sample_random(5, seed=42) print("Random samples (n=5, seed=42):") for xi, yi, vi in zip(rx, ry, rvals): print(f" ({xi:.1f}, {yi:.1f}) → {vi}")

Random samples (n=5, seed=42):
  (309.6, 326.6) → 9.0
  (175.6, 42.1) → 2.0
  (343.4, 26.6) → 4.0
  (278.9, 237.8) → 7.0
  (37.7, 58.5) → 1.0

8. Graph Operations¶

proximity — graph distance from target faces¶

In [17]:

target = np.array([0])
dist = ds.proximity(target)
print("Distance from face 0 (top-left):")
print(dist.reshape(4, 4))

target = np.array([0]) dist = ds.proximity(target) print("Distance from face 0 (top-left):") print(dist.reshape(4, 4))

Distance from face 0 (top-left):
[[0 1 2 3]
 [1 2 3 4]
 [2 3 4 5]
 [3 4 5 6]]

sieve — remove small connected regions¶

In [18]:

island_data = np.full(16, np.nan)
island_data[0] = 1.0
island_data[4:8] = [2.0, 3.0, 4.0, 5.0]
island_data[8:12] = [6.0, 7.0, 8.0, 9.0]

ds_island = ds.with_variable("depth", island_data, location="face")
print(f"Before sieve: {ds_island.get_data('depth').data}")

sieved = ds_island.sieve(threshold=3)
print(f"After sieve (threshold=3): {sieved.get_data('depth').data}")
print("(The isolated face 0 was blanked to NaN)")

island_data = np.full(16, np.nan) island_data[0] = 1.0 island_data[4:8] = [2.0, 3.0, 4.0, 5.0] island_data[8:12] = [6.0, 7.0, 8.0, 9.0] ds_island = ds.with_variable("depth", island_data, location="face") print(f"Before sieve: {ds_island.get_data('depth').data}") sieved = ds_island.sieve(threshold=3) print(f"After sieve (threshold=3): {sieved.get_data('depth').data}") print("(The isolated face 0 was blanked to NaN)")

Before sieve: [ 1. nan nan nan  2.  3.  4.  5.  6.  7.  8.  9. nan nan nan nan]
After sieve (threshold=3): [ 1. nan nan nan  2.  3.  4.  5.  6.  7.  8.  9. nan nan nan nan]
(The isolated face 0 was blanked to NaN)

9. Contouring¶

contour — extract iso-lines as GeoDataFrame¶

In [19]:

ds_clean = ds.fillna(5.0)
gdf = ds_clean.contour(levels=[2.0, 4.0, 6.0, 8.0])
print("Contour GeoDataFrame:")
print(gdf)
print(f"\nGeometry types: {gdf.geometry.geom_type.unique()}")

ds_clean = ds.fillna(5.0) gdf = ds_clean.contour(levels=[2.0, 4.0, 6.0, 8.0]) print("Contour GeoDataFrame:") print(gdf) print(f"\nGeometry types: {gdf.geometry.geom_type.unique()}")

Contour GeoDataFrame:
   level                                           geometry
0    2.0  LINESTRING (0 150, 25 125, 50 100, 66.667 66.6...
1    4.0  LINESTRING (0 400, 50 350, 100 300, 125 225, 1...
2    6.0  LINESTRING (225 400, 225 325, 225 300, 233.333...
3    8.0  LINESTRING (333.333 400, 333.333 333.333, 350 ...

Geometry types: <ArrowStringArray>
['LineString']
Length: 1, dtype: str

10. Plotting¶

plot_histogram — value distribution¶

In [20]:

fig, ax, hist = ds.plot_histogram(bins=8)
ax.set_title("Depth distribution")
ax.set_xlabel("Depth")
ax.set_ylabel("Count")

fig, ax, hist = ds.plot_histogram(bins=8) ax.set_title("Depth distribution") ax.set_xlabel("Depth") ax.set_ylabel("Count")

Out[20]:

Text(0, 0.5, 'Count')

11. Serialization¶

to_bytes — in-memory NetCDF round-trip¶

In [21]:

blob = ds.to_bytes()
print(f"Serialized size: {len(blob)} bytes")

ds_roundtrip = UgridDataset.from_bytes(blob)
print(f"Round-trip equals: {ds.equals(ds_roundtrip)}")
print(f"Round-trip CRS:    {ds_roundtrip.epsg}")

blob = ds.to_bytes() print(f"Serialized size: {len(blob)} bytes") ds_roundtrip = UgridDataset.from_bytes(blob) print(f"Round-trip equals: {ds.equals(ds_roundtrip)}") print(f"Round-trip CRS: {ds_roundtrip.epsg}")

Serialized size: 15856 bytes
Round-trip equals: True
Round-trip CRS:    28992

Summary¶

This notebook demonstrated 17 GIS-parity methods added to UgridDataset:

Category	Methods
CRS	set_crs, to_crs, to_epsg
No-data	fillna, isnull, notnull
Value transforms	clip_values, round, astype, isin
Domain metrics	domain_area, count_domain_cells
Comparison	equals, identical
Point extraction	extract, sample_random
Graph ops	proximity, sieve
Contouring	contour
Plotting	plot_histogram
I/O	to_bytes

For the complete list of all ~95 methods available on UgridDataset (including inherited mixin methods), see the Method Reference.