Don't re-project projected polygon gdfs in metre-disaggregation

CHANGELOG.md

@@ -18,9 +18,10 @@ Removed:
 
 - `Impact.impact_at_reg` method for aggregating impacts per country or custom region [#642](https://github.com/CLIMADA-project/climada_python/pull/642)
 
-### Changed
+### Changed 
 
 ### Fixed
+- `util.lines_polys_handler` solve polygon disaggregation issue in metre-based projection   [#666](https://github.com/CLIMADA-project/climada_python/pull/666)
 
 ### Deprecated
 

climada/util/lines_polys_handler.py

@@ -38,8 +38,8 @@ class DisaggMethod(Enum):
     """
     Disaggregation Method for the ... function
 
-	DIV : the geometry's distributed to equal parts over all its interpolated points
-	FIX : the geometry's value is replicated over all its interpolated points
+    DIV : the geometry's distributed to equal parts over all its interpolated points
+    FIX : the geometry's value is replicated over all its interpolated points
     """
     DIV = 'div'
     FIX = 'fix'
@@ -49,7 +49,7 @@ class AggMethod(Enum):
     """
     Aggregation Method for the aggregate_impact_mat function
 
-	SUM : the impact is summed over all points in the polygon/line
+    SUM : the impact is summed over all points in the polygon/line
     """
     SUM = 'sum'
 
@@ -656,7 +656,8 @@ def _poly_to_pnts(gdf, res, to_meters):
 
     gdf_points = gdf.copy().reset_index(drop=True)
 
-    if to_meters:
+    # Check if we need to reproject
+    if to_meters and not gdf.geometry.crs.is_projected:
         gdf_points['geometry_pnt'] = gdf_points.apply(
             lambda row: _interp_one_poly_m(row.geometry, res, gdf.crs), axis=1)
     else:
@@ -969,7 +970,9 @@ def _pnts_per_line(length, res):
 
 def _swap_geom_cols(gdf, geom_to, new_geom):
     """
-    Change which column is the geometry column
+    Change which column is the geometry column.
+    Conserves the crs of the original geometry column.
+
     Parameters
     ----------
     gdf : gpd.GeoDataFrame
@@ -986,7 +989,7 @@ def _swap_geom_cols(gdf, geom_to, new_geom):
     """
     gdf_swap = gdf.rename(columns = {'geometry': geom_to})
     gdf_swap.rename(columns = {new_geom: 'geometry'}, inplace=True)
-    gdf_swap.set_geometry('geometry', inplace=True)
+    gdf_swap.set_geometry('geometry', inplace=True, crs=gdf.crs)
     return gdf_swap
 
 

climada/util/test/test_lines_polys_handler.py

@@ -20,16 +20,17 @@
 """
 
 import unittest
-from pathlib import Path
+from unittest.mock import patch, DEFAULT
 
 import numpy as np
 
+from shapely.geometry import Point
+
 from climada.entity import Exposures
 import climada.util.lines_polys_handler as u_lp
 import climada.util.coordinates as u_coord
 from climada.util.api_client import Client
 from climada.engine import Impact, ImpactCalc
-from climada.hazard import Hazard
 from climada.entity.impact_funcs import ImpactFuncSet
 from climada.entity.impact_funcs.storm_europe import ImpfStormEurope
 
@@ -144,6 +145,60 @@ def test_point_exposure_from_polygons(self):
             ])
         np.testing.assert_allclose(exp_pnt.gdf.latitude, lat)
 
+        #projected crs, to_meters=TRUE, FIX, dissag_val
+        res = 20000
+        EXP_POLY_PROJ = Exposures(GDF_POLY.to_crs(epsg=28992))
+        exp_pnt = u_lp.exp_geom_to_pnt(
+            EXP_POLY_PROJ, res=res, to_meters=True,
+            disagg_met=u_lp.DisaggMethod.FIX, disagg_val=res**2
+            )
+        self.check_unchanged_exp(EXP_POLY_PROJ, exp_pnt)
+        val = res**2
+        self.assertEqual(np.unique(exp_pnt.gdf.value)[0], val)
+        self.assertEqual(exp_pnt.gdf.crs, EXP_POLY_PROJ.gdf.crs)
+
+    @patch.multiple(
+        "climada.util.lines_polys_handler",
+        _interp_one_poly=DEFAULT,
+        _interp_one_poly_m=DEFAULT,
+    )
+    def test_point_exposure_from_polygons_reproject_call(
+        self, _interp_one_poly, _interp_one_poly_m
+    ):
+        """Verify that the correct subroutine is called for a reprojected CRS"""
+        # Just have the mocks return an empty geometry
+        _interp_one_poly.return_value = Point(1.0, 1.0)
+        _interp_one_poly_m.return_value = Point(1.0, 1.0)
+
+        # Use geographical CRS
+        EXP_POLY_PROJ = Exposures(GDF_POLY.to_crs(epsg=4326))
+        res = 1
+        u_lp.exp_geom_to_pnt(
+            EXP_POLY_PROJ,
+            res=res,
+            to_meters=True,
+            disagg_met=u_lp.DisaggMethod.FIX,
+            disagg_val=res,
+        )
+        _interp_one_poly_m.assert_called()
+        _interp_one_poly.assert_not_called()
+
+        # Reset mocks
+        _interp_one_poly.reset_mock()
+        _interp_one_poly_m.reset_mock()
+
+        # Use projected CRS
+        EXP_POLY_PROJ = Exposures(GDF_POLY.to_crs(epsg=28992))
+        u_lp.exp_geom_to_pnt(
+            EXP_POLY_PROJ,
+            res=res,
+            to_meters=True,
+            disagg_met=u_lp.DisaggMethod.FIX,
+            disagg_val=res,
+        )
+        _interp_one_poly_m.assert_not_called()
+        _interp_one_poly.assert_called()
+
     def test_point_exposure_from_polygons_on_grid(self):
         """Test disaggregation of polygons to points on grid"""
         exp_poly = EXP_POLY.copy()