Preprocessing

Catchment topography

class hydrobricks.preprocessing.catchment_topography.CatchmentTopography(catchment: Catchment)

Bases: object

A class to represent the topography of a catchment.

calculate_slope_aspect() → None

Calculate the slope and aspect of the whole DEM.

Computes slope in degrees and aspect (0-360°) from the DEM using spatial derivatives. Results stored in self.slope and self.aspect.

extract_unit_max_elevation(mask_unit: numpy.ndarray) → float

Extract the maximum elevation for a hydro unit.

Parameters:

mask_unit – Boolean mask array identifying the cells of the hydro unit.

Returns:

Maximum elevation in meters, rounded to 2 decimal places.

Return type:

float

extract_unit_mean_aspect(mask_unit: numpy.ndarray) → float

Extract the circular mean aspect for a hydro unit.

Computes the circular mean of aspect values within a unit mask, accounting for the circular nature of aspect angles (0° and 360° are equivalent).

Parameters:

mask_unit – Boolean mask array identifying the cells of the hydro unit.

Returns:

Mean aspect in degrees (0-360).

Return type:

float

extract_unit_mean_elevation(mask_unit: numpy.ndarray) → float

Extract the mean elevation for a hydro unit.

Parameters:

mask_unit – Boolean mask array identifying the cells of the hydro unit.

Returns:

Mean elevation in meters, rounded to 2 decimal places.

Return type:

float

extract_unit_mean_slope(mask_unit: numpy.ndarray) → float

Extract the mean slope for a hydro unit.

Parameters:

mask_unit – Boolean mask array identifying the cells of the hydro unit.

Returns:

Mean slope in degrees, rounded to 2 decimal places.

Return type:

float

extract_unit_min_elevation(mask_unit: numpy.ndarray) → float

Extract the minimum elevation for a hydro unit.

Parameters:

mask_unit – Boolean mask array identifying the cells of the hydro unit.

Returns:

Minimum elevation in meters, rounded to 2 decimal places.

Return type:

float

get_hillshade(azimuth: float = 315, altitude: float = 45, z_factor: float = 1) → numpy.ndarray

Create a hillshade from the DEM.

Adapted from https://github.com/royalosyin/Work-with-DEM-data-using-Python- from-Simple-to-Complicated/blob/master/ex07-Hillshade%20from%20a%20Digital %20Elevation%20Model%20(DEM).ipynb

Parameters:

• azimuth – The desired azimuth for the hillshade.

• altitude – The desired sun angle altitude for the hillshade.

• z_factor – The z factor to amplify the relief.

Returns:

2D array of hillshade values (0-255) representing shaded relief visualization.

Return type:

np.ndarray

get_mean_elevation() → float

Get the catchment mean elevation.

Computes the mean elevation across all DEM cells within the catchment, ignoring NaN values.

Returns:

The catchment mean elevation in meters.

Return type:

float

resample_dem_and_calculate_slope_aspect(resolution: float, output_path: str | Path | None = None) → tuple[rasterio.DatasetReader, np.ndarray, np.ndarray, np.ndarray]

Resample the DEM and calculate the slope and aspect of the whole DEM.

Downsamples the DEM to a specified resolution using bilinear interpolation, then computes slope and aspect for the resampled DEM. Results are saved to disk and can be used for coarser resolution calculations.

Parameters:

• resolution – Desired pixel resolution in meters. If None or matches current resolution, uses the original DEM without resampling.

• output_path – Path of the directory to save the downsampled DEM to. Required if resolution differs from current DEM resolution.

Returns:

Tuple containing: - Resampled DEM dataset - Masked DEM data array - Slope array (in degrees) - Aspect array (in degrees, 0-360)

Return type:

tuple[rasterio.DatasetReader, np.ndarray, np.ndarray, np.ndarray]

Catchment discretization

class hydrobricks.preprocessing.catchment_discretization.CatchmentDiscretization(catchment: Catchment)

Bases: object

Class to handle the discretization of catchments.

create_elevation_bands(method: str = 'equal_intervals', number: int = 100, distance: int = 50, min_elevation: int | None = None, max_elevation: int | None = None) → None

Construction of the elevation bands based on the chosen method.

Creates hydro units based on elevation bands using either equal-interval contours or quantile-based discretization. Results are stored in the catchment’s map_unit_ids and hydro_units objects.

Parameters:

• method – The method to build the elevation bands: ‘equal_intervals’ = fixed contour intervals (needs ‘distance’ parameter) ‘quantiles’ = quantiles of the catchment area (same surface; needs ‘number’ parameter)

• number – Number of bands to create for the ‘quantiles’ method.

• distance – Distance (m) between the contour lines for the ‘equal_intervals’ method.

• min_elevation – Minimum elevation of the elevation bands (to homogenize between runs).

• max_elevation – Maximum elevation of the elevation bands (to homogenize between runs).

discretize_by(criteria: str, elevation_method: str = 'equal_intervals', elevation_number: int = 100, elevation_distance: int = 100, min_elevation: int | None = None, max_elevation: int | None = None, slope_method: str = 'equal_intervals', slope_number: int = 6, slope_distance: int = 15, min_slope: int = 0, max_slope: int = 90, radiation_method: str = 'equal_intervals', radiation_number: int = 5, radiation_distance: int = 50, min_radiation: int | None = None, max_radiation: int | None = None) → None

Construction of the elevation bands based on the chosen method.

Discretizes the catchment into hydro units based on multiple criteria (elevation, slope, aspect, radiation). Creates all combinations of specified criteria and populates the map_unit_ids and hydro_units objects.

Parameters:

• criteria – The criteria to use to discretize the catchment (can be combined): ‘elevation’ = elevation bands ‘aspect’ = aspect according to the cardinal directions (4 classes) ‘radiation’ = potential radiation (Hock, 1999) ‘slope’ = slope in degrees

• elevation_method – The method to build the elevation bands: ‘equal_intervals’ = fixed contour intervals (provide the ‘elevation_distance’ parameter) ‘quantiles’ = quantiles of the catchment area (same surface; provide the ‘elevation_number’ parameter)

• elevation_number – Number of elevation bands to create for the ‘quantiles’ method.

• elevation_distance – Distance (m) between the contour lines for the ‘equal_intervals’ method.

• min_elevation – Minimum elevation of the elevation bands (to homogenize between runs).

• max_elevation – Maximum elevation of the elevation bands (to homogenize between runs).

• slope_method – The method to build the slope categories: ‘equal_intervals’ = fixed slope intervals (needs ‘slope_distance’ parameter) ‘quantiles’ = quantiles of the catchment area (same surface; provide the ‘slope_number’ parameter)

• slope_number – Number of slope bands to create for the ‘quantiles’ method.

• slope_distance – Distance (degrees) between the slope lines for the ‘equal_intervals’ method.

• min_slope – Minimum slope of the slope bands (to homogenize between runs).

• max_slope – Maximum slope of the slope bands (to homogenize between runs).

• radiation_method – The method to build the radiation categories: ‘equal_intervals’ = fixed radiation intervals (provide the ‘radiation_distance’ parameter) ‘quantiles’ = quantiles of the catchment area (same surface; provide the ‘radiation_number’ parameter)

• radiation_number – Number of radiation bands to create when using the ‘quantiles’ method.

• radiation_distance – Distance (W/m2) in terms of radiation for the ‘equal_intervals’ method.

• min_radiation – Minimum radiation of the radiation bands (to homogenize between runs).

• max_radiation – Maximum radiation of the radiation bands (to homogenize between runs).

Catchment connectivity

class hydrobricks.preprocessing.catchment_connectivity.CatchmentConnectivity(catchment: Catchment)

Bases: object

Class to handle connectivity of catchments.

calculate(mode: str = 'multiple', force_connectivity: bool = False, precision: int = 3) → pandas.DataFrame

Calculate the connectivity between hydro units using a flow accumulation partition. Connectivity between hydro units is forced to stay within the catchment.

Computes how water flows from each hydro unit to downstream units based on flow direction and accumulation. Results in a DataFrame with connectivity values representing the fraction of flow directed to each downstream unit.

Parameters:

• mode – The mode to calculate the connectivity: ‘single’ = keep the highest connectivity only ‘multiple’ = keep all connectivity values (multiple connections)

• force_connectivity – If True, all hydro units are forced to have a connection to another hydro unit. When there is no downstream hydro unit, the connectivity is set to the neighboring hydro units, proportionally to the length of the common border. If False, and if a hydro unit contributes mostly to surfaces out of the catchment, the connectivity will be nulled. Default is False (recommended).

• precision – The precision of the connectivity values. Default is 3. This is used to round the connectivity values to a given number of decimal places.

Returns:

DataFrame with connectivity information for each hydro unit.

Return type:

pd.DataFrame

Glacier evolution — delta-h

class hydrobricks.preprocessing.glacier_evolution_delta_h.GlacierEvolutionDeltaH(hydro_units: HydroUnits | None = None)

Bases: object

Class for glacier mass balance. It is based on the following papers:

• Seibert, J., Vis, M. J. P., Kohn, I., Weiler, M., and Stahl, K.: Technical note: Representing glacier geometry changes in a semi-distributed hydrological model, Hydrol. Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018, 2018.

• Huss, M., Jouvet, G., Farinotti, D., and Bauder, A.: Future high-mountain hydrology: a new parameterization of glacier retreat, Hydrol. Earth Syst. Sci., 14, 815–829, https://doi.org/10.5194/hess-14-815-2010, 2010.

compute_ice_thickness_change(catchment: Catchment, ice_thickness_old: str, ice_thickness_new: str, elevation_bands_distance: int = 10, smooth_window: int | None = 5, nodata: float | None = None) → pd.DataFrame

Compute the normalized ice thickness change between two glacier ice thickness rasters for the delta-h method.

Computes elevation-dependent ice thickness changes and normalizes them for use in the delta-h parameterization. Results are stored in self.norm_elevations_bands and self.norm_delta_we_bands for later use.

Parameters:

• catchment – The catchment object.

• ice_thickness_old – Path to the TIF file containing the old glacier thickness in meters.

• ice_thickness_new – Path to the TIF file containing the new glacier thickness in meters.

• elevation_bands_distance – Distance between elevation bands in meters. Default is 10 m.

• smooth_window – Optional integer window for 1D smoothing (pandas rolling) applied to the aggregated curve.

• nodata – Value to treat as nodata in inputs (replaced with np.nan).

Returns:

The glacier_change_df DataFrame containing the normalized ice thickness change data.

Return type:

pd.DataFrame

compute_initial_ice_thickness(catchment: Catchment, glacier_outline: str | None = None, ice_thickness: str | None = None, elevation_bands_distance: int = 10, pixel_based_approach: bool = True) → pd.DataFrame

Extract the initial ice thickness to be used in compute_lookup_table() for the glacier mass balance calculation.

Parameters:

• catchment – The catchment object.

• glacier_outline – Path to the SHP file containing the glacier extents. If None, the glacier extents are extracted from the ice thickness geotiff. Either this or ice_thickness should be provided.

• ice_thickness – Path to the TIF file containing the glacier thickness in meters. If None, the ice thickness is estimated based on the glacier area using the Bahr et al. (1997) formula. Either this or glacier_outline should be provided.

• elevation_bands_distance – Distance between elevation bands in meters. Default is 10 m.

• pixel_based_approach – Whether to compute the glacier area evolution from the topography. If True, the glacier area is computed from the topography and the ice thicknesses are updated accordingly. Otherwise, the glacier area is updated based on the Bahr et al. (1997) formula. Default is True.

Return type:

The glacier_df DataFrame containing the glacier data.

compute_lookup_table(catchment: Catchment | None = None, glacier_profile_csv: str | None = None, glacier_df: pd.DataFrame | None = None, observed_dh: pd.DataFrame | None = None, nb_increments: int = 200, update_width: bool = True, update_width_reference: str = 'initial')

Prepare the glacier mass balance lookup table. The glacier mass balance is calculated using the delta-h method (Huss et al., 2010) using the lookup table method (Seibert et al., 2018). The glacier mass balance is calculated for each elevation band and melt increment and is synthesized in a lookup table per hydro unit.

Parameters:

• catchment – The catchment object. Needed when pixel_based_approach is True.

• glacier_profile_csv –

  Path to the CSV file containing glacier data. An elevation band is smaller than a hydro unit, usually around 10 m high. It should contain the following columns: - elevation: Elevation (m) for each elevation band i (e.g. every 10m). - glacier_area: Glacier area (m2) for each elevation band i. - glacier_thickness: Glacier thickness (m) for each elevation band i. - hydro_unit_id: Hydro unit ID for each elevation band i.

  Format example: elevation glacier_area glacier_thickness hydro_unit_id m m2 m - 1670 0 0 2 1680 0 0 2 1690 2500 14.7 2 1700 9375 23.2 3 1710 11875 25.2 3 1720 9375 23.2 3 1730 9375 23.2 3 … … … …

  If not provided, the glacier data is assumed to be already loaded in the glacier_df attribute or to be passed as a DataFrame.

• glacier_df – DataFrame containing glacier data. Alternative to glacier_profile_csv.

• observed_dh – Optional dataframe containing observed delta-h values.

• nb_increments – Number of increments for glacier mass balance calculation. Default is 200.

• update_width – Whether to update the glacier width at each iteration (Eq. 7 Seibert et al., 2018). Default is True. Ignored if pixel_based_approach is True.

• update_width_reference – Reference for updating the glacier width (Eq. 7 Seibert et al., 2018). Default is ‘initial’. Options are: - ‘initial’: Use the initial glacier width. - ‘previous’: Use the glacier width from the previous iteration. Ignored if pixel_based_approach is True.

get_lookup_table_area() → pandas.DataFrame

Get the glacier area evolution lookup table.

Return type:

The glacier area evolution lookup table.

get_lookup_table_volume() → pandas.DataFrame

Get the glacier volume evolution lookup table.

Return type:

The glacier volume evolution lookup table.

save_as_csv(output_dir: str | Path)

Save the glacier evolution lookup table as a CSV file.

Parameters:

output_dir – Path to the directory where the CSV file will be saved.

Glacier evolution — area scaling

class hydrobricks.preprocessing.glacier_evolution_area_scaling.GlacierEvolutionAreaScaling

Bases: object

Class for a simple volume-area scaling approach for glacier evolution. It computes the area per hydro unit (e.g. elevation band) for an incremental melt per band. It does not account for the glacier dynamic and should not be used for medium to large glaciers. The delta-h method is preferred for such cases.

compute_lookup_table(catchment: Catchment, ice_thickness: str | Path, nb_increments: int = 200) → None

Extract the initial ice thickness to be used in compute_lookup_table() for the glacier mass balance calculation.

Computes glacier area and volume evolution lookup tables based on incremental melting using a simple volume-area scaling approach. Results stored in lookup_table_area and lookup_table_volume attributes.

Parameters:

• catchment – The catchment object.

• ice_thickness – Path to the TIF file containing the glacier thickness in meters.

• nb_increments – Number of increments for glacier mass balance calculation. Default is 200.

get_lookup_table_area() → pandas.DataFrame

Get the glacier area evolution lookup table.

Return type:

The glacier area evolution lookup table.

get_lookup_table_volume() → pandas.DataFrame

Get the glacier volume evolution lookup table.

Return type:

The glacier volume evolution lookup table.

save_as_csv(output_dir: str | Path) → None

Save the glacier evolution lookup table as a CSV file.

Writes three CSV files: glacier area evolution, glacier volume evolution, and optional detailed glacier area and water equivalent evolution per pixel.

Parameters:

output_dir – Path to the directory where the CSV file will be saved.

Potential solar radiation

class hydrobricks.preprocessing.potential_solar_radiation.PotentialSolarRadiation(catchment: Catchment)

Bases: object

A class to handle solar radiation data for a catchment area.

static calculate_cast_shadows(dem_dataset: rasterio.Dataset, masked_dem: np.ndarray, zenith: float, azimuth: float, xv: np.ndarray = None, yv: np.ndarray = None) → np.ndarray

Calculate the cast shadows on terrain from a given sun position.

The approach relies on tilting the DEM to obtain a horizon matching the sun rays and filling the DEM. The algorithm is applied to the whole topography before masking the areas outside the catchment.

Parameters:

• dem_dataset – Opened rasterio DEM dataset.

• masked_dem – DEM array masked for catchment extent.

• zenith – Solar zenith angle in degrees.

• azimuth – Solar azimuth angle in degrees from north.

• xv – Pre-computed x grid for the DEM, to avoid repeated allocation. If None, it will be computed inside the function. Default is None.

• yv – Pre-computed y grid for the DEM, to avoid repeated allocation. If None, it will be computed inside the function. Default is None.

Returns:

Binary mask where 1 indicates shadowed areas and 0 indicates sunlit areas.

Return type:

np.ndarray

Notes

The algorithm is applied to the whole topography before masking the areas outside the catchment.

calculate_daily_potential_radiation(output_path: str, resolution: float | None = None, atmos_transmissivity: float = 0.75, steps_per_hour: int = 4, with_cast_shadows: bool = True) → None

Compute the daily mean potential clear-sky direct solar radiation at the DEM surface [W/m²] using Hock (1999)’s equation. It is computed for each day of the year and saved in a netcdf file.

Computes radiation accounting for terrain slope, aspect, latitude, and solar declination. Optionally includes cast shadows from surrounding terrain.

Parameters:

• output_path – Path to for created daily and annual mean potential clear-sky direct solar radiation files.

• resolution – Desired pixel resolution, default is the DEM resolution.

• atmos_transmissivity – Mean clear-sky atmospheric transmissivity, default is 0.75 (value taken in Hock 1999)

• steps_per_hour – Number of steps per hour to compute the potential radiation, default is 4.

• with_cast_shadows – If True, the cast shadows are taken into account. Default is True.

Notes

This function is based on the R package TopoSol, authored by Matthew Olson.

References

• Original R package: https://github.com/mattols/TopoSol

• Hock, R. (1999). A distributed temperature-index ice-and snowmelt model including potential direct solar radiation. J. Glaciol., 45(149), 101-111.

static get_solar_azimuth_to_north(hour_angles: float | np.ndarray, lat_rad: float, solar_declination: float) → float | np.ndarray

Compute the solar azimuth relative to the north. See get_solar_azimuth_to_south() for more details.

static get_solar_azimuth_to_south(hour_angles: float | np.ndarray, lat_rad: float, solar_declination: float) → np.ndarray

Compute the solar azimuth relative to the south.

The solar azimuth is the angle between the sun and the observer’s meridian. It is typically measured in degrees. Azimuth with negative values before solar noon and positive values with positive values after solar noon. Solar noon is defined by the change in sign of the hour angle (negative in the morning, positive in the afternoon). From https://www.astrolabe-science.fr/diagramme-solaire-azimut-hauteur

Parameters:

• hour_angles – Array with the hour angles.

• lat_rad – Latitude in radians.

• solar_declination – Solar declination.

Return type:

The solar azimuth in degrees.

static get_solar_declination_rad(day_of_year: int | np.ndarray) → float

Compute the solar declination.

The solar declination is the angle between the rays of the Sun and the plane of the Earth’s equator. It represents how much the sun is tilted towards or away from the observer’s latitude. The calculation involves trigonometric functions to account for the observer’s latitude and the position of the sun in the sky.

Parameters:

day_of_year – Day of the year (1-366).

Return type:

The solar declination in radians

static get_solar_hour_angle_limit(solar_declination: float | np.ndarray, lat_rad: float) → float | np.ndarray

Compute the hour angle limit value (min/max).

The hour angle is the angular distance between the sun and the observer’s meridian. It is typically measured in degrees. The tangent of solar_declin/lat_rad represents the ratio of the opposite side (vertical component of the sun’s rays) to the adjacent side (horizontal component). It helps capture how much the sun/the observer’s location is tilted north or south relative to the equator.

Parameters:

• solar_declination – Solar declination in radians.

• lat_rad – Latitude in radians.

Return type:

The limit value of the hour angle in radians.

static get_solar_zenith(hour_angles: float | np.ndarray, lat_rad: float, solar_declination: float) → float | np.ndarray

Compute the solar zenith.

The Solar zenith (IQBAL 2012) is the angle between the sun and the vertical (zenith) position directly above the observer. The result is expressed in degrees. The calculation involves trigonometric functions to account for the observer’s latitude, the solar declination, and the position of the sun in the sky (represented by the Hour Angles).

Parameters:

• hour_angles – Hour angle(s).

• lat_rad – Latitude in radians.

• solar_declination – Solar declination in radians.

Return type:

The solar zenith in degrees.

load_mean_annual_radiation_raster(dir_path: str, filename: str = 'annual_potential_radiation.tif') → None

Load the mean annual radiation raster.

Loads a pre-computed mean annual radiation raster from a file and stores it in the mean_annual_radiation attribute.

Parameters:

• dir_path – Path to the input file directory.

• filename – Name of the input file. Default is ‘annual_potential_radiation.tif’.

upscale_and_save_mean_annual_radiation_rasters(mean_annual_radiation: np.ndarray, dem: rasterio.Dataset, output_path: str, output_filename: str = 'annual_potential_radiation.tif') → None

Save the mean annual radiation rasters (downsampled and at DEM resolution) to a file.

Upscales radiation from a coarser resolution to the DEM resolution using bilinear interpolation if necessary, and saves both the intermediate downsampled raster and the final raster at DEM resolution.

Parameters:

• mean_annual_radiation – Downsampled mean annual radiation.

• dem – The DEM at the resolution of the radiation.

• output_path – Path to the output directory.

• output_filename – Name of the output file. Default is ‘annual_potential_radiation.tif’.