Add functions to compute, plot, store the local hazard exceedence intensity and RP maps (new branch)

climada/hazard/base.py

@@ -450,6 +450,46 @@ def sanitize_event_ids(self):
             LOGGER.debug('Resetting event_id.')
             self.event_id = np.arange(1, self.event_id.size + 1)
 
+    def local_return_period(self, threshold_intensities):
+        """Compute local return periods for given hazard intensities.
+
+        Parameters
+        ----------
+        threshold_intensities : np.array
+            Hazard intensities to consider.
+
+        Returns
+        -------
+        return_periods : np.array
+            Array containing computed local return periods for given hazard intensities.
+        """
+        # Ensure threshold_intensities is a numpy array
+        threshold_intensities = np.array(threshold_intensities)
+
+        num_cen = self.intensity.shape[1]
+        return_periods = np.zeros((len(threshold_intensities), num_cen))  # Adjusted for 2D structure
+
+        # Process each centroid in chunks as in local_exceedance_inten
+        cen_step = CONFIG.max_matrix_size.int() // self.intensity.shape[0]
+        if not cen_step:
+            raise ValueError('Increase max_matrix_size configuration parameter to >'
+                             f'{self.intensity.shape[0]}')
+
+        chk = -1
+        for chk in range(int(num_cen / cen_step)):
+            self._loc_return_period(
+                threshold_intensities,
+                self.intensity[:, chk * cen_step:(chk + 1) * cen_step].toarray(),
+                return_periods[:, chk * cen_step:(chk + 1) * cen_step])
+
+        if (chk + 1) * cen_step < num_cen:  # Check if there's a remainder
+            self._loc_return_period(
+                threshold_intensities,
+                self.intensity[:, (chk + 1) * cen_step:].toarray(),
+                return_periods[:, (chk + 1) * cen_step:])
+
+        return return_periods
+
     def get_event_id(self, event_name):
         """Get an event id from its name. Several events might have the same
         name.
@@ -613,6 +653,44 @@ def _loc_return_inten(self, return_periods, inten, exc_inten):
             exc_inten[:, cen_idx] = self._cen_return_inten(
                 inten_sort[:, cen_idx], freq_sort[:, cen_idx],
                 self.intensity_thres, return_periods)
+
+    def _loc_return_period(self, threshold_intensities, inten, return_periods):
+        """Compute local return periods for given hazard intensities for a specific chunk of data.
+
+        Parameters
+        ----------
+        threshold_intensities: np.array
+            Given hazard intensities for which to calculate return periods.
+        inten: np.array
+            The intensity array for a specific chunk of data.
+        return_periods: np.array
+            Array to store computed return periods for the given hazard intensities.
+        """
+        # Assuming inten is sorted and calculating cumulative frequency
+        sort_pos = np.argsort(inten, axis=0)[::-1, :]
+        inten_sort = inten[sort_pos, np.arange(inten.shape[1])]
+        freq_sort = self.frequency[sort_pos]
+        np.cumsum(freq_sort, axis=0, out=freq_sort)
+
+        for cen_idx in range(inten.shape[1]):
+            sorted_inten_cen = inten_sort[:, cen_idx]
+            cum_freq_cen = freq_sort[:, cen_idx]
+
+            for i, intensity in enumerate(threshold_intensities):
+                # Find the first occurrence where the intensity is less than the sorted intensities
+                exceedance_index = np.searchsorted(sorted_inten_cen[::-1], intensity, side='right')
+
+                # Calculate exceedance probability
+                if exceedance_index < len(cum_freq_cen):
+                    exceedance_probability = cum_freq_cen[-exceedance_index - 1]
+                else:
+                    exceedance_probability = 0  # Or set a default minimal probability
+
+                # Calculate and store return period
+                if exceedance_probability > 0:
+                    return_periods[i, cen_idx] = 1 / exceedance_probability
+                else:
+                    return_periods[i, cen_idx] = np.nan  
 
     def _check_events(self):
         """Check that all attributes but centroids contain consistent data.

climada/hazard/io.py

@@ -27,6 +27,7 @@
 from typing import Union, Optional, Callable, Dict, Any
 
 import h5py
+import netCDF4 as nc
 import numpy as np
 import pandas as pd
 import rasterio
@@ -1038,6 +1039,128 @@ def write_hdf5(self, file_name, todense=False):
                         )
         self.centroids.write_hdf5(file_name, mode='a')
 
+    def write_raster_local_exceedance_inten(self, return_periods, filename):
+        """
+        Generates exceedance intensity data for specified return periods and 
+        saves it into a GeoTIFF file.
+
+        Parameters
+        ----------
+        return_periods : np.array or list
+            Array or list of return periods (in years) for which to calculate 
+            and store exceedance intensities.
+        filename : str
+            Path and name of the file to write in tif format.
+        """
+        inten_stats = self.local_exceedance_inten(return_periods=return_periods)
+        num_bands = inten_stats.shape[0]
+
+        if not self.centroids.get_meta():
+            raise ValueError("centroids.get_meta() is required but not set.")
+
+        ### this code is to replace pixel_geom = self.centroids.calc_pixels_polygons()
+        if abs(abs(self.centroids.get_meta()['transform'].a) -
+               abs(self.centroids.get_meta()['transform'].e)) > 1.0e-5:
+            raise ValueError('Area can not be computed for not squared pixels.')
+        pixel_geom = self.centroids.geometry.buffer(self.centroids.get_meta()['transform'].a / 2).envelope
+        ###
+        profile = self.centroids.get_meta().copy()
+        profile.update(driver='GTiff', dtype='float32', count=num_bands)
+
+        with rasterio.open(filename, 'w', **profile) as dst:
+            LOGGER.info('Writing %s', filename)
+            for band in range(num_bands):
+                raster = rasterio.features.rasterize(
+                    [(x, val) for (x, val) in zip(pixel_geom, inten_stats[band].reshape(-1))],
+                    out_shape=(profile['height'], profile['width']),
+                    transform=profile['transform'], fill=0,
+                    all_touched=True, dtype=profile['dtype'])
+                dst.write(raster, band + 1)
+
+                band_name = f"Exceedance intensity for RP {return_periods[band]} years"
+                dst.set_band_description(band + 1, band_name)   
+
+    def write_raster_local_return_periods(self, threshold_intensities, filename, output_resolution=None):
+        """Write local return periods map as GeoTIFF file.
+
+        Parameters
+        ----------
+        threshold_intensities: np.array
+            Hazard intensities to consider for calculating return periods.
+        filename: str
+            File name to write in tif format.
+        output_resolution: int
+            If not None, the data is rasterized to this resolution.
+            Default is None (resolution is estimated from the data).
+        """
+
+        # Calculate the local return periods
+        variable = self.local_return_period(threshold_intensities)
+
+        # Obtain the meta information for the raster file
+        meta = self.centroids.get_meta(resolution=output_resolution)
+        meta.update(driver='GTiff', dtype=rasterio.float32, count=len(threshold_intensities))
+        res = meta["transform"][0]  # resolution from lon coordinates
+
+        if meta['height'] * meta['width'] == self.centroids.size:
+            # centroids already in raster format
+            u_coord.write_raster(filename, variable, meta)
+        else:
+            geometry = self.centroids.get_pixel_shapes(res=res)
+            with rasterio.open(filename, 'w', **meta) as dst:
+                LOGGER.info('Writing %s', filename)
+                for i_ev in range(len(threshold_intensities)):
+                    raster = rasterio.features.rasterize(
+                        (
+                            (geom, value)
+                            for geom, value
+                            in zip(geometry, variable[i_ev].flatten())
+                        ),
+                        out_shape=(meta['height'], meta['width']),
+                        transform=meta['transform'],
+                        fill=0,
+                        all_touched=True,
+                        dtype=meta['dtype'],
+                    )
+                    dst.write(raster.astype(meta['dtype']), i_ev + 1)
+
+                    # Set the band description
+                    band_name = f"RP of intensity {threshold_intensities[i_ev]} {self.units}"
+                    dst.set_band_description(i_ev + 1, band_name)
+
+
+    def write_netcdf_local_return_periods(self, threshold_intensities, filename):
+        """Generates local return period data and saves it into a NetCDF file.
+
+        Parameters
+        ----------
+        threshold_intensities: np.array
+            Hazard intensities to consider for calculating return periods.
+        filename: str
+            Path and name of the file to write the NetCDF data.
+        """
+        return_periods = self.local_return_period(threshold_intensities)
+        coords = self.centroids.coord
+
+        with nc.Dataset(filename, 'w', format='NETCDF4') as dataset:
+            centroids_dim = dataset.createDimension('centroids', coords.shape[0])
+
+            latitudes = dataset.createVariable('latitude', 'f4', ('centroids',))
+            longitudes = dataset.createVariable('longitude', 'f4', ('centroids',))
+            latitudes[:] = coords[:, 0]
+            longitudes[:] = coords[:, 1]
+            latitudes.units = 'degrees_north'
+            longitudes.units = 'degrees_east'
+
+            for i, intensity in enumerate(threshold_intensities):
+                dataset_name = f'return_period_intensity_{int(intensity)}'
+                return_period_var = dataset.createVariable(dataset_name, 'f4', ('centroids',))
+                return_period_var[:] = return_periods[i, :]
+                return_period_var.units = 'years'
+                return_period_var.description = f'Local return period for hazard intensity {intensity} {self.units}'
+
+            dataset.description = 'Local return period data for given hazard intensities'
+
     def read_hdf5(self, *args, **kwargs):
         """This function is deprecated, use Hazard.from_hdf5."""
         LOGGER.warning("The use of Hazard.read_hdf5 is deprecated."

climada/hazard/plot.py

@@ -19,9 +19,11 @@
 Define Hazard Plotting Methods.
 """
 
+import geopandas as gpd
 import numpy as np
 import matplotlib.pyplot as plt
 
+import climada.util.coordinates as u_coord
 import climada.util.plot as u_plot
 
 
@@ -65,6 +67,45 @@ def plot_rp_intensity(self, return_periods=(25, 50, 100, 250),
                                         colbar_name, title, smooth=smooth, axes=axis,
                                         figsize=figsize, adapt_fontsize=adapt_fontsize, **kwargs)
         return axis, inten_stats
+
+    def plot_local_rp(self, threshold_intensities, smooth=True, axis=None, figsize=(9, 13), adapt_fontsize=True, cmap = 'viridis_r', **kwargs):
+        """Plot hazard local return periods for given hazard intensities.
+
+        Parameters
+        ----------
+        threshold_intensities: np.array
+            Hazard intensities to consider for calculating return periods.
+        smooth: bool, optional
+            Smooth plot to plot.RESOLUTION x plot.RESOLUTION.
+        axis: matplotlib.axes._subplots.AxesSubplot, optional
+            Axis to use.
+        figsize: tuple, optional
+            Figure size for plt.subplots.
+        kwargs: optional
+            Arguments for pcolormesh matplotlib function used in event plots.
+
+        Returns
+        -------
+        axis: matplotlib.axes._subplots.AxesSubplot
+            Matplotlib axis with the plot.
+        """
+        ### code to replace self._set_coords_centroids()
+        if self.centroids.get_meta() and not self.centroids.coord.size:
+            xgrid, ygrid = u_coord.raster_to_meshgrid(
+                self.centroids.get_meta()['transform'], self.centroids.get_meta()['width'], self.centroids.get_meta()['height'])
+            self.centroids.lon = xgrid.flatten()
+            self.centroids.lat = ygrid.flatten()
+            self.centroids.geometry = gpd.GeoSeries(crs=self.centroids.get_meta()['crs'])
+        ###
+        return_periods = self.local_return_period(threshold_intensities)
+        colbar_name = 'Return Period (years)'
+        title = list()
+        for haz_int in threshold_intensities:
+            title.append('Intensity: ' + f'{haz_int} {self.units}')
+        axis = u_plot.geo_im_from_array(return_periods, self.centroids.coord,
+                                        colbar_name, title, smooth=smooth, axes=axis,
+                                        figsize=figsize, adapt_fontsize=adapt_fontsize, cmap=cmap, **kwargs)
+        return axis, return_periods
 
     def plot_intensity(self, event=None, centr=None, smooth=True, axis=None, adapt_fontsize=True,
                        **kwargs):

climada/hazard/test/test_base.py

@@ -1022,6 +1022,25 @@ def test_ref_all_pass(self):
         self.assertAlmostEqual(inten_stats[1][66], 70.608592953031405)
         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):