CLIMADA-project · emanuel-schmid · May 13, 2024 · Feb 22, 2024 · Feb 22, 2024 · Feb 23, 2024
diff --git a/Makefile b/Makefile
@@ -16,15 +16,15 @@ help:  ## Use one of the following instructions:
 
 .PHONY : lint
 lint : ## Static code analysis with Pylint
-	pylint -ry climada_petals > pylint.log || true
+	python -m pylint -ry climada_petals > pylint.log || true
 
 .PHONY : unit_test
 unit_test : ## Unit tests execution with coverage and xml reports
-	pytest $(PYTEST_ARGS) --ignore=climada_petals/test climada_petals/
+	python -m pytest $(PYTEST_ARGS) --ignore=climada_petals/test climada_petals/
 
 .PHONY : install_test
 install_test : ## Test installation was successful
-	pytest $(PYTEST_JUNIT_ARGS) --pyargs climada.engine.test.test_cost_benefit \
+	python -m pytest $(PYTEST_JUNIT_ARGS) --pyargs climada.engine.test.test_cost_benefit \
 	climada.engine.test.test_impact
 
 .PHONY : data_test
@@ -37,11 +37,11 @@ notebook_test : ## Test notebooks in doc/tutorial
 
 .PHONY : integ_test
 integ_test : ## Integration tests execution with xml reports
-	pytest $(PYTEST_ARGS) climada_petals/test/
+	python -m pytest $(PYTEST_ARGS) climada_petals/test/
 
 .PHONY : test
 test : ## Unit and integration tests execution with coverage and xml reports
-	pytest $(PYTEST_ARGS) climada_petals/
+	python -m pytest $(PYTEST_ARGS) climada_petals/
 
 .PHONY : ci-clean
 ci-clean :

diff --git a/climada_petals/hazard/emulator/geo.py b/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)

diff --git a/climada_petals/hazard/landslide.py b/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

diff --git a/climada_petals/hazard/low_flow.py b/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()

diff --git a/climada_petals/hazard/relative_cropyield.py b/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)

diff --git a/climada_petals/hazard/river_flood.py b/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 @@
             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 @@
                                           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 @@
             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 @@
             #      (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 @@
         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)

diff --git a/climada_petals/hazard/tc_rainfield.py b/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,
         pool: Optional[pathos.pools.ProcessPool] = None,
         model: str = 'R-CLIPER',
         model_kwargs: Optional[dict] = None,
@@ -297,7 +297,7 @@ def from_tracks(
         ----------
         tracks : climada.hazard.TCTracks
             Tracks of storm events.
-        centroids : Centroids, optional
+        centroids : Centroids
             Centroids where to model TC. Default: global centroids at 360 arc-seconds resolution.
         pool : pathos.pool, optional
             Pool that will be used for parallel computation of rain fields. Default: None
@@ -397,21 +397,13 @@ def from_tracks(
         TCRain
         """
         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()
-
         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()
 

diff --git a/climada_petals/hazard/tc_surge_bathtub.py b/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)
diff --git a/climada_petals/hazard/test/test_tc_rainfield.py b/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_mat, track, track_short, haz_hdf5] = client.download_dataset(centr_ds)
+    return Centroids.from_hdf5(centr_test_mat), 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():
@@ -67,7 +66,7 @@ def test_set_one_pass(self):
         """Test from_tracks constructor with a single track."""
         tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK)
         tc_track.equal_timestep()
-        tc_haz = TCRain.from_tracks(tc_track, centroids=CENTR_TEST_BRB)
+        tc_haz = TCRain.from_tracks(tc_track, CENTR_TEST_BRB)
 
         self.assertEqual(tc_haz.haz_type, 'TR')
         self.assertEqual(tc_haz.units, 'mm')
@@ -95,7 +94,7 @@ def test_tcr(self):
         tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK)
         tc_track.equal_timestep()
 
-        tc_haz = TCRain.from_tracks(tc_track, model="TCR", centroids=CENTR_TEST_BRB)
+        tc_haz = TCRain.from_tracks(tc_track, CENTR_TEST_BRB, model="TCR")
         self.assertTrue(isinstance(tc_haz.intensity, sparse.csr_matrix))
         self.assertEqual(tc_haz.intensity.shape, (1, 296))
         self.assertEqual(tc_haz.intensity.nonzero()[0].size, 296)
@@ -106,7 +105,7 @@ def test_tcr(self):
         # This increases the results by more than 70% because the default value for saturation
         # specific humidity corresponds to a temperature of only ~267 K.
         tc_track.data[0]["t600"] = xr.full_like(tc_track.data[0]["central_pressure"], 275.0)
-        tc_haz = TCRain.from_tracks(tc_track, model="TCR", centroids=CENTR_TEST_BRB)
+        tc_haz = TCRain.from_tracks(tc_track, CENTR_TEST_BRB, model="TCR")
         self.assertTrue(isinstance(tc_haz.intensity, sparse.csr_matrix))
         self.assertEqual(tc_haz.intensity.shape, (1, 296))
         self.assertEqual(tc_haz.intensity.nonzero()[0].size, 296)
@@ -117,12 +116,11 @@ 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"])
 
-        inten = TCRain.from_tracks(tr, centroids=cen).intensity.toarray()[0, :]
+        inten = TCRain.from_tracks(tr, cen).intensity.toarray()[0, :]
 
         # Centroids 1 and 2 are identical, they just use a different normalization for lon. This
         # should not affect the result at all:
@@ -135,7 +133,7 @@ def test_cross_antimeridian(self):
     def test_from_file_pass(self):
         """Test from_tracks constructor with one input."""
         tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK_SHORT)
-        tc_haz = TCRain.from_tracks(tc_track, centroids=CENTR_TEST_BRB)
+        tc_haz = TCRain.from_tracks(tc_track, CENTR_TEST_BRB)
         tc_haz.check()
 
         self.assertEqual(tc_haz.haz_type, 'TR')
@@ -159,7 +157,7 @@ def test_from_file_pass(self):
     def test_two_files_pass(self):
         """Test from_tracks constructor with two ibtracs."""
         tc_track = TCTracks.from_processed_ibtracs_csv([TEST_TRACK_SHORT, TEST_TRACK_SHORT])
-        tc_haz = TCRain.from_tracks(tc_track, centroids=CENTR_TEST_BRB)
+        tc_haz = TCRain.from_tracks(tc_track, CENTR_TEST_BRB)
         tc_haz.remove_duplicates()
         tc_haz.check()
 
@@ -206,14 +204,15 @@ def test_compute_rain_pass(self):
     def test_rainfield_diff_time_steps(self):
         """Check that the results do not depend too much on the track's time step sizes."""
         tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK)
+        centroids = CENTR_TEST_BRB
 
-        train_org = TCRain.from_tracks(tc_track)
+        train_org = TCRain.from_tracks(tc_track, centroids)
 
         tc_track.equal_timestep(time_step_h=1)
-        train_1h = TCRain.from_tracks(tc_track)
+        train_1h = TCRain.from_tracks(tc_track, centroids)
 
         tc_track.equal_timestep(time_step_h=0.5)
-        train_05h = TCRain.from_tracks(tc_track)
+        train_05h = TCRain.from_tracks(tc_track, centroids)
 
         for train in [train_1h, train_05h]:
             np.testing.assert_allclose(