Feature/relative matching distance

CHANGELOG.md

@@ -28,13 +28,17 @@ Removed:
 - Added optional parameter to `geo_im_from_array`, `plot_from_gdf`, `plot_rp_imp`, `plot_rp_intensity`,
 `plot_intensity`, `plot_fraction`, `_event_plot` to mask plotting when regions are too far from data points [#1047](https://github.com/CLIMADA-project/climada_python/pull/1047). To recreate previous plots (no masking), the parameter can be set to None.
 - Added instructions to install Climada petals on Euler cluster in `doc.guide.Guide_Euler.ipynb` [#1029](https://github.com/CLIMADA-project/climada_python/pull/1029)
-
+- Added util methods to handle crs coordinates consistently: `is_geo_coords`, `check_if_geo_coords`, `get_crs_unit`, `estimate_matching_threshold`, `degree_to_km`, and `km_to_degree` [#1080](https://github.com/CLIMADA-project/climada_python/pull/1080)
 - `ImpactFunc` and `ImpactFuncSet` now support equality comparisons via `==` [#1027](https://github.com/CLIMADA-project/climada_python/pull/1027)
 
 ### Changed
-
+- Changed default distance threshold for nearest neighbor matching in `util.coordinates.match_coordinates` from a fixed value of 100km to twice the highest resolution of the coords_to_assign [#1080](https://github.com/CLIMADA-project/climada_python/pull/1080).
 - Changed the default mask_distance in `util.plot.geo_im_from_array` to 0.03 to avoid white gaps in gridded hazard data with comparably low resolution (>80 centroids per axis) [#1073](https://github.com/CLIMADA-project/climada_python/pull/1073)
 - Increased speed of `util.plot.add_shapes` by avoiding for loops, substantially speeding up `Hazard.plot_intensity` and other functions. [#1073](https://github.com/CLIMADA-project/climada_python/pull/1073)
+- Update `util.coordinates.match_centroids`, `util.coordinates.match_coordinates`, so that they also
+accept coordinates that are not defined in degree. [#1080](https://github.com/CLIMADA-project/climada_python/pull/1080)
+- Implement cheap test to check that input coordinates at least seem geographic for functions that require
+geographic coordinates as input (e.g. `util.coordinates.dist_to_coast`, `util.coordinates.coord_on_land`, `util.coordinates.lon_normalize`, `util.coordinates.lon_bounds`). [#1080](https://github.com/CLIMADA-project/climada_python/pull/1080)
 - `Hazard.local_exceedance_intensity`, `Hazard.local_return_period` and `Impact.local_exceedance_impact`, `Impact.local_return_period`, using the `climada.util.interpolation` module: New default (no binning), binning on decimals, and faster implementation [#1012](https://github.com/CLIMADA-project/climada_python/pull/1012)
 - World Bank indicator data is now downloaded directly from their API via the function `download_world_bank_indicator`, instead of relying on the `pandas-datareader` package [#1033](https://github.com/CLIMADA-project/climada_python/pull/1033)
 - `Exposures.write_hdf5` pickles geometry data in WKB format, which is faster and more sustainable. [#1051](https://github.com/CLIMADA-project/climada_python/pull/1051)

climada/engine/impact_calc.py

@@ -92,6 +92,12 @@
             Default: True
         assign_centroids : bool, optional
             indicates whether centroids are assigned to the self.exposures object.
+            The default centroids assignment uses Euclidean distance and a threshold equal to
+            twice the highest resolution (smallest distance between any two nearest neighbours)
+            of the centroids (assuming a regular grid).
+            Recommendation: assign the centroids directly with control
+            over the distance metric and the distance threshold using
+            `exposure.assign_centroids(hazard)`. In this case, set assign_centroids to ``False``.
             Centroids assignment is an expensive operation; set this to ``False`` to save
             computation time if the hazards' centroids are already assigned to the exposures
             object.
@@ -114,12 +120,22 @@
             >>> imp = impcalc.impact(insured=True)
             >>> imp.aai_agg
 
+            >>> haz = Hazard.from_hdf5(HAZ_DEMO_H5)  # Set hazard
+            >>> impfset = ImpactFuncSet.from_excel(ENT_TEMPLATE_XLS)
+            >>> exp = Exposures(pd.read_excel(ENT_TEMPLATE_XLS))
+            >>> #Adjust threshold and distance to centroids/exposure resolution
+            >>> exp.assign_centroids(hazard, threshold=1.5, distance='euclidean')
+            >>> impcalc = ImpactCal(exp, impfset, haz)
+            >>> imp = impcalc.impact(insured=True, assign_centroids=False)
+            >>> imp.aai_agg
+
         See also
         --------
         apply_deductible_to_mat : apply deductible to impact matrix
         apply_cover_to_mat : apply cover to impact matrix
+        climada.entity.exposures.assign_centroids : assign centroids to exposures explicitly
         """
         # TODO: consider refactoring, making use of Exposures.hazard_impf
         # check for compatibility of exposures and hazard type
         if all(
             name not in self.exposures.gdf.columns

climada/engine/test/test_impact_calc.py

@@ -210,6 +210,28 @@ def test_calc_impact_TC_pass(self):
         self.assertAlmostEqual(6.570532945599105e11, impact.tot_value)
         self.assertAlmostEqual(6.512201157564421e09 * x, impact.aai_agg, 5)
 
+    def test_calc_impact_TC_change_crs_pass(self):
+        """Test compute impact"""
+        crs = "EPSG:4087"
+        exp = ENT.exposures.copy()
+        exp.to_crs(crs)
+        haz = deepcopy(HAZ)
+        haz.centroids.to_crs(crs)
+        icalc = ImpactCalc(exp, ENT.impact_funcs, haz)
+        impact = icalc.impact()
+        self.assertEqual(icalc.n_events, len(impact.at_event))
+        self.assertEqual(0, impact.at_event[0])
+        self.assertEqual(0, impact.at_event[7225])
+        self.assertAlmostEqual(1.472482938320243e08, impact.at_event[13809], delta=1)
+        self.assertAlmostEqual(7.076504723057620e10, impact.at_event[12147], delta=1)
+        self.assertEqual(0, impact.at_event[14449])
+        self.assertEqual(icalc.n_exp_pnt, len(impact.eai_exp))
+        self.assertAlmostEqual(1.518553670803242e08, impact.eai_exp[0], delta=1)
+        self.assertAlmostEqual(1.373490457046383e08, impact.eai_exp[25], 6)
+        self.assertAlmostEqual(1.066837260150042e08, impact.eai_exp[49], 6)
+        self.assertAlmostEqual(6.570532945599105e11, impact.tot_value)
+        self.assertAlmostEqual(6.512201157564421e09, impact.aai_agg, 5)
+
     def test_calc_impact_RF_pass(self):
         haz = Hazard.from_hdf5(get_test_file("test_hazard_US_flood_random_locations"))
         exp = Exposures.from_hdf5(
@@ -574,10 +596,13 @@ def test_single_exp_zero_mdr(self):
         eai_exp = np.array([aai_agg])
         at_event = np.array([imp_evt, 0])
         exp.set_geometry_points()
+        exp.assign_centroids(haz, threshold=1)
         impf_tc = ImpfTropCyclone.from_emanuel_usa()
         impf_set = ImpactFuncSet([impf_tc])
         impf_set.check()
-        imp = ImpactCalc(exp, impf_set, haz).impact(save_mat=True)
+        imp = ImpactCalc(exp, impf_set, haz).impact(
+            assign_centroids=False, save_mat=True
+        )
         check_impact(self, imp, haz, exp, aai_agg, eai_exp, at_event, at_event)
 
 

climada/entity/exposures/base.py

@@ -591,65 +591,81 @@
         self,
         hazard,
         distance="euclidean",
-        threshold=u_coord.NEAREST_NEIGHBOR_THRESHOLD,
+        threshold=None,
         overwrite=True,
     ):
-        """Assign for each exposure coordinate closest hazard coordinate.
+        """Assign for each exposure coordinate the closest hazard coordinate.
         The Exposures ``gdf`` will be altered by this method. It will have an additional
         (or modified) column named ``centr_[hazard.HAZ_TYPE]`` after the call.
 
         Uses the utility function ``u_coord.match_centroids``. See there for details
         and parameters.
 
         The value -1 is used for distances larger than ``threshold`` in point distances.
-        In case of raster hazards the value -1 is used for centroids outside of the raster.
 
         Parameters
         ----------
         hazard : Hazard
-            Hazard to match (with raster or vector centroids).
+            Hazard to match. The hazard centroids must be in the same coordinate
+            reference system (crs) as the exposures.
         distance : str, optional
-            Distance to use in case of vector centroids.
+            Distance to use for matching exposures points to hazard centroids
             Possible values are "euclidean", "haversine" and "approx".
             Default: "euclidean"
         threshold : float
-            If the distance (in km) to the nearest neighbor exceeds `threshold`,
+            Assignement threshold in units of the exposure's crs.
+            If the distance to the nearest neighbor exceeds `threshold`,
             the index `-1` is assigned.
-            Set `threshold` to 0, to disable nearest neighbor matching.
-            Default: 100 (km)
+            Set `threshold` to 0, to disable nearest neighbor matching and enforce
+            exact matching.
+            Default: twice the highest resolution (smallest distance between any two nearest neighbours) of the hazard centroids
         overwrite: bool
             If True, overwrite centroids already present. If False, do
             not assign new centroids. Default is True.
 
         See Also
         --------
-        climada.util.coordinates.match_grid_points: method to associate centroids to
-            exposure points when centroids is a raster
         climada.util.coordinates.match_coordinates:
             method to associate centroids to exposure points
+        climada.util.coordinates.estimate_matching_threshold:
+            method to calculate the default threshold
+        climada.util.coordinates.km_to_degree:
+            method to approximately convert kilometers to degrees
+        climada.util.coordinates.degree_to_km:
+            method to approximately convert degrees to kilometers
+
         Notes
         -----
-        The default order of use is:
-
-        1. if centroid raster is defined, assign exposures points to
-           the closest raster point.
-        2. if no raster, assign centroids to the nearest neighbor using
-           euclidian metric
-
-        Both cases can introduce innacuracies for coordinates in lat/lon
-        coordinates as distances in degrees differ from distances in meters
-        on the Earth surface, in particular for higher latitude and distances
-        larger than 100km. If more accuracy is needed, please use 'haversine'
-        distance metric. This however is slower for (quasi-)gridded data,
-        and works only for non-gridded data.
+        For coordinates in lat/lon coordinates distances in degrees differ from
+        distances in meters on the Earth surface, in particular for higher
+        latitude and distances larger than 100km. If more accuracy for degree
+        coordinates is needed, please use 'haversine' distance metric,
+        which however is slower.
+
+        Caution: nearest neighbourg matching can introduce serious artefacts
+        such as:
+            - exposure and hazard centroids with shifted grids can lead
+            to systematically wrong assignements.
+            - hazard centroids covering larger areas than exposures may lead
+            to sub-optimal matching if the threshold is too large
+            - projected crs often diverge at the anti-meridian and close points
+            on either side will be at a large distance. For proper handling
+            of the anti-meridian please use degree coordinates in EPSG:4326.
+            This might be relevant for countries like the Fidji or the US that
+            cross the anti-meridian.
+
+        Users are free to implement their own matching alrogithm and save the
+        matching centroid index in the appropriate column ``centr_[hazard.HAZ_TYPE]``.
         """
         haz_type = hazard.haz_type
         centr_haz = INDICATOR_CENTR + haz_type
         if centr_haz in self.gdf:
-            LOGGER.info("Exposures matching centroids already found for %s", haz_type)
             if overwrite:
                 LOGGER.info("Existing centroids will be overwritten for %s", haz_type)
             else:
+                LOGGER.info(
+                    "Exposures matching centroids already found for %s", haz_type
+                )
                 return
 
         LOGGER.info(
@@ -659,7 +675,7 @@
         )
 
         if not u_coord.equal_crs(self.crs, hazard.centroids.crs):
-            raise ValueError("Set hazard and exposure to same CRS first!")
+            raise ValueError("Set hazard and exposure to the same CRS first!")
         # Note: equal_crs is tested here, rather than within match_centroids(),
         # because exp.gdf.crs may not be defined, but exp.crs must be defined.
 

climada/entity/exposures/test/test_base.py

@@ -103,6 +103,59 @@ def test_assign_pass(self):
             self.assertEqual(exp.gdf.shape[0], len(exp.hazard_centroids("FL")))
             np.testing.assert_array_equal(exp.hazard_centroids("FL"), expected_result)
 
+    def test_assign_meter_crs_pass(self):
+        """Check that attribute `assigned` is correctly set."""
+        np_rand = np.random.RandomState(123456789)
+
+        haz = Hazard.from_raster(
+            [HAZ_DEMO_FL], haz_type="FL", window=Window(10, 20, 50, 60)
+        )
+        ncentroids = haz.centroids.size
+        crs_meters = "EPSG:4087"
+        haz.centroids.to_crs(crs_meters, inplace=True)
+        exp = Exposures(
+            crs=haz.centroids.crs,
+            lon=np.concatenate(
+                [
+                    haz.centroids.lon,
+                    haz.centroids.lon + 0.001 * (-0.5 + np_rand.rand(ncentroids)),
+                ]
+            ),
+            lat=np.concatenate(
+                [
+                    haz.centroids.lat,
+                    haz.centroids.lat + 0.001 * (-0.5 + np_rand.rand(ncentroids)),
+                ]
+            ),
+        )
+        expected_result = np.concatenate([np.arange(ncentroids), np.arange(ncentroids)])
+
+        # make sure that it works for both float32 and float64
+        for test_dtype in [np.float64, np.float32]:
+            haz.centroids.gdf["lat"] = haz.centroids.lat.astype(test_dtype)
+            haz.centroids.gdf["lon"] = haz.centroids.lon.astype(test_dtype)
+            exp.assign_centroids(haz)
+            self.assertEqual(exp.gdf.shape[0], len(exp.hazard_centroids("FL")))
+            np.testing.assert_array_equal(exp.hazard_centroids("FL"), expected_result)
+
+    def test_assign_no_equal_crs(self):
+        haz = Hazard.from_raster(
+            [HAZ_DEMO_FL], haz_type="FL", window=Window(10, 20, 50, 60)
+        )
+        crs_km = "EPSG:4326"
+        crs_meters = "EPSG:4087"
+        exp = Exposures(
+            crs=crs_meters,
+            lon=haz.centroids.lon,
+            lat=haz.centroids.lat,
+        )
+        haz.centroids.to_crs(crs_km, inplace=True)
+        with self.assertRaises(ValueError) as cm:
+            exp.assign_centroids(haz)
+        self.assertEqual(
+            "Set hazard and exposure to the same CRS first!", str(cm.exception)
+        )
+
     def test__init__meta_type(self):
         """Check if meta of type list raises a ValueError in __init__"""
         with self.assertRaises(TypeError) as cm:
@@ -212,7 +265,7 @@ def test_assign_raster_pass(self):
                 0.0,
             ],
         )
-        exp.assign_centroids(haz)
+        exp.assign_centroids(haz, threshold=u_coord.km_to_degree(100))
 
         expected_result = [
             # constant y-value, varying x-value

climada/hazard/test/test_tc_tracks.py

@@ -1350,7 +1350,7 @@ def test_track_land_params(self):
         lon_test = np.array([170, 179.18, 180.05])
         lat_test = np.array([-60, -16.56, -16.85])
         on_land = np.array([False, True, True])
-        lon_shift = np.array([-360, 0, 360])
+        lon_shift = np.array([360, 360, 360])
         # ensure both points are considered on land as is
         np.testing.assert_array_equal(
             u_coord.coord_on_land(lat=lat_test, lon=lon_test), on_land

climada/hazard/trop_cyclone/trop_cyclone_windfields.py

@@ -902,6 +902,8 @@ def _coriolis_parameter(lat: np.ndarray) -> np.ndarray:
     cp : np.ndarray of same shape as input
         Coriolis parameter.
     """
+    u_coord.check_if_geo_coords(lat, 0)
+
     return 2 * V_ANG_EARTH * np.sin(np.radians(np.abs(lat)))