Feature/coastal flooding

climada_petals/conf/climada.conf

@@ -55,9 +55,17 @@
                 }
             }
         },
+        "coastal_flood": {
+            "local_data" : {
+                "aqueduct" : "{local_data.system}/CoastalFlood/Aqueduct"
+            },
+            "resources" : {
+                "aqueduct" : "http://wri-projects.s3.amazonaws.com/AqueductFloodTool/download/v2/"
+            },
         "copernicus": {
             "local_data": "{local_data.system}/copernicus_data",
             "seasonal_forecasts": "{local_data.system}/copernicus_data/seasonal_forecasts"
+            }
         }
     }
 }

climada_petals/hazard/__init__.py

@@ -22,3 +22,4 @@
 from .tc_rainfield import *
 from .tc_surge_bathtub import *
 from .wildfire import *
+from .coastal_flood import *

climada_petals/hazard/coastal_flood.py

@@ -0,0 +1,130 @@
+from typing import Iterable, Union, Optional
+import numpy as np
+from shapely.geometry import Polygon
+
+from climada.hazard import Hazard
+from climada.util import files_handler as u_fh
+import climada.util.coordinates as u_coord
+from climada import CONFIG
+
+AQUEDUCT_SOURCE_LINK = CONFIG.hazard.coastal_flood.resources.aqueduct.str()
+DOWNLOAD_DIRECTORY = CONFIG.hazard.coastal_flood.local_data.aqueduct.dir()
+
+HAZ_TYPE = 'CF'
+"""Hazard type acronym CoastalFlood"""
+
+__all__ = ['CoastalFlood']
+
+class CoastalFlood(Hazard):
+    """
+    Contains coastal flood events pulled by the Acqueduct project :
+        https://www.wri.org/data/aqueduct-floods-hazard-maps
+    """
+
+    def __init__(self, **kwargs):
+        """Empty constructor"""
+
+        Hazard.__init__(self, HAZ_TYPE, **kwargs)
+
+    @classmethod
+    def from_aqueduct_tif(cls,
+                          rcp: str,
+                          target_year : str,
+                          return_periods: Union[int, Iterable[int]],
+                          subsidence: Optional[str] = 'wtsub',
+                          percentile: str = '95',
+                          countries: Optional[Union[str, Iterable[str]]] = None,
+                          boundaries: Iterable[float] = None,
+                          dwd_dir=DOWNLOAD_DIRECTORY):
+        """
+        Download and load coastal flooding events from the Aqueduct dataset.
+        For more info on the data see: 
+            https://files.wri.org/d8/s3fs-public/aqueduct-floods-methodology.pdf
+
+        Parameters
+        ----------
+        rcp : str
+            RCPs scenario. Possible values are historical, 45 and 85.
+        target_year : str
+            future target year. Possible values are hist, 2030, 2050 and 2080.
+        return_periods : int or list of int
+            events' return periods.
+            Possible values are 2, 5, 10, 25, 50, 100, 250, 500, 1000.
+        subsidence : str
+            If land subsidence is simulated or not.
+            Possible values are nosub and wtsub. Default wtsub.
+        percentile : str
+            Sea level rise scenario (in percentile).
+            Possible values are 05, 50 and 95.
+        countries : str or list of str
+            countries ISO3 codes
+        boundaries : tuple of floats
+            geographical boundaries in the order:
+                minimum longitude, minimum latitude,
+                maximum longitude, maximum latitude
+        Returns
+        -------
+        haz : CoastalFlood instance
+        """
+
+        if (target_year == 'hist') & (rcp != 'historical'):
+            raise ValueError(f"RCP{rcp} cannot have hist as target year")
+
+        if (rcp == 'historical') & (subsidence == 'nosub') & (target_year != 'hist'):
+            raise ValueError("Historical without subsidence can only have hist as target year")
+
+        if (rcp == 'historical') & (percentile != '0'):
+            raise ValueError("Historical shall not have a assigned percentiles of sea level rise. Leave default values.")
+
+        if isinstance(return_periods, int):
+            return_periods = [return_periods]
+
+        return_periods.sort(reverse=True)
+
+        if isinstance(countries, str):
+            countries = [countries]
+
+        rcp_name = f"rcp{rcp[0]}p{rcp[1]}" if rcp in ['45', '85'] else rcp
+        ret_per_names = [f"{str(rp).zfill(4)}" for rp in return_periods]
+        perc_name = f"0_perc_{percentile.zfill(2)}" if percentile in ['5', '50'] else '0'
+
+        file_names = [f"inuncoast_{rcp_name}_{subsidence}_{target_year}_"
+                      f"rp{rp.zfill(4)}_{perc_name}.tif"
+                      for rp in ret_per_names]
+
+        file_paths = []
+        for file_name in file_names:
+            link_to_file = "".join([AQUEDUCT_SOURCE_LINK, file_name])
+            file_paths.append(dwd_dir / file_name)
+
+            if not file_paths[-1].exists():
+                u_fh.download_file(link_to_file, download_dir=dwd_dir)
+
+        if countries:
+            # TODO: this actually does not work with multiple countries
+            geom = u_coord.get_land_geometry(countries).geoms
+
+        elif boundaries:
+            min_lon, min_lat, max_lon, max_lat = boundaries
+            geom = [Polygon([(min_lon, min_lat),
+                            (max_lon, min_lat),
+                            (max_lon, max_lat),
+                            (min_lon, max_lat)])]
+        else:
+            geom = None
+
+        event_id = np.arange(len(file_names))
+        frequencies = np.diff(1 / np.array(return_periods), prepend=0)
+        event_names = [f"1-in-{return_periods[i]}y_{percentile}pct_"
+                       f"{rcp_name}_{target_year}_{subsidence}"
+                       for i in range(len(file_names))]
+
+        haz = cls().from_raster(files_intensity=file_paths,
+                                geometry=geom,
+                                attrs={'event_id': event_id,
+                                       'event_name': event_names,
+                                       'frequency': frequencies})
+
+        haz.units = 'm'
+
+        return haz

climada_petals/hazard/test/test_coastal_flood.py

@@ -0,0 +1,66 @@
+
+"""
+This file is part of CLIMADA.
+
+Copyright (C) 2017 ETH Zurich, CLIMADA contributors listed in AUTHORS.
+
+CLIMADA is free software: you can redistribute it and/or modify it under the
+terms of the GNU General Public License as published by the Free
+Software Foundation, version 3.
+
+CLIMADA is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with CLIMADA. If not, see <https://www.gnu.org/licenses/>.
+
+---
+
+Tests on Coastal Flood Hazard"""
+
+import unittest
+import requests
+import numpy as np
+from climada_petals.hazard.coastal_flood import AQUEDUCT_SOURCE_LINK
+
+RCPS = ['historical', 'rcp4p5', 'rcp8p5']
+TARGET_YEARS = ['hist', '2030', '2050', '2080']
+RETURN_PERIODS = ['0002', '0005', '0010', '0025',
+                  '0050', '0100', '0250', '0500', '1000']
+SUBSIDENCE = ['nosub', 'wtsub']
+PERCENTILES = ['0_perc_05', '0_perc_50', '0']
+
+class TestReader(unittest.TestCase):
+    """Test that Coastal flood data exist"""
+
+    def test_files_exist(self):
+
+        file_names = [
+            f'inuncoast_{rcp}_{sub}_{year}_rp{rp}_{perc}.tif'
+                for rcp in RCPS
+                for sub in SUBSIDENCE
+                for year in TARGET_YEARS
+                for rp in RETURN_PERIODS
+                for perc in PERCENTILES
+                # You can't have:
+                # - year historic with rcp different than historical
+                # - rcp historical, no subsidence and year different than historic
+                # - rcp historical and SLR scenarios' percentiles
+                if not (((year == 'hist') & (rcp != 'historical')) |
+                        ((rcp == 'historical') & (sub == 'nosub') & (year != 'hist')) |
+                        ((rcp == 'historical') & (perc != '0')))
+                ]
+
+        test_files_pos = np.random.choice(range(len(file_names)),
+                                          size=10,
+                                          replace=False)
+        for i in test_files_pos:
+            file_path = "".join([AQUEDUCT_SOURCE_LINK, file_names[i]])
+            request_code = requests.get(file_path).status_code
+            self.assertTrue(request_code == 200)
+
+# Execute Tests
+if __name__ == "__main__":
+    TESTS = unittest.TestLoader().loadTestsFromTestCase(TestReader)
+    unittest.TextTestRunner(verbosity=2).run(TESTS)