plans.geo

Standalone geoprocessing routines with minimal dependencies.

Overview

# todo [docstring] – overview Mauris gravida ex quam, in porttitor lacus lobortis vitae. In a lacinia nisl. Pellentesque habitant morbi tristique senectus et netus et malesuada fames ac turpis egestas.

Example

# todo [docstring] – examples Lorem ipsum dolor sit amet, consectetur adipiscing elit. Nulla mollis tincidunt erat eget iaculis. Mauris gravida ex quam, in porttitor lacus lobortis vitae. In a lacinia nisl.

import numpy as np
print("Hello World!")

Mauris gravida ex quam, in porttitor lacus lobortis vitae. In a lacinia nisl. Mauris gravida ex quam, in porttitor lacus lobortis vitae. In a lacinia nisl.

Functions

buffer(grd_input, n_radius)	Calculate a buffer mask
carve_dem(grd_dem, grd_rivers[, ...])	Burn a DEM map with river lines.
classify(array, upvalues, classes)	Classify array based on list of upper values and list of classes values
convert(array, old_values, new_values)	Convert values in array from old to new values.
distance_to_outlet(grd_ldd[, n_res, ...])	Compute the distance to outlet DTO raster of a given basin.
donwscale_parameter(upscaled_value, basin, ...)
downscale_linear(scalar, covariate[, mode, ...])	Applies linear downscaling to a covariate array based on a scalar value.
downscale_parameter_to_units(upscaled_value, ...)
downscale_variance(mean, covariate, scale_factor)	Applies a downscaling operation to a mean array based on a covariance array and a scaling factor.
downscaling_mask(covariate, scale)	Compute a downscaling mask by using a covariate array
downstream_coordinates(n_dir, i, j[, ...])	Compute i and j downstream cell coordinates based on cell flow direction.
euclidean_distance(grd_input)	Calculate the Euclidean distance from pixels with value 1
extents_to_wkt_box(xmin, ymin, xmax, ymax)	Returns a WKT box from the given extents (xmin, ymin, xmax, ymax).
fuzzify(array[, min_value, max_value])	Fuzzify array between min and max values with linear membership function
get_array_bbox(array)	Finds the bounding box for the content (1s) in a 2D pseudo-boolean array.
get_basins_by_gauges(basins, gauges, ...[, ...])	Retrieve a list of upstream basins for a list of gauges.
get_upstream_features(features, field_id, ...)	Retrieves all upstream features from a given starting point.
julian_day(date)	Compute Julian day (day of year).
normalize(array[, min_value, max_value])	Normalize array between min and max values
prune(array[, min_value, max_value])	Inclusive prune values in array
rivers_wedge(grd_rivers[, wedge_width, ...])	Generate a wedge-like trench along the river lines.
shalstab_wetness(flowacc, slope, cellsize, ...)	Calculate the SHALSTAB wetness model
slope(dem, cellsize[, degree])	Calculate slope using gradient-based algorithms
soils(slope, hand[, slope_threshold, ...])	Calculates a basic soil classification grid based on slope and HAND (Height Above Nearest Drainage) data.
solar_altitude(latitude_deg, hour, date)	Compute solar vertical angle (sun elevation or altitude) for a given latitude, date, and local solar time.
solar_azimuth(latitude_deg, hour, date)	Compute solar azimuth angle (degrees clockwise from North).
solar_declination(date)	Solar declination based on Cooper (1969) widely used approximation
solar_hour(hour)	Solar hour angle (0° at solar noon, ±15° per hour)
solar_illumination(year, latitude_deg[, ...])	Simulate hourly solar altitude and azimuth for a full year.
twi(slope, flowacc, cellsize)	Calculate the Topographic Wetness Index TWI
upscale(array[, weights, mode])	Upscales an array using a specified mode.
usle_l(slope, cellsize)	Calculate the USLE L factor of Wischmeier & Smith (1978)
usle_m_a(q, prec, r, k, l, s, c, p[, cellsize])	USLE-M Annual Soil Loss (Kinnell & Risse, 1998)
usle_s(slope)	Calculates the USLE S-factor (slope steepness factor) based on the Wischmeier & Smith (1978) equation.

plans.geo.extents_to_wkt_box(xmin, ymin, xmax, ymax)

Returns a WKT box from the given extents (xmin, ymin, xmax, ymax).

Parameters:

• xmin (float) – Minimum x-coordinate (left)

• ymin – Minimum y-coordinate (bottom)

• xmax (float) – Maximum x-coordinate (right)

• ymax (float) – Maximum y-coordinate (top)

Returns:

A string representing the box in WKT format

Return type:

str

plans.geo.get_basins_by_gauges(basins, gauges, field_gauge, field_basin, field_basin_down, field_geometry='geometry', dissolve=True)

Retrieve a list of upstream basins for a list of gauges.

Parameters:

• basins (geopandas.GeoDataFrame) – A GeoDataFrame containing the basin features.

• gauges (geopandas.GeoDataFrame) – A GeoDataFrame containing the gauge locations.

• field_gauge (str) – The name of the column in gauges that contains the unique identifier for each gauge.

• field_basin (str) – The name of the column in basins that contains the unique identifier for each basin.

• field_basin_down (str) – The name of the column in basins that contains the unique identifier of the downstream basin.

• field_geometry (str) – The name of the geometry column. Default value = “geometry”

• dissolve (bool) – If True, dissolves the upstream basins into a single feature for each gauge. Default value = True

Returns:

A GeoDataFrame with the upstream basins for each gauge, including the gauge identifier. Returns None if no gauges are found within any basin.

Return type:

geopandas.GeoDataFrame

Notes

This function identifies the basins in which each gauge is located and then finds all upstream basins for each of these basins. The resulting GeoDataFrame can optionally be dissolved into a single feature for each gauge.

plans.geo.get_upstream_features(features, field_id, field_id_down, start_id, include_start=True, field_geometry='geometry', is_starting=True)

Retrieves all upstream features from a given starting point.

Parameters:

• features (pandas.DataFrame or geopandas.GeoDataFrame) – A DataFrame containing all features.

• field_id (str) – The name of the column containing the feature’s unique ID.

• field_id_down (str) – The name of the column containing the ID of the downstream feature.

• start_id (int or str) – The ID of the starting feature.

• include_start (bool) – [optional] Whether to include the starting feature in the output. Default value = True

• field_geometry (str) – [optional] The name of the column containing the geometry of the feature. Default value = “geometry”

• is_starting (bool) – [optional] A flag to indicate if it’s the initial call of the function. Default value = True

Returns:

A DataFrame with the upstream features, or None if no upstream features are found.

Return type:

pandas.DataFrame or geopandas.GeoDataFrame or None

plans.geo.get_array_bbox(array)

Finds the bounding box for the content (1s) in a 2D pseudo-boolean array.

Parameters:

array (numpy.ndarray) – 2D numpy array with values 1 (content) and 0 (background)

Returns:

Dictionary with keys “i_min”, “i_max”, “j_min”, “j_max” representing the row and column bounds of the content

Return type:

dict

plans.geo.convert(array, old_values, new_values)

Convert values in array from old to new values.

Parameters:

• array (numpy.ndarray) – Numpy array to convert values

• old_values (numpy.ndarray) – iterable of old values

• new_values (numpy.ndarray) – iterable of new values

Returns:

converted array

Return type:

numpy.ndarray

plans.geo.classify(array, upvalues, classes)

Classify array based on list of upper values and list of classes values

Parameters:

• array – numpy array to classify

• upvalues – 1d numpy array of upper values

• classes – 1d array of classes values

Returns:

numpy array reclassified

plans.geo.normalize(array, min_value=0, max_value=100)

Normalize array between min and max values

Parameters:

• array (float or numpy.ndarray) – Input array or float

• min_value (float) – minimum vale

• max_value (float) – maximum value

Returns:

Output array or float

Return type:

numpy.ndarray

plans.geo.fuzzify(array, min_value=None, max_value=None)

Fuzzify array between min and max values with linear membership function

Parameters:

• array (float or numpy.ndarray) – Input array or float

• min_value (float) – minimum value threshold

• max_value (float) – maximum value threshold

Returns:

Output array or float

Return type:

numpy.ndarray

plans.geo.prune(array, min_value=None, max_value=None)

Inclusive prune values in array

Parameters:

• array (numpy.ndarray) – Numpy array to prune values

• min_value (float) – minimum values to prune (inclusive)

• max_value (float) – maximum values to prude (inclusive)

Returns:

pruned array

Return type:

numpy.ndarray

plans.geo.upscale(array, weights=None, mode='mean')

Upscales an array using a specified mode.

Parameters:

• array (numpy.ndarray) – The input array to be upscaled.

• weights (numpy.ndarray) – [optional] An array of weights to apply. Default value = None

• mode (str) – The mode to use for upscaling, either “mean” or “sum”. Default value = “mean”

Returns:

The upscaled value.

Return type:

int or float

plans.geo.downscale_linear(scalar, covariate, mode='mean', scalar_region=None)

Applies linear downscaling to a covariate array based on a scalar value.

Parameters:

• scalar (float) – The scalar value used for downscaling.

• covariate (numpy.ndarray) – The covariate array to be downscaled.

• scalar_region (numpy.ndarray or None) – Boolean array setting a valid region for the scalar value.

• mode (str) – The aggregation mode to use, “mean” or “sum”. Default value = “mean”

Returns:

The linearly downscaled array.

Return type:

numpy.ndarray

plans.geo.downscale_variance(mean, covariate, scale_factor, reverse=False, mode='simple', max_value=None, min_value=None)

Applies a downscaling operation to a mean array based on a covariance array and a scaling factor.

Parameters:

• mean (numpy.ndarray) – The input mean array.

• covariate (numpy.ndarray) – The covariance array used for downscaling.

• scale_factor (float) – The factor by which to downscale.

• reverse (bool) – Whether to mirror the downscaling effect. Default value = False

• mode (str) – The post-processing mode to apply. Can be “simple”, “prune”, “truncate”, or “compensate”. Default value = “simple”

• max_value (float) – [optional] The maximum value for pruning or truncation.

• min_value (float) – [optional] The minimum value for pruning or truncation.

Returns:

The downscaled array.

Return type:

numpy.ndarray

plans.geo.downscaling_mask(covariate, scale)

Compute a downscaling mask by using a covariate array

Parameters:

• covariate (numpy.ndarray) – Numpy array of a covariate

• scale (float) – scale factor (expected to be positive)

Returns:

downscaling mask array

Return type:

numpy.ndarray

plans.geo.buffer(grd_input, n_radius)

Calculate a buffer mask

Parameters:

• grd_input (numpy.ndarray) – Pseudo-boolean 2D numpy array where pixels with value 1 represent the foreground.

• n_radius (int) – number of pixels for buffer (inclusive)

Returns:

Pseudo-boolean 2D numpy array representing buffer area.

Return type:

numpy.ndarray

plans.geo.euclidean_distance(grd_input)

Calculate the Euclidean distance from pixels with value 1

Parameters:

grd_input (numpy.ndarray) – Pseudo-boolean 2D numpy array where pixels with value 1 represent the foreground.

Returns:

2D numpy array representing the Euclidean distance.

Return type:

numpy.ndarray

Note

• The function uses the distance_transform_edt from scipy.ndimage to compute the Euclidean distance.

• The inputs array is treated as a binary mask, and the distance is calculated from foreground pixels (value 1).

plans.geo.donwscale_parameter(upscaled_value, basin, covariate, mode='mean')

plans.geo.downscale_parameter_to_units(upscaled_value, units, basin, covariate_table, mode='mean')

plans.geo.soils(slope, hand, slope_threshold=25, hand_threshold=5)

Calculates a basic soil classification grid based on slope and HAND (Height Above Nearest Drainage) data.

Parameters:

• slope (numpy.ndarray) – The slope array.

• hand (numpy.ndarray) – The HAND (Height Above Nearest Drainage) array.

• slope_threshold (int) – The slope threshold for identifying neosols. Default value = 25

• hand_threshold (int) – The HAND threshold for identifying alluvial soils. Default value = 5

Returns:

A grid representing soil types (1: default, 2: neosols, 3: alluvial).

Return type:

numpy.ndarray

plans.geo.slope(dem, cellsize, degree=True)

Calculate slope using gradient-based algorithms

Parameters:

• dem (numpy.ndarray) – 2D numpy array representing the Digital Elevation Model (DEM).

• cellsize (float) – The size of a grid cell in the DEM (both in x and y directions).

• degree (bool) – If True (default), the output slope values are in degrees. If False, output units are in radians.

Returns:

numpy array representing the slope values.

Return type:

numpy.ndarray

plans.geo.twi(slope, flowacc, cellsize)

Calculate the Topographic Wetness Index TWI

Parameters:

• slope (numpy.ndarray) – 2D numpy array representing the slope in degrees.

• flowacc (numpy.ndarray) – 2D numpy array representing the flow accumulation.

• cellsize (float) – The size of a grid cell (delta x = delta y).

Returns:

2D numpy array representing the Topographic Wetness Index.

Return type:

numpy.ndarray

Note

• The function uses the formula: TWI = ln( A / cellsize tan(S)), where A is flow accumulation, cellsize is the cell resolution and S is slope in radians.

• The inputs arrays slope and flowacc should have the same dimensions.

• The formula includes a small value (0.01) to prevent issues with tangent calculations for non-NaN values.

plans.geo.rivers_wedge(grd_rivers, wedge_width=3, wedge_depth=3)

Generate a wedge-like trench along the river lines.

Parameters:

• grd_rivers (numpy.ndarray) – Pseudo-boolean grid indicating the presence of rivers.

• wedge_width (int) – Width (single-sided) in pixels. Default is 3.

• wedge_depth (float) – Depth in meters. Default is 3.

Returns:

Grid of the wedge (positive).

Return type:

numpy.ndarray

Notes

• The function generates a wedge-like trench along the river lines based on distance transform.

• The inputs array grd_rivers should be a pseudo-boolean grid where rivers are represented as 1 and others as 0.

• The width w controls the width of the trench, and h controls its height.

plans.geo.carve_dem(grd_dem, grd_rivers, wedge_width=3, wedge_depth=10)

Burn a DEM map with river lines.

Parameters:

• grd_dem (numpy.ndarray) – DEM map.

• grd_rivers (numpy.ndarray) – River map (pseudo-boolean).

• wedge_width (int) – Width parameter in cells. Default is 3.

• wedge_depth (float) – Depth parameter. Default is 10.

Returns:

Burned DEM.

Return type:

numpy.ndarray

Note

• The function burns a DEM map with river lines, creating a wedge-like trench along the rivers.

• The inputs array grd_rivers should be a pseudo-boolean grid where rivers are represented as 1 and others as 0.

• The width w controls the width of the trench, and h controls its height.

plans.geo.downstream_coordinates(n_dir, i, j, ldd_convention='wbt')

Compute i and j downstream cell coordinates based on cell flow direction.

Parameters:

• n_dir (int) – Flow direction code.

• i (int) – i (row) array index.

• j (int) – j (column) array index.

• ldd_convention (str) – String of flow direction convention. Options: wbt, d8, pcraster. Default is wbt.

Returns:

Dictionary of downstream i, j, and distance factor.

Return type:

dict

Notes

• Assumes a specific flow direction convention (wbt, d8, pcraster).

• The output dictionary contains keys ‘i’, ‘j’, and ‘distance’.

• The ‘i’ and ‘j’ values represent downstream cell coordinates.

• The ‘distance’ value is the Euclidean distance to the downstream cell.

LDD conventions:

wbt convention derived from the WhiteboxTool.

wbt convention

64	128	1
32	0	2
16	8	4

d8 convention derived from the SAGA tool.

d8 convention

6	7	8
5	0	1
4	3	2

pcraster convention derived from the PC-Raster tool.

pcraster convention

1	2	3
4	5	6
7	8	9

Examples

>>> downstream_coordinates(n_dir=1, i=3, j=4, ldd_convention='pcraster')
{'i': 4, 'j': 3, 'distance': np.float64(1.4142135623730951)}

Change convention:

>>> downstream_coordinates(n_dir=5, i=10, j=15, ldd_convention='d8')
{'i': 10, 'j': 14, 'distance': np.float64(1.0)}

plans.geo.distance_to_outlet(grd_ldd, n_res=30, ldd_convention='wbt')

Compute the distance to outlet DTO raster of a given basin.

Parameters:

• grd_ldd (numpy.ndarray) – 2d numpy array of flow direction LDD

• n_res (int) – Resolution factor for the output distance raster. Default is 30.

• ldd_convention (str) – String of flow direction convention. Options: wbt, d8, pcraster. Default is wbt.

Returns:

2d numpy array distance

Return type:

numpy.ndarray

plans.geo.shalstab_wetness(flowacc, slope, cellsize, soil_phi, soil_z, soil_c, soil_p, water_p=997, g=9.8, degree=True, kPa=True)

Calculate the SHALSTAB wetness model

Parameters:

• flowacc (numpy.ndarray) – flow accumulation map (square meters)

• slope (numpy.ndarray) – slope map (degrees or radians)

• cellsize (float) – grid cell size (meters)

• soil_phi (numpy.ndarray or float) – soil angle of internal friction (degrees or radians)

• soil_z (numpy.ndarray or float) – soil depth (meters)

• soil_c (numpy.ndarray or float) – soil cohesion (Pa or N/m0² or kg/ms²)

• soil_p (numpy.ndarray or float) – soil density (kg/m0³)

• water_p (float) – water density (kg/m0³)

• g (float) – gravity acceleration (m0 / s²)

• degree (bool) – flag to note if slope and soil_phi are in degrees

• kPa (bool) – flag to note if soil_c in kPa

Returns:

map of q/T ratio

Return type:

numpy.ndarray

plans.geo.usle_l(slope, cellsize)

Calculate the USLE L factor of Wischmeier & Smith (1978)

# todo [docstring] improve notes sphinx design L = (x / 22.13) ^ m0

where:

m0 = 0.2 when sinθ < 0.01; m0 = 0.3 when 0.01 ≤ sinθ ≤ 0.03; m0 = 0.4 when 0.03 < sinθ < 0.05; m0 = 0.5 when sinθ ≥ 0.05

x is the plot lenght taken as 1.4142 * cellsize (diagonal length of cell)

Parameters:

• slope (numpy.ndarray) – slope in degrees of terrain 2d array

• cellsize (float) – cell size in meters

Returns:

Wischmeier & Smith (1978) L factor 2d array

Return type:

numpy.ndarray

Notes

This implementation follows the slope length (L) factor formulation proposed by Wischmeier & Smith (1978) as part of the Universal Soil Loss Equation (USLE). The L factor accounts for the effect of slope length on erosion by scaling soil loss relative to a standard plot length of 22.13 m.

The slope length factor is computed as:

\[L = {\left( \frac{x}{22.13} \right)}^{m}\]

where:

• \(x\) is the effective slope length (m), here approximated as \(x = 1.4142 \times \text{cellsize}\), corresponding to the diagonal length of a raster cell.

• \(m\) is an empirical exponent dependent on slope steepness, defined as a function of \(\sin \theta\).

The exponent \(m\) is assigned according to the following criteria:

• \(m = 0.2\) when \(\sin \theta < 0.01\)

• \(m = 0.3\) when \(0.01 \le \sin \theta \le 0.03\)

• \(m = 0.4\) when \(0.03 < \sin \theta < 0.05\)

• \(m = 0.5\) when \(\sin \theta \ge 0.05\)

This formulation is widely used in raster-based implementations of USLE where slope length is approximated at the grid-cell scale.

References

Wischmeier, W. H., & Smith, D. D. (1978). Predicting rainfall erosion losses: A guide to conservation planning. USDA Agriculture Handbook No. 537, U.S. Department of Agriculture, Washington, DC.

plans.geo.usle_s(slope)

Calculates the USLE S-factor (slope steepness factor) based on the Wischmeier & Smith (1978) equation.

Parameters:

slope (numpy.ndarray) – Slope in degrees of terrain.

Returns:

The USLE S-factor.

Return type:

numpy.ndarray

Notes

This implementation follows the slope steepness (S) factor formulation proposed by Wischmeier & Smith (1978), where the slope angle is expressed in degrees and internally converted to radians prior to evaluation.

The factor is defined as a function of the sine of the slope angle (\(\theta\)), making it suitable for raster-based terrain analysis where slope is commonly derived in angular units.

\[S = 65.41 \, (\sin \theta)^2 + 4.56 \, \sin \theta + 0.065\]

where:

• \(\theta\) is the slope angle in radians.

• \(S\) is a dimensionless slope steepness factor used in the Universal Soil Loss Equation (USLE).

References

Wischmeier, W. H., & Smith, D. D. (1978). Predicting rainfall erosion losses: A guide to conservation planning. USDA Agriculture Handbook No. 537, U.S. Department of Agriculture, Washington, DC.

plans.geo.usle_m_a(q, prec, r, k, l, s, c, p, cellsize=30)

USLE-M Annual Soil Loss (Kinnell & Risse, 1998)

# todo [docstring] improve notes sphinx design

Parameters:

• q (numpy.ndarray) – Annual runoff in mm/year.

• prec (numpy.ndarray or float) – Annual precipitation in mm/year.

• r (numpy.ndarray or float) – Rain erosivity in MJ mm h-1 ha-1 year-1.

• k (numpy.ndarray or float) – Erodibility K factor in ton h MJ-1 mm-1.

• l (numpy.ndarray or float) – USLE/RUSLE L factor.

• s (numpy.ndarray or float) – USLE/RUSLE S factor.

• c (numpy.ndarray or float) – C_UM factor.

• p (numpy.ndarray or float) – USLE/RUSLE P factor.

• cellsize (float) – Grid cell size in meters. Default value = 30

Returns:

Annual Soil Loss in ton/year.

Return type:

numpy.ndarray

plans.geo.julian_day(date)

Compute Julian day (day of year).

Parameters:

date (str) – Calendar date in YYYY-MM-DD format

Returns:

Julian day

Return type:

int

plans.geo.solar_declination(date)

Solar declination based on Cooper (1969) widely used approximation

Parameters:

date (str) – Calendar date in YYYY-MM-DD format

Returns:

Solar declination angle in degrees

Return type:

float

plans.geo.solar_hour(hour)

Solar hour angle (0° at solar noon, ±15° per hour)

Parameters:

hour (float) – Local solar time [0–24]

Returns:

Hour angle in degrees

Return type:

float

plans.geo.solar_altitude(latitude_deg, hour, date)

Compute solar vertical angle (sun elevation or altitude) for a given latitude, date, and local solar time.

Parameters:

• latitude_deg (float) – Latitude in degrees (positive north, negative south)

• hour (float) – Local solar time [0–24]

• date (str) – Calendar date in YYYY-MM-DD format

Returns:

Solar vertical angle in degrees above the horizon (aka solar altitude)

Return type:

float

plans.geo.solar_azimuth(latitude_deg, hour, date)

Compute solar azimuth angle (degrees clockwise from North).

Parameters:

• latitude_deg (float) – Latitude in degrees (positive north, negative south)

• hour (float) – Local solar time [0–24]

• date (str) – Calendar date in YYYY-MM-DD format

Returns:

Solar azimuth angle in degrees

Return type:

float

plans.geo.solar_illumination(year, latitude_deg, frequency='1h')

Simulate hourly solar altitude and azimuth for a full year.

Parameters:

• year (int) – Year to simulate

• latitude_deg (float) – Latitude in degrees

• frequency (str) – Frequency flag for step simulation based on pandas

Returns:

Simulation data

Return type:

pandas.`DataFrame`