added v1 exception to streamflow, and statistics for summed Q_prime

marquette/__main__.py

@@ -36,13 +36,13 @@ def main(cfg: DictConfig) -> None:
 
         log.info(f"Mapping {cfg.zone} Streamflow to TMs")
         create_TMs(cfg, edges)
+
+        log.info("Mapping Lake Pour Points to Edges")
+        map_lake_points(cfg, edges)
 
         log.info("Converting Streamflow to zarr")
         write_streamflow(cfg, edges)
 
-        log.info("Mapping Lake Pour Points to Edges")
-        map_lake_points(cfg, edges)
-
         log.info("Running Data Post-Processing Extensions")
         run_extensions(cfg, edges)
 
@@ -103,6 +103,15 @@ def run_extensions(cfg: DictConfig, edges: zarr.Group) -> None:
             log.info("q_prime_sum already exists in zarr format")
         else:
             calculate_q_prime_summation(cfg, edges)
+
+    if "q_prime_sum_stats" in cfg.extensions:
+        from marquette.merit.extensions import calculate_q_prime_sum_stats
+
+        log.info("Adding q_prime_sum statistics to your MERIT River Graph")
+        if "summed_q_prime_median" in edges:
+            log.info("q_prime_sum statistics already exists in zarr format")
+        else:
+            calculate_q_prime_sum_stats(cfg, edges)
 
 
 if __name__ == "__main__":

marquette/conf/config.yaml

@@ -1,6 +1,6 @@
 name: MERIT
 data_path: /projects/mhpi/data/${name}
-zone: 73
+zone: 71
 create_edges:
   buffer: 0.3334
   dx: 2000
@@ -22,10 +22,10 @@ create_TMs:
     TM: ${data_path}/zarr/TMs/sparse_MERIT_FLOWLINES_${zone}
     shp_files: ${data_path}/raw/basins/cat_pfaf_${zone}_MERIT_Hydro_v07_Basins_v01_bugfix1.shp
 create_streamflow:
-  version: merit_conus_v6.14
+  version: merit_conus_v1.0
   data_store: ${data_path}/streamflow/zarr/${create_streamflow.version}/${zone}
   obs_attributes: ${data_path}/gage_information/MERIT_basin_area_info
-  predictions: /projects/mhpi/yxs275/DM_output/dPL_local_daymet_new_attr_water_loss_v6v10_random_batch_all_merit_forward/
+  predictions: /projects/mhpi/yxs275/DM_output/water_loss_model/dPL_local_daymet_new_attr_RMSEloss_with_log_2800
   start_date: 01-01-1980
   end_date: 12-31-2020
 map_lake_points:
@@ -37,6 +37,7 @@ extensions:
   - global_dhbv_static_inputs
   - incremental_drainage_area
   - q_prime_sum
+  - q_prime_sum_stats
 # Hydra Config ------------------------------------------------------------------------#
 hydra:
   help:

marquette/merit/_map_lake_points.py

@@ -33,6 +33,10 @@ def _map_lake_points(cfg: DictConfig, edges: zarr.Group) -> None:
     lake_area = np.full_like(edges_comids, -1, dtype=np.float32)
     vol_total = np.full_like(edges_comids, -1, dtype=np.float32)
     depth_avg = np.full_like(edges_comids, -1, dtype=np.float32)
+    vol_res = np.full(len(edges_comids), -1, dtype=np.int32)
+    elevation = np.full_like(edges_comids, -1, dtype=np.int32)
+    shore_dev = np.full_like(edges_comids, -1, dtype=np.float32)
+    dis_avg = np.full_like(edges_comids, -1, dtype=np.float32)
 
     for idx, lake_id in enumerate(tqdm(
         lake_comids,
@@ -52,11 +56,20 @@ def _map_lake_points(cfg: DictConfig, edges: zarr.Group) -> None:
             lake_area[edge_idx] = lake_row["Lake_area"]
             vol_total[edge_idx] = lake_row["Vol_total"]
             depth_avg[edge_idx] = lake_row["Depth_avg"]
+
+            vol_res[edge_idx] = lake_row["Vol_res"]
+            elevation[edge_idx] = lake_row["Elevation"]
+            shore_dev[edge_idx] = lake_row["Shore_dev"]
+            dis_avg[edge_idx] = lake_row["Dis_avg"]
 
     edges.array("hylak_id", data=hylak_id)
     edges.array("grand_id", data=grand_id)
     edges.array("lake_area", data=lake_area)
     edges.array("vol_total", data=vol_total)
     edges.array("depth_avg", data=depth_avg)
+    edges.array("vol_res", data=vol_res)
+    edges.array("elevation", data=elevation)
+    edges.array("shore_dev", data=shore_dev)
+    edges.array("dis_avg", data=dis_avg)
 
     log.info("Wrote Lake data for zones to zarr")

marquette/merit/_streamflow_conversion_functions.py

@@ -130,7 +130,7 @@ def calculate_merit_flow(cfg: DictConfig, edges: zarr.hierarchy.Group) -> None:
     version = int(
         cfg.create_streamflow.version.lower().split("_v")[1][0]
     )  # getting the version number
-    if version >= 3:
+    if version >= 3 or version == 1:
         log.info(msg="Reading Zarr Store")
         zone_keys = [
             key for key in streamflow_predictions_root.keys() if str(cfg.zone) in key
@@ -144,7 +144,6 @@ def calculate_merit_flow(cfg: DictConfig, edges: zarr.hierarchy.Group) -> None:
         file_runoff = np.transpose(np.concatenate(zone_runoff))
         del zone_comids
         del zone_runoff
-
     else:
         log.info("Reading Zarr Store")
         file_runoff = np.transpose(streamflow_predictions_root.Runoff)

marquette/merit/extensions.py

@@ -318,3 +318,35 @@ def calculate_q_prime_summation(cfg: DictConfig, edges: zarr.Group) -> None:
         name="summed_q_prime",
         data=q_prime_np.transpose(1, 0),
     )
+
+
+def calculate_q_prime_sum_stats(cfg: DictConfig, edges: zarr.Group) -> None:
+    """Creates Q` summed data for all edges in a given MERIT zone
+
+    Parameters:
+    ----------
+    cfg: DictConfig
+        The configuration object.
+    edges: zarr.Group
+        The edges group in the MERIT zone
+    """
+    try:
+        summed_q_prime: np.ndarray = edges.summed_q_prime[:]
+    except AttributeError:
+        raise AttributeError("summed_q_prime data not found")
+    edges.array(
+        name="summed_q_prime_median",
+        data=np.median(summed_q_prime, axis=0),
+    )    
+    edges.array(
+        name="summed_q_prime_std",
+        data=np.std(summed_q_prime, axis=0),
+    )  
+    edges.array(
+        name="summed_q_prime_p90",
+        data=np.percentile(summed_q_prime, 90, axis=0),
+    )  
+    edges.array(
+        name="summed_q_prime_p10",
+        data=np.percentile(summed_q_prime, 10, axis=0),
+    )