Maurice_paper_2_main

examples/build.yml

@@ -53,20 +53,20 @@ setup_irrigation_sources:
 
 setup_create_farms_simple:
 
+setup_household_characteristics:
+
 setup_crops_from_source:
   source: MIRCA2000
+
+determine_crop_area_fractions:
+
+setup_farmer_crop_calendar:
+  year: 2000
+  reduce_crops: true
+  replace_base: false
 
 setup_farmer_characteristics_simple:
-  risk_aversion_mean: 2
-  risk_aversion_standard_deviation: 1
-  discount_rate: 0.05
   interest_rate: 0.05
-  irrigation_choice:
-    canal: 80
-    well: 20
-
-setup_cultivation_costs:
-  cultivation_costs: 0
 
 setup_crop_prices:
   crop_prices: FAO_stat

geb/HRUs.py

@@ -13,13 +13,33 @@
 import numpy as np
 import zarr.convenience
 from geb.workflows import AsyncForcingReader
+from scipy.spatial import cKDTree
 
 try:
     import cupy as cp
 except (ModuleNotFoundError, ImportError):
     pass
 
 
+def determine_nearest_river_cell(river_grid, HRU_to_grid):
+    threshold = 15
+    valid_mask = river_grid != -9999
+
+    valid_indices = np.argwhere(valid_mask)
+    valid_values = river_grid[valid_mask]
+
+    above_threshold_mask = valid_values > threshold
+    above_threshold_indices = valid_indices[above_threshold_mask]
+    above_threshold_indices_in_valid = np.flatnonzero(above_threshold_mask)
+
+    tree = cKDTree(above_threshold_indices)
+    distances, indices_in_above = tree.query(valid_indices)
+
+    nearest_indices_in_valid = above_threshold_indices_in_valid[indices_in_above]
+
+    return nearest_indices_in_valid[HRU_to_grid]
+
+
 def load_grid(filepath, layer=1, return_transform_and_crs=False):
     if filepath.suffix == ".tif":
         warnings.warn("tif files are now deprecated. Consider rebuilding the model.")
@@ -568,6 +588,11 @@ def __init__(self, data, model) -> None:
                     self.data.get_save_state_path(), "HRU.unmerged_HRU_indices.npz"
                 )
             )["data"]
+            self.nearest_river_grid_cell = np.load(
+                os.path.join(
+                    self.data.get_save_state_path(), "HRU.nearest_river_grid_cell.npz"
+                )
+            )["data"]
         else:
             (
                 self.land_use_type,
@@ -577,12 +602,22 @@ def __init__(self, data, model) -> None:
                 self.grid_to_HRU,
                 self.unmerged_HRU_indices,
             ) = self.create_HRUs()
+
+            river_grid = load_grid(
+                self.model.files["grid"]["routing/kinematic/upstream_area"]
+            )
+
+            self.nearest_river_grid_cell = determine_nearest_river_cell(
+                river_grid, self.HRU_to_grid
+            )
+
             self.register_initial_data("HRU.land_use_type")
             self.register_initial_data("HRU.land_use_ratio")
             self.register_initial_data("HRU.land_owners")
             self.register_initial_data("HRU.HRU_to_grid")
             self.register_initial_data("HRU.grid_to_HRU")
             self.register_initial_data("HRU.unmerged_HRU_indices")
+            self.register_initial_data("HRU.nearest_river_grid_cell")
         if self.model.use_gpu:
             self.land_use_type = cp.array(self.land_use_type)
         BaseVariables.__init__(self)