Add functions to compute, plot, store the local hazard exceedence intensity and RP maps

climada/hazard/test/test_base.py

@@ -1011,6 +1011,24 @@ def test_ref_all_pass(self):
         self.assertAlmostEqual(inten_stats[3][33], 88.510983305123631)
         self.assertAlmostEqual(inten_stats[2][99], 79.717518054203623)
 
+    def test_local_return_period(self):
+        """Compare local return periods against reference."""
+        haz = dummy_hazard()
+        haz.intensity = sparse.csr_matrix([[1.1, 2.1, 3.1],
+                                   [0.1, 3.1, 2.1],
+                                   [4.1, 1.1, 0.1],
+                                   [2.1, 3.1, 3.1]])
+        haz.frequency = np.full(4, .5)
+        threshold_intensities = np.array([1., 2., 4.])
+        return_stats = haz.local_return_period(threshold_intensities)
+        np.testing.assert_allclose(
+            return_stats,
+            np.array([
+                [0.66666667, 0.5, 0.66666667],
+                [1., 0.66666667, 0.66666667],
+                [2., np.nan, np.nan]
+                ])
+        )
 
 class TestYearset(unittest.TestCase):
     """Test return period statistics"""

climada/test/test_hazard.py

@@ -22,8 +22,12 @@
 import unittest
 import numpy as np
 import datetime as dt
+import os
+from tempfile import TemporaryDirectory
 from pathlib import Path
 from scipy import sparse
+from rioxarray import open_rasterio
+from xarray import load_dataset
 
 from climada import CONFIG
 from climada.hazard import tc_tracks as tc
@@ -188,6 +192,47 @@ def test_read_write_vector_fraction_pass(self):
 
         DATA_DIR.joinpath(intensity_file).unlink()
         DATA_DIR.joinpath(fraction_file).unlink()
+
+    def test_write_raster_local_exceedance_inten(self):
+        """Test write TIFF file from local exceedance intensity"""
+        haz = Hazard.from_hdf5(HAZ_TEST_TC)
+        haz.write_raster_local_exceedance_inten([10, 20, 30, 40], filename = f'{CONFIG.test_data.dir()}/test_write_raster_local_exceedance_inten.tif')
+        dataarray = open_rasterio(f'{CONFIG.test_data.dir()}/test_write_raster_local_exceedance_inten.tif')
+        np.testing.assert_array_almost_equal(
+            dataarray.data[[0, 1, 2, 3], [0, 8, 3, 6], [2, 5, 7, 1]],
+            np.array([35.829193, 37.64515 , 61.62324 , 58.595936])
+        )
+        np.testing.assert_array_equal(dataarray.data.shape, (4, 10, 10))
+
+    def test_write_raster_local_return_periods(self):
+        """Test write TIFF file from local exceedance intensity"""
+        self.temp_dir = TemporaryDirectory()
+        self.test_file_path = os.path.join(self.temp_dir.name, 'test_file.tif')
+
+        haz = Hazard.from_hdf5(HAZ_TEST_TC)
+        haz.write_raster_local_return_periods([10., 20., 30.], filename = self.test_file_path)
+        #haz.write_raster_local_return_periods([10., 20., 30.], filename = f'{CONFIG.test_data.dir()}/test_write_raster_local_return_periods.tif')
+        #dataarray = open_rasterio(f'{CONFIG.test_data.dir()}/test_write_raster_local_return_periods.tif')
+        dataarray = open_rasterio(self.test_file_path)
+
+        np.testing.assert_array_almost_equal(
+            dataarray.data[[0, 1, 2], [2, 5, 8], [1, 6, 7]],
+            np.array([2.42105263, 2.59677419, 4.01496259])
+        )
+        np.testing.assert_array_equal(dataarray.data.shape, (3, 10, 10))
+        self.temp_dir.cleanup()
+
+    def test_write_netcdf_local_return_periods(self):
+        """Test write netcdf file from local exceedance intensity"""
+        haz = Hazard.from_hdf5(HAZ_TEST_TC)
+        haz.write_netcdf_local_return_periods([10., 20., 30.], filename = f'{CONFIG.test_data.dir()}/test_write_netcdf_local_return_periods.nc')
+        dataset = load_dataset(f'{CONFIG.test_data.dir()}/test_write_netcdf_local_return_periods.nc')
+        np.testing.assert_array_almost_equal(
+            dataset.to_array().data[[0, 1, 2, 3, 4], [21, 45, 86, 65, 7]],
+            np.array([22., -80. ,1.483871, 3.108108, 5.366667])
+        )
+        np.testing.assert_array_equal(dataset.to_array().data.shape, (5, 100))
+
 
 
 

climada/test/test_plot.py

@@ -105,6 +105,14 @@ def test_hazard_rp_intensity(self):
         (axis1, axis2), _ = hazard.plot_rp_intensity([25, 50])
         self.assertEqual('Return period: 25 years', axis1.get_title())
         self.assertEqual('Return period: 50 years', axis2.get_title())
+
+    def test_plot_local_rp(self):
+        """"Plot local return period maps for different hazard intensities"""
+        hazard = Hazard.from_hdf5(HAZ_TEST_TC)
+        (axis1, axis2), return_stats = hazard.plot_local_rp([10., 20.])
+        self.assertEqual('Intensity: 10.0 m/s', axis1.get_title())
+        self.assertEqual('Intensity: 20.0 m/s', axis2.get_title())
+        np.testing.assert_array_equal(return_stats.shape, (2, 100))
 
     def test_exposures_value_pass(self):
         """Plot exposures values."""