dynamic land use changes for damage calculation

geb/HRUs.py

@@ -216,31 +216,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"]
@@ -503,14 +503,10 @@ 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
@@ -796,6 +792,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/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
@@ -330,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/model.py

@@ -84,12 +84,13 @@ def __init__(
         config: dict,
         files: dict,
         spinup: bool = False,
+        crs=4326,
         use_gpu: bool = False,
         gpu_device=0,
         timing=False,
-        crs=4326,
         mode="w",
     ):
+        self.crs = crs
         self.timing = timing
         assert mode in ("w", "r")
         self.mode = mode