Topology Detection and I/O

UGRID topology detection from NetCDF files using GDAL's Multidimensional API, and UGRID-compliant NetCDF writing.

On read, parse_ugrid_topology scans the root group for mesh topology variables and returns one MeshTopologyInfo per mesh, which Mesh2d.from_gdal_group turns into topology. On write, the topology and each data variable are emitted back as CF/UGRID-compliant arrays:

Hold "Ctrl" to enable pan & zoom

flowchart LR
    F[(".nc root group")] -->|parse_ugrid_topology| TI["list[MeshTopologyInfo]"]
    TI -->|"Mesh2d.from_gdal_group"| M["Mesh2d"]
    M -->|write_ugrid_topology| F
    V["MeshVariable"] -->|write_ugrid_data_variable| F

pyramids.netcdf.ugrid.io.parse_ugrid_topology(rg)

Detect and parse all UGRID mesh topologies from a GDAL root group.

Detection strategy (handles diverse real-world files):

1. Primary: Scan all MDArrays for cf_role = "mesh_topology" attribute.
2. Fallback: Scan for variables with topology_dimension AND node_coordinates attributes (older files without cf_role).
3. Scalar check: GDAL may filter 0-dimensional arrays from GetMDArrayNames(). Explicitly try OpenMDArray(name) for variable names found in other variables' mesh attributes.

Parameters:

Name	Type	Description	Default
rg	Group	GDAL root group from a multidimensional NetCDF file.	required

Returns:

Type	Description
list[MeshTopologyInfo]	List of MeshTopologyInfo, one per mesh found in the file.
list[MeshTopologyInfo]	Most files have exactly one mesh. Empty list if no UGRID topology.

Source code in src/pyramids/netcdf/ugrid/io.py

def parse_ugrid_topology(rg: gdal.Group) -> list[MeshTopologyInfo]:
    """Detect and parse all UGRID mesh topologies from a GDAL root group.

    Detection strategy (handles diverse real-world files):

    1. Primary: Scan all MDArrays for `cf_role = "mesh_topology"` attribute.
    2. Fallback: Scan for variables with `topology_dimension` AND
       `node_coordinates` attributes (older files without cf_role).
    3. Scalar check: GDAL may filter 0-dimensional arrays from
       `GetMDArrayNames()`. Explicitly try `OpenMDArray(name)` for
       variable names found in other variables' `mesh` attributes.

    Args:
        rg: GDAL root group from a multidimensional NetCDF file.

    Returns:
        List of MeshTopologyInfo, one per mesh found in the file.
        Most files have exactly one mesh. Empty list if no UGRID topology.
    """
    topologies: list[MeshTopologyInfo] = []
    scan = _MeshArrayScan(rg)

    for name, md_arr in _find_mesh_topology_arrays(scan).items():
        topo = _parse_single_topology(scan, name, md_arr)
        if topo is not None:
            topologies.append(topo)

    # Fallback: UXARRAY-style files declare topology purely via cf_role on the connectivity variable
    # (no central mesh_topology variable). Infer a face/node mesh from it (#589).
    if not topologies:
        inferred = _infer_topology_from_connectivity(scan)
        if inferred is not None:
            topologies.append(inferred)

    return topologies

pyramids.netcdf.ugrid.io.write_ugrid_topology(rg, mesh, mesh_name=DEFAULT_MESH_NAME, crs_wkt=None)

Write UGRID topology to a GDAL root group.

Creates the topology variable, node coordinate arrays, connectivity arrays, and optional face/edge center coordinates. Uses cf.write_attributes_to_md_array for all attribute writing.

Parameters:

Name	Type	Description	Default
rg	Group	GDAL root group to write into.	required
mesh	Mesh2d	Mesh2d instance.	required
mesh_name	str	Name for the topology variable.	DEFAULT_MESH_NAME
crs_wkt	str | None	WKT CRS string (optional).	None

Returns:

Type	Description
dict[str, Any]	Dict mapping dimension names to GDAL dimension objects,
dict[str, Any]	for use when writing data variables.

Source code in src/pyramids/netcdf/ugrid/io.py

def write_ugrid_topology(
    rg: gdal.Group,
    mesh: Mesh2d,
    mesh_name: str = DEFAULT_MESH_NAME,
    crs_wkt: str | None = None,
) -> dict[str, Any]:
    """Write UGRID topology to a GDAL root group.

    Creates the topology variable, node coordinate arrays,
    connectivity arrays, and optional face/edge center coordinates.
    Uses cf.write_attributes_to_md_array for all attribute writing.

    Args:
        rg: GDAL root group to write into.
        mesh: Mesh2d instance.
        mesh_name: Name for the topology variable.
        crs_wkt: WKT CRS string (optional).

    Returns:
        Dict mapping dimension names to GDAL dimension objects,
        for use when writing data variables.
    """
    dims: dict[str, Any] = {}

    n_node_dim = rg.CreateDimension(f"{mesh_name}_nNodes", None, None, mesh.n_node)
    n_face_dim = rg.CreateDimension(f"{mesh_name}_nFaces", None, None, mesh.n_face)
    dims[f"{mesh_name}_nNodes"] = n_node_dim
    dims[f"{mesh_name}_nFaces"] = n_face_dim

    fnc = mesh.face_node_connectivity
    max_fn_dim = rg.CreateDimension(
        f"{mesh_name}_nMaxFaceNodes",
        None,
        None,
        fnc.max_nodes_per_element,
    )
    two_dim = rg.CreateDimension("Two", None, None, 2)

    topo_dim = rg.CreateDimension(f"{mesh_name}_scalar", None, None, 1)
    topo_arr = rg.CreateMDArray(
        mesh_name,
        [topo_dim],
        gdal.ExtendedDataType.Create(gdal.GDT_Int32),
    )
    topo_attrs = {
        "cf_role": "mesh_topology",
        "topology_dimension": 2,
        "node_coordinates": f"{mesh_name}_node_x {mesh_name}_node_y",
        "face_node_connectivity": f"{mesh_name}_face_nodes",
    }
    if mesh.edge_node_connectivity is not None:
        topo_attrs["edge_node_connectivity"] = f"{mesh_name}_edge_nodes"
    if mesh.has_face_coords:
        topo_attrs["face_coordinates"] = f"{mesh_name}_face_x {mesh_name}_face_y"
    write_attributes_to_md_array(topo_arr, topo_attrs)
    topo_arr.Write(np.array([0], dtype=np.int32))

    _write_coord_array(rg, f"{mesh_name}_node_x", [n_node_dim], mesh.node_x)
    _write_coord_array(rg, f"{mesh_name}_node_y", [n_node_dim], mesh.node_y)

    _write_connectivity_array(
        rg,
        f"{mesh_name}_face_nodes",
        [n_face_dim, max_fn_dim],
        fnc,
    )

    if mesh.edge_node_connectivity is not None:
        enc = mesh.edge_node_connectivity
        n_edge_dim = rg.CreateDimension(
            f"{mesh_name}_nEdges",
            None,
            None,
            enc.n_elements,
        )
        dims[f"{mesh_name}_nEdges"] = n_edge_dim
        _write_connectivity_array(
            rg,
            f"{mesh_name}_edge_nodes",
            [n_edge_dim, two_dim],
            enc,
        )

    if mesh.has_face_coords:
        _write_coord_array(
            rg, f"{mesh_name}_face_x", [n_face_dim], cast("np.typing.NDArray", mesh.face_x)
        )
        _write_coord_array(
            rg, f"{mesh_name}_face_y", [n_face_dim], cast("np.typing.NDArray", mesh.face_y)
        )

    if crs_wkt is not None:
        _write_crs_variable(rg, crs_wkt, topo_dim)

    return dims

pyramids.netcdf.ugrid.io.write_ugrid_data_variable(rg, var, mesh_name, dims)

Write a single data variable to the GDAL group.

Parameters:

Name	Type	Description	Default
rg	Group	GDAL root group.	required
var	MeshVariable	MeshVariable instance.	required
mesh_name	str	Name of the mesh topology variable.	required
dims	dict[str, Any]	Dict mapping dimension names to GDAL dimension objects.	required

Source code in src/pyramids/netcdf/ugrid/io.py

def write_ugrid_data_variable(
    rg: gdal.Group,
    var: MeshVariable,
    mesh_name: str,
    dims: dict[str, Any],
) -> None:
    """Write a single data variable to the GDAL group.

    Args:
        rg: GDAL root group.
        var: MeshVariable instance.
        mesh_name: Name of the mesh topology variable.
        dims: Dict mapping dimension names to GDAL dimension objects.
    """
    if var.data is None:
        return

    dim_list = []
    if var.has_time and "time" in dims:
        dim_list.append(dims["time"])

    loc_dim_name = f"{mesh_name}_n{var.location.capitalize()}s"
    if loc_dim_name in dims:
        dim_list.append(dims[loc_dim_name])
    else:
        loc_dim = rg.CreateDimension(loc_dim_name, None, None, var.n_elements)
        dims[loc_dim_name] = loc_dim
        dim_list.append(loc_dim)

    dtype_map = {
        np.dtype("float64"): gdal.GDT_Float64,
        np.dtype("float32"): gdal.GDT_Float32,
        np.dtype("int64"): gdal.GDT_Int64,
        np.dtype("int32"): gdal.GDT_Int32,
        np.dtype("int16"): gdal.GDT_Int16,
        np.dtype("int8"): gdal.GDT_Int16,
        np.dtype("uint8"): gdal.GDT_Byte,
        np.dtype("uint16"): gdal.GDT_UInt16,
        np.dtype("uint32"): gdal.GDT_UInt32,
    }
    gdal_dt = dtype_map.get(var.dtype, gdal.GDT_Float64)

    md_arr = rg.CreateMDArray(
        var.name,
        dim_list,
        gdal.ExtendedDataType.Create(gdal_dt),
    )
    md_arr.Write(var.data)

    var_attrs: dict[str, Any] = {"mesh": mesh_name, "location": var.location}
    if var.units:
        var_attrs["units"] = var.units
    if var.standard_name:
        var_attrs["standard_name"] = var.standard_name
    if var.nodata is not None:
        var_attrs["_FillValue"] = var.nodata

    write_attributes_to_md_array(md_arr, var_attrs)