Merge branch 'develop' into cross-calibrate-impact-functions

Author: emanuel-schmid

.github/pull_request_template.md

@@ -19,10 +19,10 @@ This PR fixes #
 ### PR Reviewer Checklist
 
 - [ ] Read the [Contribution Guide][contrib]
-- [ ] [CLIMADA Reviewer Checklist](https://climada-python.readthedocs.io/en/latest/guide/Guide_Reviewer_Checklist.html) passed
+- [ ] [CLIMADA Reviewer Checklist](https://climada-python.readthedocs.io/en/latest/development/Guide_Review.html) passed
 - [ ] [Tests][testing] passing
 - [ ] No new [linter issues][linter]
 
 [contrib]: https://github.com/CLIMADA-project/climada_python/blob/main/CONTRIBUTING.md
-[testing]: https://climada-python.readthedocs.io/en/latest/guide/Guide_Continuous_Integration_and_Testing.html
-[linter]: https://climada-python.readthedocs.io/en/stable/guide/Guide_Continuous_Integration_and_Testing.html#3.C.--Static-Code-Analysis
+[testing]: https://climada-python.readthedocs.io/en/latest/development/Guide_continuous_integration_GitHub_actions.html
+[linter]: https://climada-python.readthedocs.io/en/latest/development/Guide_continuous_integration_GitHub_actions.html#static-code-analysis

CHANGELOG.md

@@ -28,19 +28,25 @@ 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)
 - Calibration of impact function ensembles in `climada.util.calibrate` [#1048](https://github.com/CLIMADA-project/climada_python/pull/1048)
 
 ### 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)
 - The online documentation has been completely overhauled, now uses PyData theme: [#977](https://github.com/CLIMADA-project/climada_python/pull/977)
 - `Input` to impact function calibration tasks now supports adding weights to the data [#1048](https://github.com/CLIMADA-project/climada_python/pull/1048)
 - Add `climada.hazard.xarray` module with helper structures for reading Hazard objects from `xarray` data [#1063](https://github.com/CLIMADA-project/climada_python/pull/1063)
+- The output of the `impact_yearset` was changed to only contain attributes corresponding to the yearly impact set. The application of the correction factor and the frequency of the resulting yearly impact object are corrected. [#1075](https://github.com/CLIMADA-project/climada_python/pull/1075)
 
 ### Fixed
 
@@ -106,8 +112,12 @@ Removed:
 
 ### Added
 
+- `climada.entity.impact_funcs.trop_cyclone.ImpfSetTropCyclone.get_impf_id_regions_per_countries` function [#1034](https://github.com/CLIMADA-project/climada_python/pull/1034)
+- `climada.hazard.tc_tracks.BasinBoundsStorm` Enum class `climada.hazard.tc_tracks.subset_by_basin` function [#1031](https://github.com/CLIMADA-project/climada_python/pull/1031)
 - `climada.hazard.tc_tracks.TCTracks.subset_years` function [#1023](https://github.com/CLIMADA-project/climada_python/pull/1023)
-- `climada.hazard.tc_tracks.TCTracks.from_FAST` function, add Australia basin (AU) [#993](https://github.com/CLIMADA-project/climada_python/pull/993)
+-`climada.hazard.tc_tracks.compute_track_density` function, `climada.hazard.tc_tracks.compute_genesis_density` function, `climada.hazard.plot.plot_track_density` function
+ [#1003](https://github.com/CLIMADA-project/climada_python/pull/1003)
+-`climada.hazard.tc_tracks.TCTracks.from_FAST` function, add Australia basin (AU) [#993](https://github.com/CLIMADA-project/climada_python/pull/993)
 - Add `osm-flex` package to CLIMADA core [#981](https://github.com/CLIMADA-project/climada_python/pull/981)
 - `doc.tutorial.climada_entity_Exposures_osm.ipynb` tutorial explaining how to use `osm-flex` with CLIMADA
 - `climada.util.coordinates.bounding_box_global` function [#980](https://github.com/CLIMADA-project/climada_python/pull/980)

climada/engine/impact_calc.py

@@ -92,6 +92,12 @@ def impact(
             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,10 +120,20 @@ def impact(
             >>> 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

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 @@ def assign_centroids(
         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 @@ def assign_centroids(
         )
 
         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