Topology Detection and I/O

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

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] = []
    all_array_names = rg.GetMDArrayNames() or []

    mesh_var_names: list[str] = []
    for name in all_array_names:
        md_arr = rg.OpenMDArray(name)
        if md_arr is None:
            continue
        attrs = _read_attributes(md_arr)
        cf_role = attrs.get("cf_role", "")
        has_topo = "topology_dimension" in attrs and "node_coordinates" in attrs
        if cf_role == "mesh_topology" or has_topo:
            mesh_var_names.append(name)

    referenced_meshes: set[str] = set()
    for name in all_array_names:
        md_arr = rg.OpenMDArray(name)
        if md_arr is None:
            continue
        attrs = _read_attributes(md_arr)
        mesh_ref = attrs.get("mesh")
        if isinstance(mesh_ref, str) and mesh_ref not in all_array_names:
            referenced_meshes.add(mesh_ref)

    for mesh_name in referenced_meshes:
        if mesh_name not in mesh_var_names:
            md_arr = rg.OpenMDArray(mesh_name)
            if md_arr is not None:
                attrs = _read_attributes(md_arr)
                if "node_coordinates" in attrs:
                    mesh_var_names.append(mesh_name)

    for name in mesh_var_names:
        md_arr = rg.OpenMDArray(name)
        if md_arr is None:
            continue
        topo = _parse_single_topology(rg, name, md_arr)
        if topo is not None:
            topologies.append(topo)

    return topologies

pyramids.netcdf.ugrid.io.write_ugrid_topology(rg, mesh, mesh_name='mesh2d', 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.	'mesh2d'
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 = "mesh2d",
    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._face_x is not None and mesh._face_y is not None:
        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._face_x is not None and mesh._face_y is not None:
        _write_coord_array(rg, f"{mesh_name}_face_x", [n_face_dim], mesh._face_x)
        _write_coord_array(rg, f"{mesh_name}_face_y", [n_face_dim], 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 = {"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)