Merge pull request #56 from nismod/feature/damages

Interpolate damage curves, fix polygon intersection extent

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name="nismod-snail"
-version="0.4.0"
+version="0.4.1"
 license={file = "LICENSE"}
 description="The spatial networks impact assessment library"
 readme="README.md"

src/snail/damages.py

@@ -94,7 +94,7 @@ def from_csv(
                 damage_col: "damage",
             }
         )
-        return PiecewiseLinearDamageCurve(curve_data)
+        return cls(curve_data)
 
     @classmethod
     def from_excel(
@@ -142,7 +142,47 @@ def from_excel(
                 damage_col: "damage",
             }
         )
-        return PiecewiseLinearDamageCurve(curve_data)
+        return cls(curve_data)
+
+    @classmethod
+    def interpolate(
+        cls,
+        a,
+        b,
+        factor: float,
+    ):
+        """Interpolate damage values between two curves
+
+        ```
+        new_curve_damage = a_damage + ((b_damage - a_damage) * factor)
+        ```
+
+        Parameters
+        ----------
+        a: PiecewiseLinearDamageCurve
+        b: PiecewiseLinearDamageCurve
+        factor: float
+            Interpolation factor, used to calculate the new curve
+
+        Returns
+        -------
+        PiecewiseLinearDamageCurve
+        """
+        # find the sorted, unique values of intensity used by both curves
+        intensity = numpy.unique(numpy.concatenate((a.intensity, b.intensity)))
+        # calculate the damage fraction at each intensity
+        a_damage = a.damage_fraction(intensity)
+        b_damage = b.damage_fraction(intensity)
+        # interpolate at each point
+        damage = a_damage + ((b_damage - a_damage) * factor)
+        return cls(
+            pandas.DataFrame(
+                {
+                    "intensity": intensity,
+                    "damage": damage,
+                }
+            )
+        )
 
     def damage_fraction(self, exposure: numpy.array) -> numpy.array:
         """Evaluate damage fraction for exposure to a given hazard intensity"""

src/snail/intersection.py

@@ -208,13 +208,13 @@ def split_polygons(
 def generate_grid_boxes(grid):
     a, b, c, d, e, f = grid.transform
     idx = numpy.arange(grid.width * grid.height)
-    i, j = numpy.unravel_index(idx, (grid.height, grid.width))
-    ulx = i * a + j * b + c
-    uly = i * d + j * e + f
-    lrx = (i + 1) * a + (j + 1) * b + c
-    lry = (i + 1) * d + (j + 1) * e + f
+    i, j = numpy.unravel_index(idx, (grid.width, grid.height))
+    xmin = i * a + j * b + c
+    ymax = i * d + j * e + f
+    xmax = (i + 1) * a + (j + 1) * b + c
+    ymin = (i + 1) * d + (j + 1) * e + f
     return geopandas.GeoDataFrame(
-        data={}, geometry=box(ulx, lry, lrx, uly), crs=grid.crs
+        data={}, geometry=box(xmin, ymin, xmax, ymax), crs=grid.crs
     )
 
 
@@ -327,6 +327,14 @@ def apply_indices(
     index_i="index_i",
     index_j="index_j",
 ) -> geopandas.GeoDataFrame:
+    if features.empty:
+        logging.info("Returning empty dataframe")
+        # return an empty dataframe with the expected columns
+        empty_df = features.copy()
+        empty_df[index_i] = numpy.array([], dtype="int64")
+        empty_df[index_j] = numpy.array([], dtype="int64")
+        return empty_df
+
     def f(geom, *args, **kwargs):
         return get_indices(geom, grid, index_i, index_j)
 

tests/test_damages.py

@@ -192,3 +192,40 @@ def test_linear_curve_scale_x_down(curve):
 
     expected = numpy.array([0, 1, 2, 3])
     assert_allclose(decreased.intensity, expected)
+
+
+def test_interpolation(curve):
+    curve_050 = curve.scale_y(0.5)
+    curve_100 = curve
+    expected = curve.scale_y(0.75)
+    actual = PiecewiseLinearDamageCurve.interpolate(curve_050, curve_100, 0.5)
+    assert actual == expected
+
+
+def test_interpolation_unaligned():
+    a = PiecewiseLinearDamageCurve(
+        pandas.DataFrame(
+            {
+                "intensity": [0.0, 10.0],
+                "damage": [0.0, 1.0],
+            }
+        )
+    )
+    b = PiecewiseLinearDamageCurve(
+        pandas.DataFrame(
+            {
+                "intensity": [0.0, 5.0, 10.0],
+                "damage": [0.0, 0.1, 0.8],
+            }
+        )
+    )
+    expected = PiecewiseLinearDamageCurve(
+        pandas.DataFrame(
+            {
+                "intensity": [0.0, 5.0, 10.0],
+                "damage": [0.0, 0.3, 0.9],
+            }
+        )
+    )
+    actual = PiecewiseLinearDamageCurve.interpolate(a, b, 0.5)
+    assert actual == expected

tests/test_intersection.py

@@ -13,6 +13,7 @@
     GridDefinition,
     split_linestrings,
     split_polygons,
+    generate_grid_boxes,
 )
 
 
@@ -196,6 +197,32 @@ def test_split_polygons(self, grid, polygon, polygon_split):
         assert_array_equal(actual["col1"].values, expected["col1"].values)
 
 
+def test_box_geom_bounds():
+    """Values take from tests/integration/range.tif"""
+    grid = GridDefinition(
+        crs=CRS.from_epsg(4326),
+        width=23,
+        height=14,
+        transform=(
+            0.008333333347826087,
+            0.0,
+            -1.341666667,
+            0.0,
+            -0.008333333285714315,
+            51.808333333,
+        ),
+    )
+    box_geoms = generate_grid_boxes(grid)
+    minb = box_geoms.bounds.min()
+    maxb = box_geoms.bounds.max()
+
+    atol = 1e-4
+    assert abs(minb.minx - -1.3416667) < atol
+    assert abs(minb.miny - 51.6916667) < atol
+    assert abs(maxb.maxx - -1.1500000) < atol
+    assert abs(maxb.maxy - 51.8083333) < atol
+
+
 def sort_polygons(df):
     iterations = 6  # all coords must be <= (2**p - 1) ; 2**6 - 1 == 63
     ndimensions = 2