riverarchitect.raster

Raster operations for River Architect, built on rasterio, GDAL, numpy and scipy.

This module is the open-source replacement for the arcpy / Spatial Analyst raster algebra that the original River Architect was written against. Every function here has a named arcpy counterpart, documented in Migrating from arcpy.

Two conventions run through the whole module and matter more than any individual function:

NoData is always ``numpy.nan`` in memory.

Rasters are read with the declared NoData mask applied and converted to float64. Never compare cell values against a sentinel such as -999; third-party inputs carry -3.4e38, 3.4e38, 0 or -9999 interchangeably. On write, the canonical value riverarchitect.config.NODATA is stamped back in.

Alignment is explicit.

arcpy silently reconciled operands of differing extent through arcpy.env.extent and arcpy.env.snapRaster. numpy does not: it either raises on a shape mismatch or, worse, broadcasts into a spatially meaningless result. Use align() before combining any two rasters that do not provably share a grid.

Example

Compute a water surface elevation from a DEM and a depth raster whose extents differ:

from riverarchitect import raster

dem, dem_prof = raster.read("dem.tif")
h, h_prof = raster.read("h000098_750.tif")
h = raster.align(h, h_prof, dem_prof)          # onto the DEM grid
wse = raster.con(h > 0, dem + h)               # false branch becomes NoData
raster.write("wse.tif", wse, dem_prof)

riverarchitect.raster.read(path, band=1)

Read a raster into a float array with NoData as numpy.nan.

Replaces arcpy.Raster(path) and arcpy.RasterToNumPyArray.

Parameters:

• path (str) – path to a GDAL-readable raster.

• band (int) – 1-based band index.

Returns:

(numpy.ndarray, dict) - the array and the rasterio profile.

Return type:

tuple

riverarchitect.raster.write(path, array, prof, dtype='float32', nodata=None, compress='lzw')

Write array to path, stamping in the canonical NoData value.

Replaces arcpy.Raster.save() and arcpy.NumPyArrayToRaster.

Parameters:

• path (str) – output path (GeoTIFF).

• array (numpy.ndarray) – 2-D array; numpy.nan marks NoData.

• prof (dict) – rasterio profile supplying the grid and CRS.

• dtype (str) – output data type.

• nodata (float) – NoData value; defaults to config.NODATA.

• compress (str) – GeoTIFF compression.

Returns:

path.

Return type:

str

riverarchitect.raster.profile_of(path)

Return the rasterio profile of path without reading pixel data.

riverarchitect.raster.cell_size(prof)

Return (dx, dy) in map units.

Replaces GetRasterProperties_management(ras, "CELLSIZEX"), which returned a localised string and so broke arithmetic on systems using a comma decimal separator.

riverarchitect.raster.align(array, src_prof, ref_prof, resampling=rasterio.enums.Resampling.nearest)

Resample array onto the grid described by ref_prof.

This is the explicit form of arcpy.env.extent + env.snapRaster + env.cellSize. Use Resampling.nearest (the default, and arcpy’s default) for anything categorical or where original cell values must survive; Resampling.bilinear only for genuinely continuous resampling.

Parameters:

• array (numpy.ndarray) – source array with numpy.nan NoData.

• src_prof (dict) – profile describing array.

• ref_prof (dict) – profile describing the target grid.

• resampling (rasterio.enums.Resampling) – resampling method.

Returns:

array on the reference grid, NoData as numpy.nan.

Return type:

numpy.ndarray

riverarchitect.raster.resample(path, factor, resampling=rasterio.enums.Resampling.bilinear)

Read path at a coarser or finer resolution.

Replaces arcpy.Resample_management. factor > 1 coarsens.

riverarchitect.raster.con(condition, true_value, false_value=numpy.nan)

Conditional evaluation. Replaces arcpy.sa.Con.

Note the default: the two-argument form of arcpy’s Con yields NoData, not zero, where the condition is false. Passing false_value=0 instead is a common and hard-to-spot source of wetted-area artefacts.

riverarchitect.raster.is_null(array)

Replaces arcpy.sa.IsNull.

riverarchitect.raster.set_null(condition, array)

Replaces arcpy.sa.SetNull: sets cells to NoData where condition holds.

riverarchitect.raster.cell_statistics(arrays, statistic='MAXIMUM', ignore_nodata=True)

Per-cell statistic across a stack of aligned rasters.

Replaces arcpy.sa.CellStatistics(rasters, statistic, "DATA"). ignore_nodata=True corresponds to "DATA", which is what every call site in the original code used.

Parameters:

• arrays (list) – aligned 2-D arrays of identical shape.

• statistic (str) – MAXIMUM, MINIMUM, MEAN, SUM, MEDIAN, STD or RANGE.

• ignore_nodata (bool) – ignore NoData cells instead of propagating them.

Returns:

numpy.ndarray

riverarchitect.raster.reclassify(array, breaks, values, right=True)

Map continuous values onto classes. Replaces arcpy.sa.Reclassify/RemapRange.

Parameters:

• array (numpy.ndarray) – input.

• breaks (list) – ascending class boundaries.

• values (list) – output value per class; len(values) == len(breaks) + 1.

• right (bool) – whether intervals are closed on the right.

Returns:

numpy.ndarray

riverarchitect.raster.polygonize(array, prof, mask=None, connectivity=4)

Vectorise a raster into polygons. Replaces arcpy.RasterToPolygon_conversion.

Parameters:

• array (numpy.ndarray) – values to vectorise (cast to int32, as arcpy requires).

• prof (dict) – rasterio profile supplying transform and CRS.

• mask (numpy.ndarray) – optional boolean mask of cells to include.

• connectivity (int) – 4 (arcpy’s default) or 8.

Returns:

with a gridcode column, mirroring arcpy’s output field.

Return type:

geopandas.GeoDataFrame

riverarchitect.raster.rasterize(gdf, prof, column='gridcode', fill=0, all_touched=False, dtype='int32')

Burn polygons into a raster grid. Replaces arcpy.PolygonToRaster_conversion.

riverarchitect.raster.extract_by_mask(path, geometries, crop=True)

Clip a raster to polygon geometries. Replaces arcpy.sa.ExtractByMask.

riverarchitect.raster.raster_to_points(array, prof, step=1)

Cell centres and values of every valid cell.

Replaces arcpy.RasterToPoint_conversion plus the arcpy.sa.Sample call the original code needed to recover full float precision - the values here are exact by construction.

Parameters:

• array (numpy.ndarray) – source with numpy.nan NoData.

• prof (dict) – rasterio profile.

• step (int) – take every step-th valid cell (subsampling for speed).

Returns:

(points, values) where points is (n, 2) of x/y coordinates.

Return type:

tuple

riverarchitect.raster.idw(points, values, prof, k=12, power=2.0)

Inverse-distance weighted interpolation onto the grid of prof.

Replaces arcpy.Idw_3d(..., search_radius="Variable 12") with its default power of 2.

Parameters:

• points (numpy.ndarray) – (n, 2) source coordinates.

• values (numpy.ndarray) – (n,) source values.

• prof (dict) – profile describing the output grid.

• k (int) – number of nearest neighbours.

• power (float) – inverse-distance exponent.

Returns:

interpolated grid.

Return type:

numpy.ndarray

riverarchitect.raster.nearest_neighbour(points, values, prof)

Nearest-neighbour interpolation.

Matches the original’s “Nearest Neighbor” option, which was Idw_3d with a single neighbour. Values are taken straight from the source array rather than through the inverse-distance weighting, so the output reproduces input values exactly instead of carrying the rounding error of a weight division.

riverarchitect.raster.kriging(points, values, prof, model='spherical', nlags=12, n_closest=12, return_variance=False)

Ordinary kriging onto the grid of prof.

Replaces arcpy.sa.Kriging(..., KrigingModelOrdinary("Spherical", ...)).

Unlike arcpy, the estimation variance is returned on request. The original code had the variance raster commented out because of an arcpy limitation; here it is free.

Kriging cost grows steeply with the number of source points. n_closest keeps the solve local, which is the equivalent of arcpy’s search_radius="Variable 12".

Parameters:

• points (numpy.ndarray) – (n, 2) source coordinates.

• values (numpy.ndarray) – (n,) source values.

• prof (dict) – profile describing the output grid.

• model (str) – variogram model - spherical, exponential, gaussian, linear or power.

• nlags (int) – number of variogram lag bins.

• n_closest (int) – neighbours used per estimate.

• return_variance (bool) – also return the estimation variance grid.

Returns:

the grid, or (grid, variance).

Return type:

numpy.ndarray or tuple

riverarchitect.raster.label_regions(binary, connectivity=4)

Label connected components. Replaces arcpy.sa.RegionGroup.

Parameters:

• binary (numpy.ndarray) – boolean or 0/1 array.

• connectivity (int) – 4 (arcpy’s default) or 8.

Returns:

(labels, count).

Return type:

tuple

riverarchitect.raster.disconnected_mask(binary, connectivity=4, target=None)

Boolean mask of every wetted region that does not reach the main channel.

This is the rule StrandingRisk applies to find pools that have become detached from the main wetted area as discharge drops, and therefore represent fish stranding risk.

Parameters:

• binary (numpy.ndarray) – boolean or 0/1 wetted mask.

• connectivity (int) – 4 (arcpy’s RegionGroup default) or 8.

• target (numpy.ndarray) – boolean mask of the main channel. A region counts as connected when it overlaps it. Without one the largest region is taken as the main channel, which is what it is at a single discharge but not necessarily across a recession - see riverarchitect.stranding.StrandingRisk. A target that no region overlaps is ignored rather than stranding the whole reach.

Returns:

(mask, n_pools).

Return type:

tuple

riverarchitect.raster.within_radius(mask, radius, dx=1.0, dy=1.0)

Cells lying within radius of any True cell of mask.

Replaces the original’s spatial join round trip: raster to points, then SpatialJoin_analysis(..., match_option="CLOSEST", search_radius=r), then points back to raster. That produced exactly this - the set of cells with a source cell inside the radius - and SHArC used it to spread a cover element’s influence over the area it shelters.

Parameters:

• mask (numpy.ndarray) – boolean source cells.

• radius (float) – influence radius in map units.

• dx (float) – cell size, from cell_size().

• dy (float) – cell size, from cell_size().

Returns:

boolean mask of the source cells and everything within the radius.

Return type:

numpy.ndarray

riverarchitect.raster.zonal_statistics(zones, raster_path, stats=('count', 'mean', 'min', 'max', 'sum'), nodata=None)

Statistics of raster_path within each zone polygon.

Replaces arcpy.sa.ZonalStatisticsAsTable.

Pass nodata explicitly for rasters whose declared NoData is an extreme float such as 3.4e38; otherwise those cells contaminate the statistics.

riverarchitect.raster.tabulate_area(array, prof, breaks)

Area per class. Replaces arcpy.sa.TabulateArea.

Returns:

{class_index: area_in_map_units_squared}.

Return type:

dict

riverarchitect.raster.slope(dem, dx, dy, units='DEGREE')

Surface slope. Replaces arcpy.sa.Slope.

Returns degrees by default, or percent rise when units="PERCENT".

riverarchitect.raster.list_rasters(directory, pattern='*.tif')

Sorted list of raster paths in directory. Replaces arcpy.ListRasters().