Merge branch 'develop' into feature/climada_boario

CHANGELOG.md

@@ -18,9 +18,12 @@ Code freeze date: YYYY-MM-DD
 
 ### Changed
 
+- Adaptations to refactoring of the `climada.hazard.Centroids` class, to be compatible with `climada>=5.0` [#122](https://github.com/CLIMADA-project/climada_petals/pull/122)
+
 ### Fixed
 
 - Fix `climada.hazard.tc_rainfield` for TC tracks crossing the antimeridian [#105](https://github.com/CLIMADA-project/climada_petals/pull/105)
+- Update the table of content for the tutorials [#125](https://github.com/CLIMADA-project/climada_petals/pull/125)
 
 ### Deprecated
 

climada_petals/hazard/emulator/geo.py

@@ -27,6 +27,7 @@
 from shapely.geometry import Polygon
 
 from climada.hazard import Centroids
+from climada.util.constants import NATEARTH_CENTROIDS
 import climada.util.coordinates as u_coord
 import climada_petals.hazard.emulator.const as const
 
@@ -102,8 +103,8 @@ def centroids(self, latlon=None, res_as=360):
         centroids : climada.hazard.Centroids object
         """
         if latlon is None:
-            centroids = Centroids.from_base_grid(res_as=res_as)
-            centroids.set_meta_to_lat_lon()
+            # default centroids: Natural Earth data at given resolution
+            centroids = Centroids.from_hdf5(NATEARTH_CENTROIDS[res_as])
         else:
             centroids = Centroids.from_lat_lon(*latlon)
         msk = shapely.vectorized.contains(self.shape, centroids.lon, centroids.lat)

climada_petals/hazard/landslide.py

@@ -299,8 +299,9 @@ def from_prob(cls, bbox, path_sourcefile, corr_fact=10e6, n_years=500,
 
         haz = cls()
         # raster with occurrence probs
-        haz.centroids.meta, prob_matrix = \
+        meta, prob_matrix = \
             u_coord.read_raster(path_sourcefile, geometry=[shapely.geometry.box(*bbox, ccw=True)])
+        haz.centroids = Centroids.from_meta(meta)
         prob_matrix = prob_matrix.squeeze()/corr_fact
 
         # sample events from probabilities
@@ -314,11 +315,6 @@ def from_prob(cls, bbox, path_sourcefile, corr_fact=10e6, n_years=500,
         haz.event_name = np.array(range(n_years))
         haz.event_id = np.array(range(n_years))
 
-        if not haz.centroids.meta['crs'].is_epsg_code:
-            haz.centroids.meta['crs'] = haz.centroids.meta['crs'
-                               ].from_user_input(DEF_CRS)
-        haz.centroids.set_geometry_points()
-
         haz.check()
 
         return haz

climada_petals/hazard/low_flow.py

@@ -539,11 +539,7 @@ def _centroids_resolution(centroids):
         -------
         float
         """
-        if centroids.meta:
-            res_centr = abs(centroids.meta['transform'][4]), \
-                        centroids.meta['transform'][0]
-        else:
-            res_centr = np.abs(u_coord.get_resolution(centroids.lat, centroids.lon))
+        res_centr = np.abs(u_coord.get_resolution(centroids.lat, centroids.lon))
         if np.abs(res_centr[0] - res_centr[1]) > 1.0e-6:
             LOGGER.warning('Centroids do not represent regular pixels %s.', str(res_centr))
             return (res_centr[0] + res_centr[1]) / 2
@@ -645,7 +641,6 @@ def _init_centroids(dis_xarray, centr_res_factor=1):
         (dis_xarray.lon.values.min(), dis_xarray.lat.values.min(),
          dis_xarray.lon.values.max(), dis_xarray.lat.values.max()),
         res=res_data / centr_res_factor)
-    centroids.set_meta_to_lat_lon()
     centroids.set_area_approx()
     centroids.set_on_land()
     centroids.empty_geometry_points()

climada_petals/hazard/relative_cropyield.py

@@ -246,8 +246,7 @@ def from_isimip_netcdf(cls, input_dir=None, filename=None, bbox=None,
         haz.units = 't / y / ha'
         haz.date = np.array(dt.str_to_date(
             [event_ + '-01-01' for event_ in haz.event_name]))
-        haz.centroids.set_meta_to_lat_lon()
-        haz.centroids.region_id = (
+        haz.centroids.gdf['region_id'] = (
             coord.coord_on_land(haz.centroids.lat, haz.centroids.lon)).astype(dtype=int)
         haz.check()
         return haz
@@ -371,8 +370,6 @@ def plot_time_series(self, event=None):
         else:
             event = [str(n) for n in range(event[0], event[1] + 1)]
 
-        self.centroids.set_meta_to_lat_lon()
-
         # definition of plot extents
         len_lat = abs(self.centroids.lat[0] - self.centroids.lat[-1]) * (2.5 / 13.5)
         len_lon = abs(self.centroids.lon[0] - self.centroids.lon[-1]) * (5 / 26)

climada_petals/hazard/river_flood.py

@@ -23,7 +23,6 @@
 
 import logging
 import datetime as dt
-import copy
 from collections.abc import Iterable
 from pathlib import Path
 
@@ -137,22 +136,22 @@ def from_nc(cls, dph_path=None, frc_path=None, origin=False,
             if ISINatIDGrid:
 
                 dest_centroids = RiverFlood._select_exact_area(countries, reg)[0]
-                meta_centroids = copy.copy(dest_centroids)
-                meta_centroids.set_lat_lon_to_meta()
+                centroids_meta = dest_centroids.get_meta()
 
                 haz = cls.from_raster(files_intensity=[dph_path],
                                       files_fraction=[frc_path], band=bands.tolist(),
-                                      transform=meta_centroids.meta['transform'],
-                                      width=meta_centroids.meta['width'],
-                                      height=meta_centroids.meta['height'],
+                                      transform=centroids_meta['transform'],
+                                      width=centroids_meta['width'],
+                                      height=centroids_meta['height'],
                                       resampling=Resampling.nearest)
-                x_i = ((dest_centroids.lon - haz.centroids.meta['transform'][2]) /
-                       haz.centroids.meta['transform'][0]).astype(int)
-                y_i = ((dest_centroids.lat - haz.centroids.meta['transform'][5]) /
-                       haz.centroids.meta['transform'][4]).astype(int)
+                haz_centroids_meta = haz.centroids.get_meta()
+                x_i = ((dest_centroids.lon - haz_centroids_meta['transform'][2]) /
+                       haz_centroids_meta['transform'][0]).astype(int)
+                y_i = ((dest_centroids.lat - haz_centroids_meta['transform'][5]) /
+                       haz_centroids_meta['transform'][4]).astype(int)
 
-                fraction = haz.fraction[:, y_i * haz.centroids.meta['width'] + x_i]
-                intensity = haz.intensity[:, y_i * haz.centroids.meta['width'] + x_i]
+                fraction = haz.fraction[:, y_i * haz_centroids_meta['width'] + x_i]
+                intensity = haz.intensity[:, y_i * haz_centroids_meta['width'] + x_i]
 
                 haz.centroids = dest_centroids
                 haz.intensity = sp.sparse.csr_matrix(intensity)
@@ -166,14 +165,12 @@ def from_nc(cls, dph_path=None, frc_path=None, origin=False,
                                           files_fraction=[frc_path],
                                           band=bands.tolist(),
                                           geometry=cntry_geom.geoms)
-                    # self.centroids.set_meta_to_lat_lon()
                 else:
                     cntry_geom = u_coord.get_land_geometry(countries)
                     haz = cls.from_raster(files_intensity=[dph_path],
                                           files_fraction=[frc_path],
                                           band=bands.tolist(),
                                           geometry=cntry_geom.geoms)
-                    # self.centroids.set_meta_to_lat_lon()
 
         elif shape:
             shapes = gpd.read_file(shape)
@@ -195,7 +192,6 @@ def from_nc(cls, dph_path=None, frc_path=None, origin=False,
             haz = cls.from_raster(files_intensity=[dph_path],
                                   files_fraction=[frc_path],
                                   band=bands.tolist())
-            # self.centroids.set_meta_to_lat_lon()
 
         else:  # use given centroids
             # if centroids.meta or grid_is_regular(centroids)[0]:
@@ -206,8 +202,6 @@ def from_nc(cls, dph_path=None, frc_path=None, origin=False,
             #      (transform)
             #      reprojection change resampling"""
             # else:
-            if centroids.meta:
-                centroids.set_meta_to_lat_lon()
             metafrc, fraction = u_coord.read_raster(frc_path, band=bands.tolist())
             metaint, intensity = u_coord.read_raster(dph_path, band=bands.tolist())
             x_i = ((centroids.lon - metafrc['transform'][2]) /
@@ -321,7 +315,7 @@ def set_flooded_area(self, save_centr=False):
         MemoryError
         """
         self.centroids.set_area_pixel()
-        area_centr = self.centroids.area_pixel
+        area_centr = self.centroids.get_area_pixel()
         event_years = np.array([dt.date.fromordinal(self.date[i]).year
                                 for i in range(len(self.date))])
         years = np.unique(event_years)

climada_petals/hazard/tc_rainfield.py

@@ -229,7 +229,7 @@ def set_from_tracks(self, *args, **kwargs):
     def from_tracks(
         cls,
         tracks: TCTracks,
-        centroids: Optional[Centroids] = None,
+        centroids: Centroids = None,
         pool: Optional[pathos.pools.ProcessPool] = None,
         model: str = 'R-CLIPER',
         model_kwargs: Optional[dict] = None,
@@ -298,7 +298,8 @@ def from_tracks(
         tracks : climada.hazard.TCTracks
             Tracks of storm events.
         centroids : Centroids, optional
-            Centroids where to model TC. Default: global centroids at 360 arc-seconds resolution.
+            Centroids where to model TC. Default: centroids at 360 arc-seconds resolution within
+            tracks' bounds.
         pool : pathos.pool, optional
             Pool that will be used for parallel computation of rain fields. Default: None
         model : str, optional
@@ -398,20 +399,15 @@ def from_tracks(
         """
         num_tracks = tracks.size
         if centroids is None:
-            centroids = Centroids.from_base_grid(res_as=360, land=True)
-
-        if not centroids.coord.size:
-            centroids.set_meta_to_lat_lon()
+            centroids = Centroids.from_pnt_bounds(tracks.get_bounds(), res=0.1)
 
         if ignore_distance_to_coast:
             # Select centroids with lat <= max_latitude
             [idx_centr_filter] = (np.abs(centroids.lat) <= max_latitude).nonzero()
         else:
             # Select centroids which are inside max_dist_inland_km and lat <= max_latitude
-            if not centroids.dist_coast.size:
-                centroids.set_dist_coast()
             [idx_centr_filter] = (
-                (centroids.dist_coast <= max_dist_inland_km * 1000)
+                (centroids.get_dist_coast() <= max_dist_inland_km * 1000)
                 & (np.abs(centroids.lat) <= max_latitude)
             ).nonzero()
 

climada_petals/hazard/tc_surge_bathtub.py

@@ -72,15 +72,11 @@ def from_tc_winds(wind_haz, topo_path, inland_decay_rate=0.2, add_sea_level_rise
         """
         centroids = copy.deepcopy(wind_haz.centroids)
 
-        if not centroids.coord.size:
-            centroids.set_meta_to_lat_lon()
-
         # Select wind-affected centroids which are inside MAX_DIST_COAST and |lat| < 61
-        if not centroids.dist_coast.size or np.all(centroids.dist_coast >= 0):
-            centroids.set_dist_coast(signed=True, precomputed=True)
+        centroids_dist_coast = centroids.get_dist_coast(signed=True)
         coastal_msk = (wind_haz.intensity > 0).sum(axis=0).A1 > 0
-        coastal_msk &= (centroids.dist_coast < 0)
-        coastal_msk &= (centroids.dist_coast >= -MAX_DIST_COAST * 1000)
+        coastal_msk &= (centroids_dist_coast < 0)
+        coastal_msk &= (centroids_dist_coast >= -MAX_DIST_COAST * 1000)
         coastal_msk &= (np.abs(centroids.lat) <= MAX_LATITUDE)
 
         # Load elevation at coastal centroids
@@ -117,7 +113,7 @@ def from_tc_winds(wind_haz, topo_path, inland_decay_rate=0.2, add_sea_level_rise
 
         if inland_decay_rate != 0:
             # Add decay according to distance from coast
-            dist_coast_km = np.abs(centroids.dist_coast[coastal_idx]) / 1000
+            dist_coast_km = np.abs(centroids_dist_coast[coastal_idx]) / 1000
             coastal_centroids_h += inland_decay_rate * dist_coast_km
         coastal_centroids_h -= add_sea_level_rise
 
@@ -168,12 +164,9 @@ def _fraction_on_land(centroids, topo_path):
     """
     bounds = np.array(centroids.total_bounds)
     shape = [0, 0]
-    if centroids.meta:
-        shape = centroids.shape
-        cen_trans = centroids.meta['transform']
-    else:
-        shape[0], shape[1], cen_trans = u_coord.pts_to_raster_meta(
-            bounds, min(u_coord.get_resolution(centroids.lat, centroids.lon)))
+    shape[0], shape[1], cen_trans = u_coord.pts_to_raster_meta(
+        points_bounds=bounds,
+        res=min(u_coord.get_resolution(centroids.lat, centroids.lon)))
 
     read_raster_buffer = 0.5 * max(np.abs(cen_trans[0]), np.abs(cen_trans[4]))
     bounds += read_raster_buffer * np.array([-1., -1., 1., 1.])
@@ -191,9 +184,4 @@ def _fraction_on_land(centroids, topo_path):
                             resampling=rasterio.warp.Resampling.average,
                             src_nodata=dem_nodata, dst_nodata=0.0)
 
-    if not centroids.meta:
-        x_i = ((centroids.lon - cen_trans[2]) / cen_trans[0]).astype(int)
-        y_i = ((centroids.lat - cen_trans[5]) / cen_trans[4]).astype(int)
-        fractions = fractions[y_i, x_i]
-
     return fractions.reshape(-1)

climada_petals/hazard/test/test_tc_rainfield.py

@@ -20,7 +20,6 @@
 """
 
 import datetime as dt
-from pathlib import Path
 import unittest
 
 import numpy as np
@@ -45,12 +44,12 @@
 
 def getTestData():
     client = Client()
-    centr_ds = client.get_dataset_info(name='test_tc_rainfield', status='test_dataset')
-    _, [centr_test_mat, track, track_short, haz_mat] = client.download_dataset(centr_ds)
-    return Centroids.from_mat(centr_test_mat), track, track_short, haz_mat
+    centr_ds = client.get_dataset_info(name='tc_rainfield_test', status='test_dataset')
+    _, [centr_test_hdf5, track, track_short, haz_hdf5] = client.download_dataset(centr_ds)
+    return Centroids.from_hdf5(centr_test_hdf5), track, track_short, haz_hdf5
 
 
-CENTR_TEST_BRB, TEST_TRACK, TEST_TRACK_SHORT, HAZ_TEST_MAT = getTestData()
+CENTR_TEST_BRB, TEST_TRACK, TEST_TRACK_SHORT, HAZ_TEST_HDF5 = getTestData()
 
 
 def tcrain_examples():
@@ -117,7 +116,6 @@ def test_cross_antimeridian(self):
         # Two locations on the island Taveuni (Fiji), one west and one east of 180° longitude.
         # We list the second point twice, with different lon-normalization:
         cen = Centroids.from_lat_lon([-16.95, -16.8, -16.8], [179.9, 180.1, -179.9])
-        cen.set_dist_coast(precomputed=True)
 
         # Cyclone YASA (2020) passed directly over Fiji
         tr = TCTracks.from_ibtracs_netcdf(storm_id=["2020346S13168"])

climada_petals/hazard/test/test_tc_surge_bathtub.py

@@ -49,8 +49,7 @@ def __init__(self, bounds, res_deg, slope=0.006):
         self.shape = (lat.size, lon.size)
         self.transform = rasterio.Affine(res_deg, 0, bounds[0], 0, -res_deg, bounds[3])
         centroids = Centroids.from_lat_lon(*[ar.ravel() for ar in np.meshgrid(lon, lat)][::-1])
-        centroids.set_dist_coast(signed=True, precomputed=True)
-        self.dist_coast = centroids.dist_coast
+        self.dist_coast = centroids.get_dist_coast(signed=True)
         self.slope = slope
 
     def __enter__(self):
@@ -95,14 +94,13 @@ def test_fraction_on_land(self):
         shape = (lat.size, lon.size)
         lon, lat = [ar.ravel() for ar in np.meshgrid(lon, lat)]
         centroids = Centroids.from_lat_lon(lat, lon)
-        centroids.set_dist_coast(signed=True, precomputed=True)
-
+
         dem_bounds = (bounds[0] - 1, bounds[1] - 1, bounds[2] + 1, bounds[3] + 1)
         dem_res = 3 / (60 * 60)
         with tmp_artifical_topo(dem_bounds, dem_res) as topo_path:
             fraction = _fraction_on_land(centroids, topo_path)
         fraction = fraction.reshape(shape)
-        dist_coast = centroids.dist_coast.reshape(shape)
+        dist_coast = centroids.get_dist_coast(signed=True).reshape(shape)
 
         # check valid range and order of magnitude
         self.assertTrue(np.all((fraction >= 0) & (fraction <= 1)))
@@ -155,7 +153,6 @@ def test_surge_from_track(self):
         shape = (lat.size, lon.size)
         lon, lat = [ar.ravel() for ar in np.meshgrid(lon, lat)]
         centroids = Centroids.from_lat_lon(lat, lon)
-        centroids.set_dist_coast(signed=True, precomputed=True)
 
         wind_haz = TropCyclone.from_tracks(tc_track, centroids=centroids)
 
@@ -182,7 +179,6 @@ def test_cross_antimeridian(self):
         # Two locations on the island Taveuni (Fiji), one west and one east of 180° longitude.
         # We list the second point twice, with different lon-normalization:
         cen = Centroids.from_lat_lon([-16.95, -16.8, -16.8], [179.9, 180.1, -179.9])
-        cen.set_dist_coast(precomputed=True)
 
         # Cyclone YASA (2020) passed directly over Fiji
         tr = TCTracks.from_ibtracs_netcdf(storm_id=["2020346S13168"])

climada_petals/hazard/test/test_wildfire.py

@@ -273,21 +273,26 @@ def test_centroids_pass(self):
         """ Test _firms_centroids_creation """
         wf = WildFire()
         firms = wf._clean_firms_df(TEST_FIRMS)
-        centroids = wf._firms_centroids_creation(firms, 0.375/ONE_LAT_KM, 1/2)
-        self.assertEqual(centroids.meta['width'], 144)
-        self.assertEqual(centroids.meta['height'], 138)
-        self.assertAlmostEqual(centroids.meta['transform'][0], 0.006749460043196544)
-        self.assertAlmostEqual(centroids.meta['transform'][1], 0.0)
-        self.assertTrue(centroids.meta['transform'][2]<= firms.longitude.min())
-        self.assertAlmostEqual(centroids.meta['transform'][3], 0.0)
-        self.assertAlmostEqual(centroids.meta['transform'][4], -centroids.meta['transform'][0])
-        self.assertTrue(centroids.meta['transform'][5] >= firms.latitude.max())
-        self.assertTrue(firms.latitude.max() <= centroids.total_bounds[3])
-        self.assertTrue(firms.latitude.min() >= centroids.total_bounds[1])
-        self.assertTrue(firms.longitude.max() <= centroids.total_bounds[2])
-        self.assertTrue(firms.longitude.min() >= centroids.total_bounds[0])
+        res_data = 0.375/ONE_LAT_KM
+        centr_res_factor = 1/2
+        centroids = wf._firms_centroids_creation(firms, res_data, centr_res_factor)
+        centroids_meta = centroids.get_meta()
+        self.assertEqual(centroids_meta['width'], 144)
+        self.assertEqual(centroids_meta['height'], 138)
+        self.assertAlmostEqual(centroids_meta['transform'][0], 0.006749460043196544)
+        self.assertAlmostEqual(centroids_meta['transform'][1], 0.0)
+        self.assertLessEqual(centroids_meta['transform'][2], firms.longitude.min())
+        self.assertAlmostEqual(centroids_meta['transform'][3], 0.0)
+        self.assertAlmostEqual(centroids_meta['transform'][4], -centroids_meta['transform'][0])
+        self.assertGreaterEqual(centroids_meta['transform'][5], firms.latitude.max())
+        # Meta pixel coordinates are shifted by half resolution to the pixel center lat/lon
+        pixel_center_shift = (res_data / centr_res_factor) / 2
+        self.assertLessEqual(firms.latitude.max(), centroids.total_bounds[3] + pixel_center_shift)
+        self.assertGreaterEqual(firms.latitude.min(), centroids.total_bounds[1] - pixel_center_shift)
+        self.assertLessEqual(firms.longitude.max(), centroids.total_bounds[2] + pixel_center_shift)
+        self.assertGreaterEqual(firms.longitude.min(), centroids.total_bounds[0] - pixel_center_shift)
         self.assertTrue(centroids.lat.size)
-        self.assertTrue(centroids.area_pixel.size)
+        self.assertTrue(centroids.get_area_pixel().size)
         self.assertTrue(centroids.on_land.size)
 
     def test_firms_resolution_pass(self):