Merge pull request #276 from Dewberry/feature/terrain_metrics

Feature/terrain metrics
Dewberry · Dec 20, 2024 · 4f05684 · 4f05684
2 parents 71c5122 + d19b5d1
commit 4f05684
Show file tree

Hide file tree

Showing 22 changed files with 533 additions and 17 deletions.
diff --git a/docs/images/terrain_agreement/1.png b/docs/images/terrain_agreement/1.png
diff --git a/docs/images/terrain_agreement/2.png b/docs/images/terrain_agreement/2.png
diff --git a/docs/images/terrain_agreement/3.png b/docs/images/terrain_agreement/3.png
diff --git a/docs/images/terrain_agreement/4.png b/docs/images/terrain_agreement/4.png
diff --git a/docs/images/terrain_agreement/5.png b/docs/images/terrain_agreement/5.png
diff --git a/docs/images/terrain_agreement/6.png b/docs/images/terrain_agreement/6.png
diff --git a/docs/images/terrain_agreement/7.png b/docs/images/terrain_agreement/7.png
diff --git a/docs/images/terrain_agreement/8.png b/docs/images/terrain_agreement/8.png
diff --git a/docs/source/images/terrain_agreement/1.png b/docs/source/images/terrain_agreement/1.png
diff --git a/docs/source/images/terrain_agreement/2.png b/docs/source/images/terrain_agreement/2.png
diff --git a/docs/source/images/terrain_agreement/3.png b/docs/source/images/terrain_agreement/3.png
diff --git a/docs/source/images/terrain_agreement/4.png b/docs/source/images/terrain_agreement/4.png
diff --git a/docs/source/images/terrain_agreement/5.png b/docs/source/images/terrain_agreement/5.png
diff --git a/docs/source/images/terrain_agreement/6.png b/docs/source/images/terrain_agreement/6.png
diff --git a/docs/source/images/terrain_agreement/7.png b/docs/source/images/terrain_agreement/7.png
diff --git a/docs/source/images/terrain_agreement/8.png b/docs/source/images/terrain_agreement/8.png
diff --git 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 <endpoints/conflate_model>`, 
+(relevant endpoints: :doc:`conflate_model <endpoints/conflate_model>`,
 :doc:`compute_conflation_metrics <endpoints/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 <endpoints/extract_submodel>`, 
+(relevant endpoints: :doc:`extract_submodel <endpoints/extract_submodel>`,
 :doc:`create_ras_terrain <endpoints/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
 <https://data.usgs.gov/datacatalog/data/USGS:3a81321b-c153-416f-98b7-cc8e5f0e17c3>`_
 
+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 <endpoints/create_model_run_normal_depth>`, 
-:doc:`run_incremental_normal_depth <endpoints/run_incremental_normal_depth>`, 
-:doc:`run_known_wse <endpoints/run_known_wse>`, 
+(relevant endpoints:
+:doc:`create_model_run_normal_depth <endpoints/create_model_run_normal_depth>`,
+:doc:`run_incremental_normal_depth <endpoints/run_incremental_normal_depth>`,
+:doc:`run_known_wse <endpoints/run_known_wse>`,
 :doc:`create_fim_lib <endpoints/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.
+as soon as a reach forecast is released.
diff --git 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
@@ -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: