Merge remote-tracking branch 'origin/main' into maurice_paper_2_main

.github/workflows/docs_test.yml

@@ -12,7 +12,7 @@ env:
 jobs:
   build-docs:
     name: Build Sphinx Documentation
-    runs-on: ubuntu-24.04
+    runs-on: ubuntu-latest
     steps:
       # Check out the repository
       - uses: actions/checkout@v4

.github/workflows/publish.yml

@@ -7,7 +7,7 @@ on:
 jobs:
   pypi-publish:
     name: Publish release to PyPI
-    runs-on: ubuntu-24.04
+    runs-on: ubuntu-latest
     environment:
       name: "release"
       url: https://pypi.org/project/geb/

.github/workflows/pytest.yml

@@ -8,19 +8,17 @@ on:
 
 jobs:
   test:
-    runs-on: ubuntu-24.04
+    runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v4
     - name: "Set up Python"
       uses: actions/setup-python@v5
       with:
         python-version-file: "pyproject.toml"
-    - name: Update TBB
-      run: sudo apt-get update && sudo apt-get install -y libtbb-dev
     - name: Set up uv
       # Install latest uv version using the installer
       run: curl -LsSf https://astral.sh/uv/install.sh | sh
     - name: Install the project
       run: uv sync --extra dev
     - name: Run tests
-      run: uv run pytest
+      run: uv run pytest -s

.github/workflows/ruff-check.yml

@@ -2,7 +2,7 @@ name: Ruff check
 on: [push, pull_request]
 jobs:
   ruff:
-    runs-on: ubuntu-24.04
+    runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
       - uses: chartboost/ruff-action@v1

.github/workflows/ruff-format.yml

@@ -2,7 +2,7 @@ name: Ruff format
 on: [push, pull_request]
 jobs:
   ruff:
-    runs-on: ubuntu-24.04
+    runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
       - uses: chartboost/ruff-action@v1

README.md

@@ -1,13 +1,13 @@
-Hello! Welcome to GEB! You can find full documentation [here](https://geb-model.github.io/GEB/).
-
 ## Installation
 
-The model can be installed with pip as follows:
+GEB can be installed with pip, including all dependencies on Windows, Linux and Mac OS X.
 
-`pip install geb`
+```bash
+pip install geb
+```
 
 ## Overview
-GEB (Geographic Environmental and Behavioural model) aims to simulate both the environment, the individual behaviour of people and their interactions at small and large scale. The model does so through a "deep" coupling of an agent-based model which simulates millions of individual people or households, a hydrological model, a vegetation model and a hydrodynamic model.
+GEB (Geographical Environmental and Behavioural model) simulates the environment (e.g., hydrology, floods), the individual people, households and orginizations as well as their interactions at both small and large scale. The model does so through a "deep" coupling of an agent-based model a hydrological model, a vegetation model and a hydrodynamic model. You can find full documentation [here](https://geb-model.github.io/GEB/).
 
 The figure below shows a schematic overview of the model agent-based and hydrological model.
 
@@ -17,11 +17,11 @@ The figure below shows a schematic overview of the model agent-based and hydrolo
 
 ### Model framework
 
-de Bruijn, J. A., Smilovic, M., Burek, P., Guillaumot, L., Wada, Y., and Aerts, J. C. J. H.: GEB v0.1: a large-scale agent-based socio-hydrological model – simulating 10 million individual farming households in a fully distributed hydrological model, Geosci. Model Dev., 16, 2437–2454, [https://doi.org/10.5194/gmd-16-2437-2023](https://doi.org/10.5194/gmd-16-2437-2023), 2023.
+> de Bruijn, J. A., Smilovic, M., Burek, P., Guillaumot, L., Wada, Y., and Aerts, J. C. J. H.: GEB v0.1: a large-scale agent-based socio-hydrological model – simulating 10 million individual farming households in a fully distributed hydrological model, Geosci. Model Dev., 16, 2437–2454, [https://doi.org/10.5194/gmd-16-2437-2023](https://doi.org/10.5194/gmd-16-2437-2023), 2023.
 
 ### Applications
 
-Kalthof, M. W. M. L., de Bruijn, J., de Moel, H., Kreibich, H., and Aerts, J. C. J. H.: Adaptive Behavior of Over a Million Individual Farmers Under Consecutive Droughts: A Large-Scale Agent-Based Modeling Analysis in the Bhima Basin, India, EGUsphere preprint, [https://doi.org/10.5194/egusphere-2024-1588](https://doi.org/10.5194/egusphere-2024-1588), 2024.
+> Kalthof, M. W. M. L., de Bruijn, J., de Moel, H., Kreibich, H., and Aerts, J. C. J. H.: Adaptive Behavior of Over a Million Individual Farmers Under Consecutive Droughts: A Large-Scale Agent-Based Modeling Analysis in the Bhima Basin, India, EGUsphere preprint, [https://doi.org/10.5194/egusphere-2024-1588](https://doi.org/10.5194/egusphere-2024-1588), 2024.
 
 ## Building on the shoulders of giants
 

examples/model.yml

@@ -1,6 +1,6 @@
 # The {GEB_PACKAGE_DIR} prefix is a special prefix that is replaced by the path to the geb package.
 # You can also find the reasonable_defaults/model.yml file relative to this file.
-inherits: "{GEB_PACKAGE_DIR}/examples/reasonable_defaults/model.yml"
+inherits: "{GEB_PACKAGE_DIR}/reasonable_default_config.yml"
 
 ###
 # This section contains several general model settings.

geb/HRUs.py

@@ -236,31 +236,31 @@ def __init__(self, data, model):
         self.data = data
         self.model = model
         self.scaling = 1
-        mask, transform, self.crs = load_grid(
+        mask, self.transform, self.crs = load_grid(
             self.model.files["grid"]["areamaps/grid_mask"],
             return_transform_and_crs=True,
         )
         self.mask = mask.astype(bool)
-        self.gt = transform.to_gdal()
+        self.gt = self.transform.to_gdal()
         self.bounds = (
-            transform.c,
-            transform.f + transform.e * mask.shape[0],
-            transform.c + transform.a * mask.shape[1],
-            transform.f,
+            self.transform.c,
+            self.transform.f + self.transform.e * mask.shape[0],
+            self.transform.c + self.transform.a * mask.shape[1],
+            self.transform.f,
         )
         self.lon = np.linspace(
-            transform.c + transform.a / 2,
-            transform.c + transform.a * mask.shape[1] - transform.a / 2,
+            self.transform.c + self.transform.a / 2,
+            self.transform.c + self.transform.a * mask.shape[1] - self.transform.a / 2,
             mask.shape[1],
         )
         self.lat = np.linspace(
-            transform.f + transform.e / 2,
-            transform.f + transform.e * mask.shape[0] - transform.e / 2,
+            self.transform.f + self.transform.e / 2,
+            self.transform.f + self.transform.e * mask.shape[0] - self.transform.e / 2,
             mask.shape[0],
         )
 
-        assert math.isclose(transform.a, -transform.e)
-        self.cell_size = transform.a
+        assert math.isclose(self.transform.a, -self.transform.e)
+        self.cell_size = self.transform.a
 
         self.cell_area_uncompressed = load_grid(
             self.model.files["grid"]["areamaps/cell_area"]
@@ -523,19 +523,28 @@ def __init__(self, data, model) -> None:
 
         self.scaling = submask_height // self.data.grid.shape[0]
         assert submask_width // self.data.grid.shape[1] == self.scaling
-        self.gt = (
-            self.data.grid.gt[0],
-            self.data.grid.gt[1] / self.scaling,
-            self.data.grid.gt[2],
-            self.data.grid.gt[3],
-            self.data.grid.gt[4],
-            self.data.grid.gt[5] / self.scaling,
-        )
+
+        self.transform = self.data.grid.transform * Affine.scale(1 / self.scaling)
+
+        self.gt = self.transform.to_gdal()
 
         self.mask = self.data.grid.mask.repeat(self.scaling, axis=0).repeat(
             self.scaling, axis=1
         )
         self.cell_size = self.data.grid.cell_size / self.scaling
+
+        # get lats and lons for subgrid
+        self.lon = np.linspace(
+            self.gt[0] + self.cell_size / 2,
+            self.gt[0] + self.cell_size * submask_width - self.cell_size / 2,
+            submask_width,
+        )
+        self.lat = np.linspace(
+            self.gt[3] + self.cell_size / 2,
+            self.gt[3] + self.cell_size * submask_height - self.cell_size / 2,
+            submask_height,
+        )
+
         if self.model.load_initial_data:
             self.land_use_type = np.load(
                 os.path.join(self.data.get_save_state_path(), "HRU.land_use_type.npz")
@@ -818,6 +827,8 @@ def decompress(self, HRU_array: np.ndarray) -> np.ndarray:
             HRU_array = HRU_array.get()
         if np.issubdtype(HRU_array.dtype, np.integer):
             nanvalue = -1
+        elif np.issubdtype(HRU_array.dtype, bool):
+            nanvalue = False
         else:
             nanvalue = np.nan
         outarray = HRU_array[self.unmerged_HRU_indices]

geb/__init__.py

@@ -1,3 +1,27 @@
 """GEB simulates the environment, the individual behaviour of people, households and organizations - including their interactions - at small and large scale."""
 
-__version__ = "1.0.0b2"
+__version__ = "1.0.0b4"
+
+import os
+import platform
+import sys
+from pathlib import Path
+from numba import config
+import faulthandler
+
+
+if platform.system() != "Windows":
+    # Modify LD_LIBRARY_PATH on Unix-like systems (Linux, macOS)
+    import tbb  # noqa: F401
+
+    tbb_path = Path(sys.prefix) / "lib" / "libtbb.so"
+    assert tbb_path.exists(), f"tbb shared library not found at {tbb_path}"
+    os.environ["LD_LIBRARY_PATH"] = str(tbb_path)
+
+# set threading layer to tbb, this is much faster than other threading layers
+config.THREADING_LAYER = "tbb"
+
+# set environment variable for GEB package directory
+os.environ["GEB_PACKAGE_DIR"] = str(Path(__file__).parent)
+
+faulthandler.enable()

geb/agents/households.py

@@ -2,15 +2,16 @@
 import geopandas as gpd
 import calendar
 from .general import AgentArray, downscale_volume, AgentBaseClass
-from ..hydrology.landcover import SEALED
+from ..hydrology.landcover import (
+    SEALED,
+    FOREST,
+)
 import pandas as pd
 from os.path import join
 from damagescanner.core import object_scanner
 import json
-import xarray as xr
 import rioxarray
 from rasterio.features import shapes
-import rasterio
 from shapely.geometry import shape
 
 try:
@@ -19,7 +20,7 @@
     pass
 
 
-def from_landuse_raster_to_polygon(rasterdata, landuse_category):
+def from_landuse_raster_to_polygon(mask, transform, crs):
     """
     Convert raster data into separate GeoDataFrames for specified land use values.
 
@@ -30,32 +31,14 @@ def from_landuse_raster_to_polygon(rasterdata, landuse_category):
     Returns:
     - Geodataframe
     """
-    data = rasterdata["data"].values
-    data = data.astype(np.uint8)
 
-    y_coords = rasterdata.coords["y"].values
-    x_coords = rasterdata.coords["x"].values
-
-    transform = rasterio.transform.from_origin(
-        x_coords[0],
-        y_coords[0],
-        abs(x_coords[1] - x_coords[0]),
-        abs(y_coords[1] - y_coords[0]),
-    )
-
-    mask = data == landuse_category
-
-    shapes_gen = shapes(data, mask=mask, transform=transform)
+    shapes_gen = shapes(mask.astype(np.uint8), mask=mask, transform=transform)
 
     polygons = []
-    for geom, value in shapes_gen:
-        if value == landuse_category:
-            polygons.append(shape(geom))
+    for geom, _ in shapes_gen:
+        polygons.append(shape(geom))
 
-    gdf = gpd.GeoDataFrame(
-        {"value": [landuse_category] * len(polygons), "geometry": polygons}
-    )
-    gdf.set_crs(epsg=4326, inplace=True)
+    gdf = gpd.GeoDataFrame({"geometry": polygons}, crs=crs)
 
     return gdf
 
@@ -81,15 +64,18 @@ def __init__(self, model, agents, reduncancy: float) -> None:
         self.rail["object_type"] = "rail"
 
         # Load landuse and make turn into polygons
-        self.landuse = xr.open_zarr(
-            self.model.files["region_subgrid"][
-                "landsurface/full_region_cultivated_land"
-            ]
+        self.forest = from_landuse_raster_to_polygon(
+            self.model.data.HRU.decompress(self.model.data.HRU.land_use_type == FOREST),
+            self.model.data.HRU.transform,
+            self.model.crs,
         )
-        self.forest = from_landuse_raster_to_polygon(self.landuse, 0)
         self.forest["object_type"] = "forest"
 
-        self.agriculture = from_landuse_raster_to_polygon(self.landuse, 1)
+        self.agriculture = from_landuse_raster_to_polygon(
+            self.model.data.HRU.decompress(self.model.data.HRU.land_owners != -1),
+            self.model.data.HRU.transform,
+            self.model.crs,
+        )
         self.agriculture["object_type"] = "agriculture"
 
         # Load maximum damages
@@ -270,8 +256,6 @@ def initiate(self) -> None:
             n=self.n, max_n=self.max_n, fill_value=1, dtype=np.float32
         )
 
-        self.buildings = gpd.read_file(self.model.files["geoms"]["assets/buildings"])
-
     def flood(self, flood_map, simulation_root, return_period=None):
         if return_period is not None:
             flood_path = join(simulation_root, f"hmax RP {int(return_period)}.tif")
@@ -332,7 +316,7 @@ def flood(self, flood_map, simulation_root, return_period=None):
         print(f"damages to rail are: {total_damages_rail}")
 
         total_flood_damages = (
-            +total_damage_structure
+            total_damage_structure
             + total_damages_content
             + total_damages_roads
             + total_damages_rail

geb/agents/market.py

@@ -170,8 +170,18 @@ def crop_prices(self) -> np.ndarray:
         ):
             return self.get_modelled_crop_prices()
         else:
-            index = self._crop_prices[0].get(self.model.current_time)
-            return self._crop_prices[1][index]
+            if self._crop_prices[0] is None:
+                print("WARNING: Using static crop prices")
+                return np.full(
+                    (
+                        len(self.cumulative_inflation_per_region),
+                        len(self.agents.crop_farmers.crop_ids),
+                    ),
+                    self._crop_prices[1],
+                )
+            else:
+                index = self._crop_prices[0].get(self.model.current_time)
+                return self._crop_prices[1][index]
 
     @property
     def year_index(self) -> int:

geb/agents/reservoir_operators.py

@@ -152,3 +152,7 @@ def get_available_water_reservoir_command_areas(self, reservoir_storage_m3):
 
     def step(self) -> None:
         return None
+
+    @property
+    def storage(self):
+        return self.model.data.grid.storage