Feature/centroids as gdf

Makefile

@@ -2,6 +2,9 @@
 # test, coverage and lint
 ###
 
+# Default pytest command
+PYTEST_CMD := pytest
+
 PYTEST_JUNIT_ARGS = --junitxml=tests_xml/tests.xml
 
 PYTEST_COV_ARGS = \
@@ -16,15 +19,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/
+	$(PYTEST_CMD) $(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 \
+	$(PYTEST_CMD) $(PYTEST_JUNIT_ARGS) --pyargs climada.engine.test.test_cost_benefit \
 	climada.engine.test.test_impact
 
 .PHONY : data_test
@@ -37,11 +40,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/
+	$(PYTEST_CMD) $(PYTEST_ARGS) climada_petals/test/
 
 .PHONY : test
 test : ## Unit and integration tests execution with coverage and xml reports
-	pytest $(PYTEST_ARGS) climada_petals/
+	$(PYTEST_CMD) $(PYTEST_ARGS) climada_petals/
 
 .PHONY : ci-clean
 ci-clean :

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 @@
             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)

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()
 

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)