diff --git a/docs/images/terrain_agreement/1.png b/docs/images/terrain_agreement/1.png new file mode 100644 index 00000000..1c6f2b90 Binary files /dev/null and b/docs/images/terrain_agreement/1.png differ diff --git a/docs/images/terrain_agreement/2.png b/docs/images/terrain_agreement/2.png new file mode 100644 index 00000000..b79d2f31 Binary files /dev/null and b/docs/images/terrain_agreement/2.png differ diff --git a/docs/images/terrain_agreement/3.png b/docs/images/terrain_agreement/3.png new file mode 100644 index 00000000..e50ec782 Binary files /dev/null and b/docs/images/terrain_agreement/3.png differ diff --git a/docs/images/terrain_agreement/4.png b/docs/images/terrain_agreement/4.png new file mode 100644 index 00000000..c7f6e09b Binary files /dev/null and b/docs/images/terrain_agreement/4.png differ diff --git a/docs/images/terrain_agreement/5.png b/docs/images/terrain_agreement/5.png new file mode 100644 index 00000000..35451caa Binary files /dev/null and b/docs/images/terrain_agreement/5.png differ diff --git a/docs/images/terrain_agreement/6.png b/docs/images/terrain_agreement/6.png new file mode 100644 index 00000000..f588caf4 Binary files /dev/null and b/docs/images/terrain_agreement/6.png differ diff --git a/docs/images/terrain_agreement/7.png b/docs/images/terrain_agreement/7.png new file mode 100644 index 00000000..85954b71 Binary files /dev/null and b/docs/images/terrain_agreement/7.png differ diff --git a/docs/images/terrain_agreement/8.png b/docs/images/terrain_agreement/8.png new file mode 100644 index 00000000..b2772d1b Binary files /dev/null and b/docs/images/terrain_agreement/8.png differ diff --git a/docs/source/images/terrain_agreement/1.png b/docs/source/images/terrain_agreement/1.png new file mode 100644 index 00000000..1c6f2b90 Binary files /dev/null and b/docs/source/images/terrain_agreement/1.png differ diff --git a/docs/source/images/terrain_agreement/2.png b/docs/source/images/terrain_agreement/2.png new file mode 100644 index 00000000..b79d2f31 Binary files /dev/null and b/docs/source/images/terrain_agreement/2.png differ diff --git a/docs/source/images/terrain_agreement/3.png b/docs/source/images/terrain_agreement/3.png new file mode 100644 index 00000000..e50ec782 Binary files /dev/null and b/docs/source/images/terrain_agreement/3.png differ diff --git a/docs/source/images/terrain_agreement/4.png b/docs/source/images/terrain_agreement/4.png new file mode 100644 index 00000000..c7f6e09b Binary files /dev/null and b/docs/source/images/terrain_agreement/4.png differ diff --git a/docs/source/images/terrain_agreement/5.png b/docs/source/images/terrain_agreement/5.png new file mode 100644 index 00000000..35451caa Binary files /dev/null and b/docs/source/images/terrain_agreement/5.png differ diff --git a/docs/source/images/terrain_agreement/6.png b/docs/source/images/terrain_agreement/6.png new file mode 100644 index 00000000..f588caf4 Binary files /dev/null and b/docs/source/images/terrain_agreement/6.png differ diff --git a/docs/source/images/terrain_agreement/7.png b/docs/source/images/terrain_agreement/7.png new file mode 100644 index 00000000..85954b71 Binary files /dev/null and b/docs/source/images/terrain_agreement/7.png differ diff --git a/docs/source/images/terrain_agreement/8.png b/docs/source/images/terrain_agreement/8.png new file mode 100644 index 00000000..b2772d1b Binary files /dev/null and b/docs/source/images/terrain_agreement/8.png differ diff --git a/docs/source/tech_summary.rst b/docs/source/tech_summary.rst index 49c0aac3..280f00c2 100644 --- a/docs/source/tech_summary.rst +++ b/docs/source/tech_summary.rst @@ -59,7 +59,7 @@ a geopackage file. 2 - Model and NWM network conflation ------------------------------------ -(relevant endpoints: :doc:`conflate_model `, +(relevant endpoints: :doc:`conflate_model `, :doc:`compute_conflation_metrics `) .. image:: images/source_w_nwm.png @@ -95,7 +95,7 @@ definitions for each of the JSON fields are provided below. are measured at each HEC-RAS cross-section and summary statistics are reported in the conflation metrics output. - * **centerline_offset** measures the straightline distance between RAS centerline + * **centerline_offset** measures the straightline distance between RAS centerline and NWM reach line * **thalweg_offset** measures the straightline distance between lowest point @@ -114,9 +114,9 @@ definitions for each of the JSON fields are provided below. * **network** is the distance along the NWM reach between upstream and downstream cross-section - + * **network_to_ras_ratio** is the network length divided by ras length - + .. image:: images/length_metrics.png :width: 400 :alt: Length conflation metrics @@ -139,7 +139,7 @@ definitions for each of the JSON fields are provided below. 3 - Sub model creation ---------------------- -(relevant endpoints: :doc:`extract_submodel `, +(relevant endpoints: :doc:`extract_submodel `, :doc:`create_ras_terrain `) Once NWM reaches have been associated with relevant parts of the HEC-RAS model, @@ -152,18 +152,105 @@ from any virtual raster source, but by default, ripple1d will download a `1/3 arcsecond DEM from USGS `_ +As part of terrain generation, a suite of metrics are generated to quantify the agreement of the newly generated DEM terrain and the source model cross-section geometry. Metrics are first generated for each cross-section at a set of water surface elevations ranging from the section invert to the lower of the two source model section endpoints. All metrics (except for residual summary statistics) are aggregated to the cross-section level by averaging across all measured stages. Another set of shape metrics as well as residual summary statistics are computed for the whole cross-section. All cross-section metrics (except for residual summary statistics) are aggregated to the model level by averaging across all cross-sections. + +**Example Cross-Sections and Their Metrics** + +**Perfectly Aligned** + +.. image:: images/terrain_agreement/1.png + :width: 1000 + :alt: Perfectly aligned data + :align: center + +**Noisy** + +.. image:: images/terrain_agreement/2.png + :width: 1000 + :alt: Noisy data + :align: center + +**Vertically Offset** + +.. image:: images/terrain_agreement/3.png + :width: 1000 + :alt: Vertically offset data + :align: center + +**Horizontally Offset** + +.. image:: images/terrain_agreement/4.png + :width: 1000 + :alt: Horizontally offset data + :align: center + +**Squeezed** + +.. image:: images/terrain_agreement/5.png + :width: 1000 + :alt: Squeezed data + :align: center + +**Truncated** + +.. image:: images/terrain_agreement/6.png + :width: 1000 + :alt: Truncated data + :align: center + +**Low Fidelity** + +.. image:: images/terrain_agreement/7.png + :width: 1000 + :alt: Low-fidelity data + :align: center + +**Complete Misalignment** + +.. image:: images/terrain_agreement/8.png + :width: 1000 + :alt: Completely misaligned data + :align: center + +**Metric Descriptions and Interpretations** + +* **Residual Summary Statistics** These statistics summarize the difference between source model and DEM elevations at each cross-section vertex. These metrics can be used to assess the magnitude of difference between the two sections, however, since they are not scaled, acceptable ranges will vary from river to river. (Note: normalized RMSE is RMSE divided by the interquartile range and attempts to be a scaled error metric) + +* **Inundation Overlap** The intersection of the wetted top widths divided by the union of the wetted top widths (closer to 1 is better). This metric can be used to determine spatially explicit agreement of inundation. A good example is shown in the horizontally offset example above. + +* **Top-Width Agreement** Calculated as one minus the symmetric mean absolute percentage error (sMAPE) of the source model wetted top-width and the DEM wetted top-width (closer to 1 is better). This metric is a non-spatially explicit version of inundation overlap. A good example is shown in the horizontally offset example above as well as the squeezed example. + +* **Flow Area Overlap** The intersection of the flow areas divided by the union of the flow areas (closer to 1 is better). This metric can be used to determine spatially explicit agreement of the cross-section area. A good example is shown in the horizontally offset example above. + +* **Flow Area Agreement** Calculated as one minus the sMAPE of the source model flow area and the DEM flow area (closer to 1 is better). This metric is a non-spatially explicit version of flow area overlap. A good example is shown in the horizontally offset example above as well as the squeezed example. + +* **Hydraulic Radius Agreement** Calculated as one minus the sMAPE of the source model hydraulic radius and the DEM hydarulic radius (closer to 1 is better). This metric captures some of how well the hydarulic characteristics of the sections agree. + +* **Correlation** Pearson's correlation between the source model and DEM cross-sections (closer to 1 is better). This metric captures how well the shape of the two sections match. + +* **Max Cross-Correlation** The maximum Pearson's correlation between the source model and DEM cross-sections across all horizontal shifts of the DEM section (closer to 1 is better). This metric captures how well the shape of the two sections match, however, it is insensitive to horizontal shifts in elevations. Compare to correlation in the horizontal shift example above. + +* **Spectral Correlation** Spectral correlation between source model and DEM cross-sections, as defined by the HydroErr library (https://github.com/BYU-Hydroinformatics/HydroErr/blob/42a84f3e006044f450edc7393ed54d59f27ef35b/HydroErr/HydroErr.py#L3615). Furthermore the metric has been transformed to range 0-1 and so that values closer to 1 are better. This metric captures how well the shape of the two sections match. + +* **Spectral Angle** Spectral angle between source model and DEM cross-sections, as defined by the HydroErr library (https://github.com/BYU-Hydroinformatics/HydroErr/blob/42a84f3e006044f450edc7393ed54d59f27ef35b/HydroErr/HydroErr.py#L3538). Furthermore the metric has been transformed to range 0-1 and so that values closer to 1 are better. This metric captures how well the shape of the two sections match. + +* **R-Squared** Coefficient of determination between the source model and DEM elevation series (closer to 1 is better). This metric captures how well the shape of the two sections match. + +* **Thalweg Elevation Difference** Source model invert minus the DEM invert/ Values closer to 0 are better, negative values reflect a higher DEM invert, and positive values reflect a higher source model invert. Since this metric is not scaled, acceptable ranges will vary from river to river. + + 4 - SRC development and FIM pre-processing ------------------------------------------ -(relevant endpoints: -:doc:`create_model_run_normal_depth `, -:doc:`run_incremental_normal_depth `, -:doc:`run_known_wse `, +(relevant endpoints: +:doc:`create_model_run_normal_depth `, +:doc:`run_incremental_normal_depth `, +:doc:`run_known_wse `, :doc:`create_fim_lib `) Once submodel geometry has been set up, you can run various discharges through the model and record the results. Ripple1d has several tools to develop -SRCs for a NWM reach. +SRCs for a NWM reach. * **Initial Normal Depth Run.** Discharges ranging from 1.2 times the reach high flow threshold to the reach 1% AEP discharge will be incrementally run @@ -187,4 +274,4 @@ SRCs for a NWM reach. Ripple1d generates HEC-RAS inundation depth grids for each of the known water surface elevation runs. These grids are cached along with their associated discharges and downstream conditions so that reach-scale FIM may be retrieved -as soon as a reach forecast is released. \ No newline at end of file +as soon as a reach forecast is released. diff --git a/ripple1d/consts.py b/ripple1d/consts.py index f232a458..735af8b3 100644 --- a/ripple1d/consts.py +++ b/ripple1d/consts.py @@ -51,3 +51,26 @@ SHOW_RAS = False HYDROFABRIC_CRS = 5070 + +TERRAIN_AGREEMENT_PRECISION = { + "inundation_overlap": 3, + "flow_area_overlap": 3, + "top_width_agreement": 3, + "flow_area_agreement": 3, + "hydraulic_radius_agreement": 3, + "mean": 2, + "std": 2, + "max": 2, + "min": 2, + "p_25": 2, + "p_50": 2, + "p_75": 2, + "rmse": 2, + "normalized_rmse": 3, + "r_squared": 3, + "spectral_angle": 3, + "spectral_correlation": 3, + "correlation": 3, + "max_cross_correlation": 3, + "thalweg_elevation_difference": 2, +} diff --git a/ripple1d/data_model.py b/ripple1d/data_model.py index 869c097b..180a0993 100644 --- a/ripple1d/data_model.py +++ b/ripple1d/data_model.py @@ -262,6 +262,15 @@ def ras_terrain_hdf(self): """RAS Terrain HDF file.""" return str(Path(self.terrain_directory) / f"{self.model_name}.hdf") + @property + def terrain_file(self) -> str | None: + """Terrain file path.""" + terrain_name = self.ripple1d_parameters.get("source_terrain") + if terrain_name is None: + return + suffix = terrain_name.split("/")[-1].replace(".vrt", ".tif") + return str(Path(self.terrain_directory) / f"{self.model_name}.{suffix}") + @property def fim_results_directory(self): """FIM results directory.""" @@ -355,6 +364,11 @@ def model_stac_json_file(self): """STAC JSON file.""" return self.derive_path(".model.stac.json") + @property + def terrain_agreement_file(self): + """Terrain agreement JSON file.""" + return self.derive_path(".terrain_agreement.json") + @dataclass class FlowChangeLocation: diff --git a/ripple1d/ops/ras_terrain.py b/ripple1d/ops/ras_terrain.py index 19451f12..f3abf97f 100644 --- a/ripple1d/ops/ras_terrain.py +++ b/ripple1d/ops/ras_terrain.py @@ -5,22 +5,29 @@ import json import logging import os +import re +from math import ceil, comb, pi from pathlib import Path import geopandas as gpd +import numpy as np import rasterio +import rioxarray +import xarray as xr from pyproj import CRS +from shapely import LineString from ripple1d.consts import ( MAP_DEM_BUFFER_DIST_FT, MAP_DEM_UNCLIPPED_SRC_URL, MAP_DEM_VERT_UNITS, METERS_PER_FOOT, + TERRAIN_AGREEMENT_PRECISION, ) -from ripple1d.data_model import NwmReachModel -from ripple1d.ras import create_terrain +from ripple1d.data_model import XS, NwmReachModel +from ripple1d.ras import RasGeomText, create_terrain from ripple1d.utils.dg_utils import clip_raster, reproject_raster -from ripple1d.utils.ripple_utils import fix_reversed_xs, xs_concave_hull +from ripple1d.utils.ripple_utils import fix_reversed_xs, resample_vertices, xs_concave_hull def get_geometry_mask(gdf_xs_conc_hull: str, MAP_DEM_UNCLIPPED_SRC_URL: str) -> gpd.GeoDataFrame: @@ -52,6 +59,7 @@ def create_ras_terrain( vertical_units: str = MAP_DEM_VERT_UNITS, resolution: float = None, resolution_units: str = None, + terrain_agreement_resolution: float = 3, ): """Create a RAS terrain file. @@ -119,9 +127,6 @@ def create_ras_terrain( if not os.path.exists(nwm_rm.terrain_directory): os.makedirs(nwm_rm.terrain_directory, exist_ok=True) - gdf_xs = gpd.read_file(nwm_rm.ras_gpkg_file, layer="XS").explode(ignore_index=True) - crs = gdf_xs.crs - mask = get_geometry_mask(nwm_rm.xs_concave_hull, terrain_source_url) # clip dem @@ -158,4 +163,372 @@ def create_ras_terrain( os.remove(src_dem_reprojected_localfile) nwm_rm.update_write_ripple1d_parameters({"source_terrain": terrain_source_url}) logging.info(f"create_ras_terrain complete") + + # Calculate terrain agreement metrics + agreement_path = compute_terrain_agreement_metrics(submodel_directory, terrain_agreement_resolution) + result["terrain_agreement"] = agreement_path return result + + +### Terrain Agreement Metrics ### + + +def compute_terrain_agreement_metrics(submodel_directory: str, max_sample_distance: float = 3): + """Compute a suite of agreement metrics between source model XS data and mapping DEM.""" + # Load model information + nwm_rm = NwmReachModel(submodel_directory) + dem_path = nwm_rm.terrain_file + + # Add DEM data to geom object + geom = RasGeomText.from_gpkg(nwm_rm.derive_path(".gpkg"), "", "") + section_data = sample_terrain(geom, dem_path, max_interval=max_sample_distance) + + # Compute agreement metrics + metrics = geom_agreement_metrics(section_data) + + # Save results and summary + with open(nwm_rm.terrain_agreement_file, "w") as f: + json.dump(metrics, f, indent=4) + nwm_rm.update_write_ripple1d_parameters({"terrain_agreement_summary": metrics["summary"]}) + return nwm_rm.terrain_agreement_file + + +def interpolater(coords: np.ndarray, stations: np.ndarray) -> np.ndarray: + """Interpolate faster than shapely linestring.""" + x = coords[:, 0] + y = coords[:, 1] + dists = np.cumsum(np.sqrt((np.diff(x, 1) ** 2) + (np.diff(y, 1) ** 2))) + dists = np.insert(dists, 0, 0) + newx = np.interp(stations, dists, x) + newy = np.interp(stations, dists, y) + return newx, newy + + +def sample_terrain(geom: RasGeomText, dem_path: str, max_interval: float = 3): + """Add DEM station_elevations to cross-sections.""" + section_data = {} + with rioxarray.open_rasterio(dem_path) as dem: + for section in geom.cross_sections: + station_elevation_points = np.array(geom.cross_sections[section].station_elevation_points) + stations = station_elevation_points[:, 0] + stations = resample_vertices(stations, max_interval) + rel_stations = stations - stations.min() + + xs, ys = interpolater(np.array(geom.cross_sections[section].coords), rel_stations) + tgt_x = xr.DataArray(xs, dims="points") + tgt_y = xr.DataArray(ys, dims="points") + el = dem.sel(band=1, x=tgt_x, y=tgt_y, method="nearest").values + dem_xs = np.column_stack((stations, el)) + src_resampled = np.interp(stations, station_elevation_points[:, 0], station_elevation_points[:, 1]) + src_xs = np.column_stack((stations, src_resampled)) + section_data[section] = {"dem_xs": dem_xs, "src_xs": src_xs} + return section_data + + +def geom_agreement_metrics(xs_data: dict) -> dict: + """Compute a suite of agreement metrics between source model XS data and a sampled DEM.""" + metrics = {"xs": {}, "summary": {}} + for section in xs_data: + metrics["xs"][section] = xs_agreement_metrics(xs_data[section]) + + # aggregate + metrics["summary"] = summarize_dict({i: metrics["xs"][i]["summary"] for i in metrics["xs"]}) # Summarize summaries + del metrics["summary"]["residuals"] # Averages are not applicable here + + return round_values(metrics) + + +def xs_agreement_metrics(xs: XS) -> dict: + """Compute a suite of agreement metrics between source model XS data and mapping DEM.""" + metrics = {} + + # Elevation-specific metrics and their summaries + metrics["elevation"] = variable_metrics(xs["src_xs"], xs["dem_xs"]) + metrics["summary"] = summarize_dict(metrics["elevation"]) + + # Whole XS metrics + src_xs_el = xs["src_xs"][:, 1] + dem_xs_el = xs["dem_xs"][:, 1] + metrics["summary"]["r_squared"] = r_squared(src_xs_el, dem_xs_el) + metrics["summary"]["spectral_angle"] = spectral_angle(src_xs_el, dem_xs_el) + metrics["summary"]["spectral_correlation"] = spectral_correlation(src_xs_el, dem_xs_el) + metrics["summary"]["correlation"] = correlation(src_xs_el, dem_xs_el) + metrics["summary"]["max_cross_correlation"] = cross_correlation(src_xs_el, dem_xs_el) + metrics["summary"]["thalweg_elevation_difference"] = thalweg_elevation_difference(src_xs_el, dem_xs_el) + + return metrics + + +def round_values(in_dict: dict) -> dict: + """Round values to keep precision consistent (recursive).""" + for k, v in in_dict.items(): + if isinstance(v, float): + in_dict[k] = round(v, TERRAIN_AGREEMENT_PRECISION[k]) + elif isinstance(v, dict): + in_dict[k] = round_values(v) + else: + in_dict[k] = None # at the moment, we have no other value types + return in_dict + + +def summarize_dict(in_dict: dict) -> dict: + """Aggregate metrics from a dictionary.""" + summary = {} + keys = list(in_dict.keys()) + submetrics = list(in_dict[keys[0]].keys()) + for m in submetrics: + if isinstance(in_dict[keys[0]][m], float): # average floats + summary[m] = sum([in_dict[k][m] for k in keys]) / len(keys) + elif isinstance(in_dict[keys[0]][m], dict): # Provide last entry for dicts + summary[m] = in_dict[keys[-1]][m] + return summary + + +def residual_metrics(residuals: np.ndarray) -> dict: + """Summary statistics on residuals between two arrays.""" + metrics = {} + metrics["mean"] = residuals.mean() + metrics["std"] = residuals.std() + metrics["max"] = residuals.max() + metrics["min"] = residuals.min() + metrics["p_25"] = np.percentile(residuals, 25) + metrics["p_50"] = np.percentile(residuals, 50) + metrics["p_75"] = np.percentile(residuals, 75) + metrics["rmse"] = np.sqrt((residuals * residuals).mean()) + metrics["normalized_rmse"] = metrics["rmse"] / (metrics["p_75"] - metrics["p_25"]) + return metrics + + +def ss(x: np.ndarray, y: np.ndarray) -> float: + """Calculate the sum of squares.""" + return np.sum((x - np.mean(x)) * (y - np.mean(y))) + + +def r_squared(a1: np.ndarray, a2: np.ndarray) -> float: + """Calculate R-squared for two series.""" + return ss(a1, a2) ** 2 / (ss(a1, a1) * ss(a2, a2)) + + +def spectral_angle(a1: np.ndarray, a2: np.ndarray, norm: bool = True) -> float: + """Calculate the spectral angle of two vectors. + + Following the approach of https://github.com/BYU-Hydroinformatics/HydroErr/blob/42a84f3e006044f450edc7393ed54d59f27ef35b/HydroErr/HydroErr.py#L3541 + """ + a = np.dot(a1, a2) + b = np.linalg.norm(a1) * np.linalg.norm(a2) + sa = np.arccos(a / b) + if norm: + sa = 1 - (abs(sa) / (pi / 2)) + return sa + + +def spectral_correlation(a1: np.ndarray, a2: np.ndarray, norm: bool = True) -> float: + """Calculate the spectral angle of two vectors. + + Following the approach of https://github.com/BYU-Hydroinformatics/HydroErr/blob/42a84f3e006044f450edc7393ed54d59f27ef35b/HydroErr/HydroErr.py#L3541 + """ + a = np.dot(a1 - np.mean(a1), a2 - np.mean(a2)) + b = np.linalg.norm(a1 - np.mean(a1)) + c = np.linalg.norm(a2 - np.mean(a2)) + e = b * c + if abs(a / e) > 1: + sc = 0 + else: + sc = np.arccos(a / e) + if norm: + sc = 1 - (abs(sc) / (pi / 2)) + return sc + + +def correlation(a1: np.ndarray, a2: np.ndarray) -> float: + """Calculate Pearson's correlation of two series.""" + num = np.sum((a1 - a1.mean()) * (a2 - a2.mean())) + denom = np.sqrt(np.sum(np.square(a1 - a1.mean())) * np.sum(np.square(a2 - a2.mean()))) + return num / denom + + +def cross_correlation(a1: np.ndarray, a2: np.ndarray) -> float: + """Calculate the maximum cross-correlation between two time series.""" + a1_detrend = a1 - a1.mean() + a2_detrend = a2 - a2.mean() + cross_corr = np.correlate(a1_detrend, a2_detrend, mode="full") + norm_factor = np.sqrt(np.sum(a1_detrend**2) * np.sum(a2_detrend**2)) + cross_corr /= norm_factor + return cross_corr.max() + + +def thalweg_elevation_difference(a1: np.ndarray, a2: np.ndarray) -> float: + """Calculate the elevation difference between the thalweg (low point) of two series.""" + return a1.min() - a2.min() + + +def variable_metrics(src_xs: np.ndarray, dem_xs: np.ndarray) -> dict: + """Calculate metrics for WSE values every half foot.""" + all_metrics = {} + residuals = src_xs - dem_xs # Calculate here to save compute + for wse in get_wses(src_xs): + tmp_src_xs = add_intersection_pts(src_xs, wse) + tmp_dem_xs = add_intersection_pts(dem_xs, wse) + metrics = {} + metrics["inundation_overlap"] = inundation_agreement(src_xs, dem_xs, wse) + metrics["flow_area_overlap"] = flow_area_overlap(src_xs, dem_xs, wse) + metrics["top_width_agreement"] = top_width_agreement(tmp_src_xs, tmp_dem_xs, wse) + metrics["flow_area_agreement"] = flow_area_agreement(tmp_src_xs, tmp_dem_xs, wse) + metrics["hydraulic_radius_agreement"] = hydraulic_radius_agreement(tmp_src_xs, tmp_dem_xs, wse) + resid_mask = (src_xs[:, 1] < wse) | (dem_xs[:, 1] < wse) + metrics["residuals"] = residual_metrics(residuals[resid_mask]) + all_metrics[wse] = metrics + return all_metrics + + +def get_wses(xs: np.ndarray, ramp_rate: int = 2, repeats: int = 5) -> list[float]: + """Derive grid of water surface elevations from minimum el to lowest cross-section endpoint.""" + start_el = xs[:, 1].min() + start_el = ceil(start_el * 2) / 2 # Round to nearest half foot + end_el = min((xs[0, 1], xs[-1, 1])) + end_el = ceil(end_el * 2) / 2 # Round to nearest half foot + + increments = np.arange(0, 10, 1) + increments = (ramp_rate**increments) * 0.5 + increments = np.repeat(increments, repeats) + increments = np.cumsum(increments) + + series = increments + start_el + if series[-1] < end_el: + np.append(series, end_el) + else: + series = series[series <= end_el] + return np.round(series, 1) + + +def add_intersection_pts(section_pts: np.ndarray, wse: float) -> np.ndarray: + """Add points to XS where WSE hits.""" + intersection = wse_intersection_pts(section_pts, wse) + if len(intersection) == 0: + return section_pts + inds = np.searchsorted(section_pts[:, 0], intersection[:, 0]) + return np.insert(section_pts, inds, intersection, axis=0) + + +def wse_intersection_pts(section_pts: np.ndarray, wse: float) -> list[tuple[float]]: + """Find where the cross-section terrain intersects the water-surface elevation.""" + intersection_pts = [] + + # Iterate through all pairs of points and find any points where the line would cross the wse + for i in range(len(section_pts) - 1): + p1 = section_pts[i] + p2 = section_pts[i + 1] + + if p1[1] > wse and p2[1] > wse: # No intesection + continue + elif p1[1] < wse and p2[1] < wse: # Both below wse + continue + elif p1[1] == wse or p2[1] == wse: # already vertex present + continue + + # Define line + m = (p2[1] - p1[1]) / (p2[0] - p1[0]) + b = p1[1] - (m * p1[0]) + + # Find intersection point with Cramer's rule + determinant = lambda a, b: (a[0] * b[1]) - (a[1] * b[0]) + div = determinant((1, 1), (-m, 0)) + tmp_y = determinant((b, wse), (-m, 0)) / div + tmp_x = determinant((1, 1), (b, wse)) / div + + intersection_pts.append((tmp_x, tmp_y)) + return np.array(intersection_pts) + + +def smape_series(a1: np.ndarray, a2: np.ndarray) -> float: + """Return the symmetric mean absolute percentage errror of two series.""" + num = np.abs(a1 - a2) + denom = np.abs(a1) + np.abs(a2) + return (np.sum(num / denom)) / len(a1) + + +def smape_single(a1: float, a2: float) -> float: + """Return the symmetric mean absolute percentage errror of two values.""" + if a1 == a2: + return 0.0 # handles zero denominator + return abs(a1 - a2) / (abs(a1) + abs(a2)) + + +def inundation_agreement(src_el: np.ndarray, dem_el: np.ndarray, wse: float) -> float: + """Calculate the percent of the cross-section with agreeing wet/dry.""" + dx = np.diff(src_el[:, 0], 1) + + src_wet = src_el[:, 1] < wse + src_wet = src_wet[1:] | src_wet[:-1] + dem_wet = dem_el[:, 1] < wse + dem_wet = dem_wet[1:] | dem_wet[:-1] + agree = src_wet & dem_wet + total = src_wet | dem_wet + return np.sum(agree * dx) / np.sum(total * dx) + + +def flow_area_overlap(src_el: np.ndarray, dem_el: np.ndarray, wse: float) -> float: + """Calculate the percent of unioned flow area that agree.""" + dx = np.diff(src_el[:, 0], 1) + + src_depths = np.clip(wse - src_el[:, 1], 0, None) + src_areas = ((src_depths[1:] + src_depths[:-1]) / 2) * dx + + dem_depths = np.clip(wse - dem_el[:, 1], 0, None) + dem_areas = ((dem_depths[1:] + dem_depths[:-1]) / 2) * dx + + combo = np.column_stack([src_areas, dem_areas]) + agree = np.min(combo, axis=1) + max_area = np.max(combo, axis=1) + return np.sum(agree) / np.sum(max_area) + + +def top_width_agreement(src_el: np.ndarray, dem_el: np.ndarray, wse: float) -> float: + """Calculate the agreement of the wetted top widths.""" + src_tw = get_wetted_top_width(src_el, wse) + dem_tw = get_wetted_top_width(dem_el, wse) + return 1 - smape_single(src_tw, dem_tw) + + +def flow_area_agreement(src_el: np.ndarray, dem_el: np.ndarray, wse: float) -> float: + """Calculate the agreement of the wetted top widths.""" + src_area = get_flow_area(src_el, wse) + dem_area = get_flow_area(dem_el, wse) + return 1 - smape_single(src_area, dem_area) + + +def hydraulic_radius_agreement(src_el: np.ndarray, dem_el: np.ndarray, wse: float) -> float: + """Calculate the agreement of the hydraulic radii.""" + src_radius = get_hydraulic_radius(src_el, wse) + dem_radius = get_hydraulic_radius(dem_el, wse) + return 1 - smape_single(src_radius, dem_radius) + + +def get_wetted_top_width(station_elevation_series: np.ndarray, wse: float) -> float: + """Derive wetted-top-width for a given stage.""" + dx = np.diff(station_elevation_series[:, 0], 1) + wet = station_elevation_series[:, 1] < wse + wet = wet[1:] | wet[:-1] + return np.sum(dx * wet) + + +def get_flow_area(station_elevation_series: np.ndarray, wse: float) -> float: + """Derive wetted-top-width for a given stage.""" + depths = np.clip(wse - station_elevation_series[:, 1], 0, None) + return np.trapezoid(depths, station_elevation_series[:, 0]) + + +def get_wetted_perimeter(station_elevation_series: np.ndarray, wse: float) -> float: + """Derive wetted-perimeter for a given stage.""" + station_elevation_series = station_elevation_series[station_elevation_series[:, 1] < wse] + diffs = np.diff(station_elevation_series, axis=0) + return np.sum(np.sqrt(np.square(diffs[:, 0]) + np.square(diffs[:, 1]))) + + +def get_hydraulic_radius(station_elevation_series: np.ndarray, wse: float) -> float: + """Derive wetted-perimeter for a given stage.""" + wp = get_wetted_perimeter(station_elevation_series, wse) + if wp == 0: + return 0 + a = get_flow_area(station_elevation_series, wse) + return a / wp diff --git a/ripple1d/ras.py b/ripple1d/ras.py index 017fa8ce..aa861b5c 100644 --- a/ripple1d/ras.py +++ b/ripple1d/ras.py @@ -51,6 +51,7 @@ assert_no_store_all_maps_error_message, decode, replace_line_in_contents, + resample_vertices, search_contents, text_block_from_start_end_str, ) diff --git a/ripple1d/utils/ripple_utils.py b/ripple1d/utils/ripple_utils.py index af8ef15a..d48d20a3 100644 --- a/ripple1d/utils/ripple_utils.py +++ b/ripple1d/utils/ripple_utils.py @@ -9,6 +9,7 @@ import boto3 import geopandas as gpd +import numpy as np import pandas as pd from dotenv import find_dotenv, load_dotenv from pyproj import CRS @@ -434,3 +435,20 @@ def assert_no_store_all_maps_error_message(compute_message_file: str): raise RASStoreAllMapsError( f"{repr(line)} found in {compute_message_file}. Full file content:\n{content}\n^^^ERROR^^^" ) + + +def resample_vertices(stations: np.ndarray, max_interval: float) -> np.ndarray: + """Resample a set of stations so that no gaps are larger than max_interval.""" + i = 0 + max_iter = 1e10 + while i < (len(stations) - 1) and i < max_iter: + gap = stations[i + 1] - stations[i] + if gap <= max_interval: + i += 1 + continue + subdivisions = (gap // max_interval) - 1 + modulo = gap - (subdivisions * max_interval) + new_pts = np.arange(stations[i] + (modulo / 2), stations[i + 1], max_interval) + stations = np.insert(stations, i + 1, new_pts) + i += len(new_pts) + 1 + return stations