riverarchitect.volume

Triangulated-surface volume computation for earthworks quantities.

Replaces arcpy.SurfaceVolume_3d (3D Analyst) with an equivalent, licence-free implementation.

Why triangulated and not a cell sum

Summing z * cell_area treats every cell as a vertical prism, which over- or underestimates by a boundary term proportional to the perimeter of the surface. On a flat 10 x 10 m test surface the prism sum overestimates by 21 %, and for surfaces truncated at their edge that error does not vanish under grid refinement. SurfaceVolume_3d instead integrates under the triangulated surface through the cell centres. River Architect reports earthworks quantities as a deliverable, so the triangulated result is authoritative.

The surface is built exactly as a TIN over the raster cell centres: every 2x2 block of neighbouring centres forms a quad that is split into two triangles. Over each triangle the elevation varies linearly, and the volume between the triangle and a horizontal reference plane is integrated in closed form, clipping the triangle where it crosses the plane.

This module is pure numpy: no arcpy, no GDAL, testable anywhere.

riverarchitect.volume.surface_volume(z, dx, dy, plane=0.0, reference='ABOVE')[source]

Volume between a raster surface and a horizontal reference plane.

Equivalent to arcpy.SurfaceVolume_3d(in_surface, “”, reference, plane, 1.0).

Parameters:

  • z (np.ndarray) – 2-D elevation array. NoData cells must be np.nan.

  • dx (float) – cell size in x (map units).

  • dy (float) – cell size in y (map units).

  • plane (float) – reference plane elevation.

  • reference (str) – “ABOVE” for the volume between the plane and the surface where the surface lies above it, “BELOW” for the volume where it lies below.

Returns:

{“volume”: float, “area_2d”: float, “area_3d”: float}

volume - cubic map units area_2d - planimetric area contributing to the volume area_3d - surface (draped) area contributing to the volume

Return type:

dict

Triangles touching a NoData cell are omitted, matching arcpy’s behaviour of building the TIN only over cells that carry data.