diff --git a/CHANGELOG.md b/CHANGELOG.md
index c2cc8f96de..467b61b960 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,18 +10,62 @@ Code freeze date: YYYY-MM-DD
### Dependency Changes
-### Added
-
### Changed
+- Centroids complete overhaul. Most function should be backward compatible. Internal data is stored in a geodataframe attribute. Raster are now stored as points, and the meta attribute is removed. Several methds were deprecated or removed. [#787](https://github.com/CLIMADA-project/climada_python/pull/787)
- Improved error messages produced by `ImpactCalc.impact()` in case impact function in the exposures is not found in impf_set [#863](https://github.com/CLIMADA-project/climada_python/pull/863)
### Fixed
+### Added
+
+- climada.hazard.centroids.centr.Centroids.get_area_pixel
+- climada.hazard.centroids.centr.Centroids.get_dist_coast
+- climada.hazard.centroids.centr.Centroids.get_elevation
+- climada.hazard.centroids.centr.Centroids.get_meta
+- climada.hazard.centroids.centr.Centroids.get_pixel_shapes
+- climada.hazard.centroids.centr.Centroids.to_crs
+- climada.hazard.centroids.centr.Centroids.to_default_crs
+- climada.hazard.centroids.centr.Centroids.write_csv
+- climada.hazard.centroids.centr.Centroids.write_excel
+
### Deprecated
+- climada.hazard.centroids.centr.Centroids.from_lat_lon
+- climada.hazard.centroids.centr.Centroids.def set_area_pixel
+- climada.hazard.centroids.centr.Centroids.def set_area_approx
+- climada.hazard.centroids.centr.Centroids.set_dist_coast
+- climada.hazard.centroids.centr.Centroids.empty_geometry_points
+- climada.hazard.centroids.centr.Centroids.set_meta_to_lat_lon
+- climada.hazard.centroids.centr.Centroids.set_lat_lon_to_meta
+
### Removed
+- climada.hazard.base.Hazard.clear
+- climada.hazard.base.Hazard.raster_to_vector
+- climada.hazard.base.Hazard.read_mat
+- climada.hazard.base.Hazard.reproject_raster
+- climada.hazard.base.Hazard.set_vector
+- climada.hazard.base.Hazard.vector_to_raster
+- climada.hazard.centroids.centr.Centroids.calc_pixels_polygons
+- climada.hazard.centroids.centr.Centroids.check
+- climada.hazard.centroids.centr.Centroids.clear
+- climada.hazard.centroids.centr.Centroids.equal
+- climada.hazard.centroids.centr.Centroids.from_mat
+- climada.hazard.centroids.centr.Centroids.from_base_grid
+- climada.hazard.centroids.centr.Centroids.read_excel
+- climada.hazard.centroids.centr.Centroids.read_hdf5
+- climada.hazard.centroids.centr.Centroids.read_mat
+- climada.hazard.centroids.centr.Centroids.set_elevation
+- climada.hazard.centroids.centr.Centroids.set_geometry_points
+- climada.hazard.centroids.centr.Centroids.set_lat_lon
+- climada.hazard.centroids.centr.Centroids.set_raster_file
+- climada.hazard.centroids.centr.Centroids.set_raster_from_pnt_bounds
+- climada.hazard.centroids.centr.Centroids.set_vector_file
+- climada.hazard.centroids.centr.Centroids.values_from_raster_files
+- climada.hazard.centroids.centr.Centroids.values_from_vector_files
+- climada.hazard.centroids.centr.generate_nat_earth_centroids
+
## 4.1.1
Release date: 2024-02-21
diff --git a/climada/engine/test/test_cost_benefit.py b/climada/engine/test/test_cost_benefit.py
index 19dc11859d..095716fc9c 100644
--- a/climada/engine/test/test_cost_benefit.py
+++ b/climada/engine/test/test_cost_benefit.py
@@ -35,16 +35,15 @@
from climada.test import get_test_file
-HAZ_TEST_MAT = get_test_file('atl_prob_no_name', file_format='matlab')
ENT_TEST_MAT = get_test_file('demo_today', file_format='MAT-file')
-
+HAZ_TEST_TC :Path = get_test_file('test_tc_florida')
class TestSteps(unittest.TestCase):
"""Test intermediate steps"""
def test_calc_impact_measures_pass(self):
"""Test _calc_impact_measures against reference value"""
- self.assertTrue(HAZ_TEST_MAT.is_file(), "{} is not a file".format(HAZ_TEST_MAT))
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ self.assertTrue(HAZ_TEST_TC.is_file(), "{} is not a file".format(HAZ_TEST_TC))
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
self.assertTrue(ENT_TEST_MAT.is_file(), "{} is not a file".format(ENT_TEST_MAT))
entity = Entity.from_mat(ENT_TEST_MAT)
@@ -230,7 +229,7 @@ def test_cb_one_meas_fut_pass(self):
def test_calc_cb_no_change_pass(self):
"""Test _calc_cost_benefit without present value against reference value"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_mat(ENT_TEST_MAT)
entity.measures._data['TC'] = entity.measures._data.pop('XX')
for meas in entity.measures.get_measure('TC'):
@@ -267,7 +266,7 @@ def test_calc_cb_no_change_pass(self):
def test_calc_cb_change_pass(self):
"""Test _calc_cost_benefit with present value against reference value"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_mat(ENT_TEST_MAT)
entity.measures._data['TC'] = entity.measures._data.pop('XX')
for meas in entity.measures.get_measure('TC'):
@@ -438,7 +437,7 @@ def test_norm_value(self):
def test_combine_fut_pass(self):
"""Test combine_measures with present and future"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -498,7 +497,7 @@ def test_combine_fut_pass(self):
def test_combine_current_pass(self):
"""Test combine_measures with only future"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -538,7 +537,7 @@ def test_combine_current_pass(self):
def test_apply_transf_current_pass(self):
"""Test apply_risk_transfer with only future"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -588,7 +587,7 @@ def test_apply_transf_current_pass(self):
def test_apply_transf_cost_fact_pass(self):
"""Test apply_risk_transfer with only future annd cost factor"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -636,7 +635,7 @@ def test_apply_transf_cost_fact_pass(self):
def test_apply_transf_future_pass(self):
"""Test apply_risk_transfer with present and future"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -692,7 +691,7 @@ def test_apply_transf_future_pass(self):
def test_remove_measure(self):
"""Test remove_measure method"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -720,7 +719,7 @@ class TestCalc(unittest.TestCase):
def test_calc_change_pass(self):
"""Test calc with future change"""
# present
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.exposures.gdf.rename(columns={'impf_': 'impf_TC'}, inplace=True)
entity.check()
@@ -777,7 +776,7 @@ def test_calc_change_pass(self):
def test_calc_no_change_pass(self):
"""Test calc without future change"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_excel(ENT_DEMO_TODAY)
entity.check()
entity.exposures.ref_year = 2018
@@ -808,7 +807,7 @@ class TestRiskFuncs(unittest.TestCase):
def test_impact(self):
ent = Entity.from_excel(ENT_DEMO_TODAY)
ent.check()
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
impact = ImpactCalc(ent.exposures, ent.impact_funcs, hazard).impact()
return impact
diff --git a/climada/engine/test/test_forecast.py b/climada/engine/test/test_forecast.py
index b207bc4645..ef249ae203 100644
--- a/climada/engine/test/test_forecast.py
+++ b/climada/engine/test/test_forecast.py
@@ -110,14 +110,12 @@ def test_Forecast_plot(self):
HAZ_DIR.joinpath('storm_europe_cosmoe_forecast_vmax_testfile.nc'),
run_datetime=dt.datetime(2018,1,1),
event_date=dt.datetime(2018,1,3))
- haz1.centroids.lat += 0.6
- haz1.centroids.lon -= 1.2
+ haz1.centroids.gdf.geometry = haz1.centroids.gdf.geometry.translate(-1.2, 0.6)
haz2 = StormEurope.from_cosmoe_file(
HAZ_DIR.joinpath('storm_europe_cosmoe_forecast_vmax_testfile.nc'),
run_datetime=dt.datetime(2018,1,1),
event_date=dt.datetime(2018,1,3))
- haz2.centroids.lat += 0.6
- haz2.centroids.lon -= 1.2
+ haz2.centroids.gdf.geometry = haz2.centroids.gdf.geometry.translate(-1.2, 0.6)
#exposure
data = {}
data['latitude'] = haz1.centroids.lat
diff --git a/climada/entity/exposures/base.py b/climada/entity/exposures/base.py
index 5d57a90751..53504739e0 100644
--- a/climada/entity/exposures/base.py
+++ b/climada/entity/exposures/base.py
@@ -498,7 +498,6 @@ def from_raster(cls, file_name, band=1, src_crs=None, window=None,
--------
Exposures
"""
- exp = cls()
meta, value = u_coord.read_raster(file_name, [band], src_crs, window,
geometry, dst_crs, transform, width,
height, resampling)
@@ -507,14 +506,16 @@ def from_raster(cls, file_name, band=1, src_crs=None, window=None,
lry = uly + meta['height'] * yres
x_grid, y_grid = np.meshgrid(np.arange(ulx + xres / 2, lrx, xres),
np.arange(uly + yres / 2, lry, yres))
+ return cls(
+ {
+ 'longitude': x_grid.flatten(),
+ 'latitude': y_grid.flatten(),
+ 'value': value.reshape(-1),
+ },
+ meta=meta,
+ crs=meta['crs'],
+ )
- if exp.crs is None:
- exp.set_crs()
- exp.gdf['longitude'] = x_grid.flatten()
- exp.gdf['latitude'] = y_grid.flatten()
- exp.gdf['value'] = value.reshape(-1)
- exp.meta = meta
- return exp
def plot_scatter(self, mask=None, ignore_zero=False, pop_name=True,
buffer=0.0, extend='neither', axis=None, figsize=(9, 13),
diff --git a/climada/entity/exposures/test/test_base.py b/climada/entity/exposures/test/test_base.py
index 79f0f51104..4a1146e26d 100644
--- a/climada/entity/exposures/test/test_base.py
+++ b/climada/entity/exposures/test/test_base.py
@@ -60,7 +60,6 @@ def test_assign_pass(self):
np_rand = np.random.RandomState(123456789)
haz = Hazard.from_raster([HAZ_DEMO_FL], haz_type='FL', window=Window(10, 20, 50, 60))
- haz.raster_to_vector()
ncentroids = haz.centroids.size
exp = Exposures(crs=haz.centroids.crs)
@@ -74,8 +73,8 @@ def test_assign_pass(self):
# make sure that it works for both float32 and float64
for test_dtype in [np.float64, np.float32]:
- haz.centroids.lat = haz.centroids.lat.astype(test_dtype)
- haz.centroids.lon = haz.centroids.lon.astype(test_dtype)
+ haz.centroids.gdf["lat"] = haz.centroids.lat.astype(test_dtype)
+ haz.centroids.gdf["lon"] = haz.centroids.lon.astype(test_dtype)
exp.assign_centroids(haz)
self.assertEqual(exp.gdf.shape[0], len(exp.gdf[INDICATOR_CENTR + 'FL']))
np.testing.assert_array_equal(exp.gdf[INDICATOR_CENTR + 'FL'].values, expected_result)
@@ -131,7 +130,7 @@ def test_assign_raster_pass(self):
'width': 20, 'height': 10,
'transform': rasterio.Affine(1.5, 0.0, -20, 0.0, -1.4, 8)
}
- haz = Hazard('FL', centroids=Centroids(meta=meta))
+ haz = Hazard('FL', centroids=Centroids.from_meta(meta))
# explicit points with known results (see `expected_result` for details)
exp = Exposures(crs=DEF_CRS)
@@ -151,9 +150,9 @@ def test_assign_raster_pass(self):
expected_result = [
# constant y-value, varying x-value
- -1, 0, 0, 0, 0, 1,
+ 0, 0, 0, 0, 0, 1,
# constant x-value, varying y-value
- -1, 0, 0, 20,
+ 0, 0, 0, 20,
# out of bounds: topleft, top, topright, right, bottomright, bottom, bottomleft, left
-1, -1, -1, -1, -1, -1, -1, -1,
# some explicit points within the raster
@@ -171,6 +170,9 @@ def test_assign_raster_same_pass(self):
np.testing.assert_array_equal(exp.gdf[INDICATOR_CENTR + 'FL'].values,
np.arange(haz.centroids.size, dtype=int))
+ # Test fails because exposures stores the crs in the meta attribute as rasterio object,
+ # while the centroids stores the crs in the geodataframe, which is not a rasterio object.
+ # The comparison in assign_centroids then fails.
def test_assign_large_hazard_subset_pass(self):
"""Test assign_centroids with raster hazard"""
exp = Exposures.from_raster(HAZ_DEMO_FL, window=Window(10, 20, 50, 60))
@@ -178,11 +180,10 @@ def test_assign_large_hazard_subset_pass(self):
exp.gdf.longitude[[0, 1]] = exp.gdf.longitude[[1, 0]]
exp.check()
haz = Hazard.from_raster([HAZ_DEMO_FL], haz_type='FL')
- haz.raster_to_vector()
exp.assign_centroids(haz)
assigned_centroids = haz.centroids.select(sel_cen=exp.gdf[INDICATOR_CENTR + 'FL'].values)
- np.testing.assert_array_equal(assigned_centroids.lat, exp.gdf.latitude)
- np.testing.assert_array_equal(assigned_centroids.lon, exp.gdf.longitude)
+ np.testing.assert_array_equal(np.unique(assigned_centroids.lat), np.unique(exp.gdf.latitude))
+ np.testing.assert_array_equal(np.unique(assigned_centroids.lon), np.unique(exp.gdf.longitude))
def test_affected_total_value(self):
haz_type = "RF"
diff --git a/climada/entity/measures/test/test_base.py b/climada/entity/measures/test/test_base.py
index 2a25ace36f..0890fec26a 100644
--- a/climada/entity/measures/test/test_base.py
+++ b/climada/entity/measures/test/test_base.py
@@ -39,7 +39,11 @@
DATA_DIR = CONFIG.measures.test_data.dir()
-HAZ_TEST_MAT = get_test_file('atl_prob_no_name', file_format='matlab')
+HAZ_TEST_TC :Path = get_test_file('test_tc_florida', file_format='hdf5')
+"""
+Hazard test file from Data API: Hurricanes from 1851 to 2011 over Florida with 100 centroids.
+Fraction is empty. Format: HDF5.
+"""
ENT_TEST_MAT = Path(exposures_test.__file__).parent / 'data' / 'demo_today.mat'
class TestApply(unittest.TestCase):
@@ -78,7 +82,7 @@ def test_cutoff_hazard_pass(self):
meas = MeasureSet.from_mat(ENT_TEST_MAT)
act_1 = meas.get_measure(name='Seawall')[0]
- haz = Hazard.from_mat(HAZ_TEST_MAT)
+ haz = Hazard.from_hdf5(HAZ_TEST_TC)
exp = Exposures.from_mat(ENT_TEST_MAT)
exp.gdf.rename(columns={'impf': 'impf_TC'}, inplace=True)
exp.check()
@@ -112,7 +116,7 @@ def test_cutoff_hazard_region_pass(self):
act_1 = meas.get_measure(name='Seawall')[0]
act_1.exp_region_id = [1]
- haz = Hazard.from_mat(HAZ_TEST_MAT)
+ haz = Hazard.from_hdf5(HAZ_TEST_TC)
exp = Exposures.from_mat(ENT_TEST_MAT)
exp.gdf['region_id'] = np.zeros(exp.gdf.shape[0])
exp.gdf.region_id.values[10:] = 1
@@ -246,7 +250,7 @@ def test_filter_exposures_pass(self):
imp_set = ImpactFuncSet.from_mat(ENT_TEST_MAT)
- haz = Hazard.from_mat(HAZ_TEST_MAT)
+ haz = Hazard.from_hdf5(HAZ_TEST_TC)
exp.assign_centroids(haz)
new_exp = copy.deepcopy(exp)
@@ -329,7 +333,7 @@ def test_filter_exposures_pass(self):
def test_apply_ref_pass(self):
"""Test apply method: apply all measures but insurance"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_mat(ENT_TEST_MAT)
entity.measures._data['TC'] = entity.measures._data.pop('XX')
@@ -365,7 +369,7 @@ def test_apply_ref_pass(self):
def test_calc_impact_pass(self):
"""Test calc_impact method: apply all measures but insurance"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_mat(ENT_TEST_MAT)
entity.exposures.gdf.rename(columns={'impf': 'impf_TC'}, inplace=True)
@@ -399,7 +403,7 @@ def test_calc_impact_pass(self):
def test_calc_impact_transf_pass(self):
"""Test calc_impact method: apply all measures and insurance"""
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
entity = Entity.from_mat(ENT_TEST_MAT)
entity.exposures.gdf.rename(columns={'impf': 'impf_TC'}, inplace=True)
diff --git a/climada/hazard/base.py b/climada/hazard/base.py
index 12b1d0226a..908e6d01f1 100644
--- a/climada/hazard/base.py
+++ b/climada/hazard/base.py
@@ -29,28 +29,26 @@
from typing import Union, Optional, Callable, Dict, Any, List
import warnings
-import geopandas as gpd
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pathos.pools import ProcessPool as Pool
import rasterio
-from rasterio.features import rasterize
-from rasterio.warp import reproject, Resampling, calculate_default_transform
+import rasterio.features
+import rasterio.warp
import sparse as sp
from scipy import sparse
import xarray as xr
+from climada import CONFIG
from climada.hazard.centroids.centr import Centroids
-import climada.util.plot as u_plot
import climada.util.checker as u_check
+import climada.util.constants as u_const
+import climada.util.coordinates as u_coord
import climada.util.dates_times as u_dt
-from climada import CONFIG
import climada.util.hdf5_handler as u_hdf5
-import climada.util.coordinates as u_coord
-from climada.util.constants import ONE_LAT_KM, DEF_CRS, DEF_FREQ_UNIT
-from climada.util.coordinates import NEAREST_NEIGHBOR_THRESHOLD
+import climada.util.plot as u_plot
LOGGER = logging.getLogger(__name__)
@@ -66,8 +64,8 @@
},
'col_centroids': {'sheet_name': 'centroids',
'col_name': {'cen_id': 'centroid_id',
- 'lat': 'latitude',
- 'lon': 'longitude'
+ 'latitude': 'lat',
+ 'longitude': 'lon',
}
}
}
@@ -174,7 +172,7 @@ def __init__(self,
centroids: Optional[Centroids] = None,
event_id: Optional[np.ndarray] = None,
frequency: Optional[np.ndarray] = None,
- frequency_unit: str = DEF_FREQ_UNIT,
+ frequency_unit: str = u_const.DEF_FREQ_UNIT,
event_name: Optional[List[str]] = None,
date: Optional[np.ndarray] = None,
orig: Optional[np.ndarray] = None,
@@ -227,7 +225,8 @@ def __init__(self,
"""
self.haz_type = haz_type
self.units = units
- self.centroids = centroids if centroids is not None else Centroids()
+ self.centroids = centroids if centroids is not None else Centroids(
+ lat=np.empty(0), lon=np.empty(0))
# following values are defined for each event
self.event_id = event_id if event_id is not None else np.array([], int)
self.frequency = frequency if frequency is not None else np.array(
@@ -260,19 +259,9 @@ def get_default(cls, attribute):
Any
"""
return {
- 'frequency_unit': DEF_FREQ_UNIT,
+ 'frequency_unit': u_const.DEF_FREQ_UNIT,
}.get(attribute)
- def clear(self):
- """Reinitialize attributes (except the process Pool)."""
- for (var_name, var_val) in self.__dict__.items():
- if isinstance(var_val, np.ndarray) and var_val.ndim == 1:
- setattr(self, var_name, np.array([], dtype=var_val.dtype))
- elif isinstance(var_val, sparse.csr_matrix):
- setattr(self, var_name, sparse.csr_matrix(np.empty((0, 0))))
- elif not isinstance(var_val, Pool):
- setattr(self, var_name, self.get_default(var_name) or var_val.__class__())
-
def check(self):
"""Check dimension of attributes.
@@ -280,14 +269,13 @@ def check(self):
------
ValueError
"""
- self.centroids.check()
self._check_events()
@classmethod
def from_raster(cls, files_intensity, files_fraction=None, attrs=None,
band=None, haz_type=None, pool=None, src_crs=None, window=None,
geometry=None, dst_crs=None, transform=None, width=None,
- height=None, resampling=Resampling.nearest):
+ height=None, resampling=rasterio.warp.Resampling.nearest):
"""Create Hazard with intensity and fraction values from raster files
If raster files are masked, the masked values are set to 0.
@@ -323,7 +311,7 @@ def from_raster(cls, files_intensity, files_fraction=None, attrs=None,
reproject to given crs
transform : rasterio.Affine
affine transformation to apply
- wdith : float, optional
+ width : float, optional
number of lons for transform
height : float, optional
number of lats for transform
@@ -351,14 +339,17 @@ def from_raster(cls, files_intensity, files_fraction=None, attrs=None,
if haz_type is not None:
hazard_kwargs["haz_type"] = haz_type
- centroids = Centroids.from_raster_file(
- files_intensity[0], src_crs=src_crs, window=window, geometry=geometry, dst_crs=dst_crs,
- transform=transform, width=width, height=height, resampling=resampling)
+ centroids, meta = Centroids.from_raster_file(
+ files_intensity[0], src_crs=src_crs, window=window,
+ geometry=geometry, dst_crs=dst_crs, transform=transform,
+ width=width, height=height, resampling=resampling, return_meta=True,
+ )
+
if pool:
chunksize = max(min(len(files_intensity) // pool.ncpus, 1000), 1)
inten_list = pool.map(
- centroids.values_from_raster_files,
- [[f] for f in files_intensity],
+ _values_from_raster_files,
+ [[f] for f in files_intensity], itertools.repeat(meta),
itertools.repeat(band), itertools.repeat(src_crs),
itertools.repeat(window), itertools.repeat(geometry),
itertools.repeat(dst_crs), itertools.repeat(transform),
@@ -367,8 +358,8 @@ def from_raster(cls, files_intensity, files_fraction=None, attrs=None,
intensity = sparse.vstack(inten_list, format='csr')
if files_fraction is not None:
fract_list = pool.map(
- centroids.values_from_raster_files,
- [[f] for f in files_fraction],
+ _values_from_raster_files,
+ [[f] for f in files_fraction], itertools.repeat(meta),
itertools.repeat(band), itertools.repeat(src_crs),
itertools.repeat(window), itertools.repeat(geometry),
itertools.repeat(dst_crs), itertools.repeat(transform),
@@ -376,15 +367,16 @@ def from_raster(cls, files_intensity, files_fraction=None, attrs=None,
itertools.repeat(resampling), chunksize=chunksize)
fraction = sparse.vstack(fract_list, format='csr')
else:
- intensity = centroids.values_from_raster_files(
- files_intensity, band=band, src_crs=src_crs, window=window, geometry=geometry,
- dst_crs=dst_crs, transform=transform, width=width, height=height,
- resampling=resampling)
+ intensity = _values_from_raster_files(
+ files_intensity, meta=meta, band=band, src_crs=src_crs, window=window,
+ geometry=geometry, dst_crs=dst_crs, transform=transform, width=width,
+ height=height, resampling=resampling,
+ )
if files_fraction is not None:
- fraction = centroids.values_from_raster_files(
- files_fraction, band=band, src_crs=src_crs, window=window, geometry=geometry,
- dst_crs=dst_crs, transform=transform, width=width, height=height,
- resampling=resampling)
+ fraction = _values_from_raster_files(
+ files_fraction, meta=meta, band=band, src_crs=src_crs, window=window,
+ geometry=geometry, dst_crs=dst_crs, transform=transform, width=width,
+ height=height, resampling=resampling)
if files_fraction is None:
fraction = intensity.copy()
@@ -399,12 +391,6 @@ def set_raster(self, *args, **kwargs):
"Use Hazard.from_raster instead.")
self.__dict__ = Hazard.from_raster(*args, **kwargs).__dict__
- def set_vector(self, *args, **kwargs):
- """This function is deprecated, use Hazard.from_vector."""
- LOGGER.warning("The use of Hazard.set_vector is deprecated."
- "Use Hazard.from_vector instead.")
- self.__dict__ = Hazard.from_vector(*args, **kwargs).__dict__
-
@classmethod
def from_xarray_raster_file(
cls, filepath: Union[pathlib.Path, str], *args, **kwargs
@@ -458,7 +444,7 @@ def from_xarray_raster(
intensity: str = "intensity",
coordinate_vars: Optional[Dict[str, str]] = None,
data_vars: Optional[Dict[str, str]] = None,
- crs: str = DEF_CRS,
+ crs: str = u_const.DEF_CRS,
rechunk: bool = False,
):
"""Read raster-like data from an xarray Dataset
@@ -764,8 +750,10 @@ def from_xarray_raster(
)
# Transform coordinates into centroids
- centroids = Centroids.from_lat_lon(
- data[coords["latitude"]].values, data[coords["longitude"]].values, crs=crs,
+ centroids = Centroids(
+ lat=data[coords["latitude"]].values,
+ lon=data[coords["longitude"]].values,
+ crs=crs,
)
def to_csr_matrix(array: xr.DataArray) -> sparse.csr_matrix:
@@ -1017,72 +1005,6 @@ def vshape(array):
LOGGER.debug("Hazard successfully loaded. Number of events: %i", num_events)
return cls(centroids=centroids, intensity=intensity_matrix, **hazard_kwargs)
- @classmethod
- def from_vector(cls, files_intensity, files_fraction=None, attrs=None,
- inten_name=None, frac_name=None, dst_crs=None, haz_type=None):
- """Read vector files format supported by fiona. Each intensity name is
- considered an event.
-
- Parameters
- ----------
- files_intensity : list(str)
- file names containing intensity, default: ['intensity']
- files_fraction : (list(str))
- file names containing fraction,
- default: ['fraction']
- attrs : dict, optional
- name of Hazard attributes and their values
- inten_name : list(str), optional
- name of variables containing the intensities of each event
- frac_name : list(str), optional
- name of variables containing
- the fractions of each event
- dst_crs : crs, optional
- reproject to given crs
- haz_type : str, optional
- acronym of the hazard type (e.g. 'TC').
- default: None, which will use the class default ('' for vanilla
- `Hazard` objects, hard coded in some subclasses)
-
- Returns
- -------
- haz : climada.hazard.Hazard
- Hazard from vector file
- """
- if not attrs:
- attrs = {}
- if not inten_name:
- inten_name = ['intensity']
- if not frac_name:
- inten_name = ['fraction']
- if files_fraction is not None and len(files_intensity) != len(files_fraction):
- raise ValueError('Number of intensity files differs from fraction files:'
- f' {len(files_intensity)} != {len(files_fraction)}')
-
- hazard_kwargs = {}
- if haz_type is not None:
- hazard_kwargs["haz_type"] = haz_type
-
- if len(files_intensity) > 0:
- centroids = Centroids.from_vector_file(files_intensity[0], dst_crs=dst_crs)
- elif files_fraction is not None and len(files_fraction) > 0:
- centroids = Centroids.from_vector_file(files_fraction[0], dst_crs=dst_crs)
- else:
- centroids = Centroids()
-
- intensity = centroids.values_from_vector_files(
- files_intensity, val_names=inten_name, dst_crs=dst_crs)
- if files_fraction is None:
- fraction = intensity.copy()
- fraction.data.fill(1)
- else:
- fraction = centroids.values_from_vector_files(
- files_fraction, val_names=frac_name, dst_crs=dst_crs)
-
- hazard_kwargs.update(cls._attrs_to_kwargs(attrs, num_events=intensity.shape[0]))
- return cls(
- centroids=centroids, intensity=intensity, fraction=fraction, **hazard_kwargs)
-
@staticmethod
def _attrs_to_kwargs(attrs: Dict[str, Any], num_events: int) -> Dict[str, Any]:
"""Transform attributes to init kwargs or use default values
@@ -1133,169 +1055,17 @@ def _attrs_to_kwargs(attrs: Dict[str, Any], num_events: int) -> Dict[str, Any]:
return kwargs
- def reproject_raster(self, dst_crs=False, transform=None, width=None, height=None,
- resampl_inten=Resampling.nearest, resampl_fract=Resampling.nearest):
- """Change current raster data to other CRS and/or transformation
-
- Parameters
- ----------
- dst_crs: crs, optional
- reproject to given crs
- transform: rasterio.Affine
- affine transformation to apply
- wdith: float
- number of lons for transform
- height: float
- number of lats for transform
- resampl_inten: rasterio.warp,.Resampling optional
- resampling function used for reprojection to dst_crs for intensity
- resampl_fract: rasterio.warp,.Resampling optional
- resampling function used for reprojection to dst_crs for fraction
- """
- if not self.centroids.meta:
- raise ValueError('Raster not set')
- if not dst_crs:
- dst_crs = self.centroids.meta['crs']
- if transform and not width or transform and not height:
- raise ValueError('Provide width and height to given transformation.')
- if not transform:
- transform, width, height = calculate_default_transform(
- self.centroids.meta['crs'], dst_crs, self.centroids.meta['width'],
- self.centroids.meta['height'], self.centroids.meta['transform'][2],
- (self.centroids.meta['transform'][5]
- + self.centroids.meta['height'] * self.centroids.meta['transform'][4]),
- (self.centroids.meta['transform'][2]
- + self.centroids.meta['width'] * self.centroids.meta['transform'][0]),
- self.centroids.meta['transform'][5])
- dst_meta = self.centroids.meta.copy()
- dst_meta.update({'crs': dst_crs, 'transform': transform,
- 'width': width, 'height': height
- })
- intensity = np.zeros((self.size, dst_meta['height'], dst_meta['width']))
- fraction = np.zeros((self.size, dst_meta['height'], dst_meta['width']))
- kwargs = {'src_transform': self.centroids.meta['transform'],
- 'src_crs': self.centroids.meta['crs'],
- 'dst_transform': transform, 'dst_crs': dst_crs,
- 'resampling': resampl_inten}
- for idx_ev, inten in enumerate(self.intensity.toarray()):
- reproject(
- source=np.asarray(inten.reshape((self.centroids.meta['height'],
- self.centroids.meta['width']))),
- destination=intensity[idx_ev, :, :],
- **kwargs)
- kwargs.update(resampling=resampl_fract)
- for idx_ev, fract in enumerate(self.fraction.toarray()):
- reproject(
- source=np.asarray(
- fract.reshape((self.centroids.meta['height'],
- self.centroids.meta['width']))),
- destination=fraction[idx_ev, :, :],
- **kwargs)
- self.centroids.meta = dst_meta
- self.intensity = sparse.csr_matrix(
- intensity.reshape(self.size, dst_meta['height'] * dst_meta['width']))
- self.fraction = sparse.csr_matrix(
- fraction.reshape(self.size, dst_meta['height'] * dst_meta['width']))
- self.check()
-
- def reproject_vector(self, dst_crs, scheduler=None):
+ def reproject_vector(self, dst_crs):
"""Change current point data to a a given projection
Parameters
----------
dst_crs : crs
reproject to given crs
- scheduler : str, optional
- used for dask map_partitions. “threads”,
- “synchronous” or “processes”
- """
- self.centroids.set_geometry_points(scheduler)
- self.centroids.geometry = self.centroids.geometry.to_crs(dst_crs)
- self.centroids.lat = self.centroids.geometry[:].y
- self.centroids.lon = self.centroids.geometry[:].x
- self.check()
-
- def raster_to_vector(self):
- """Change current raster to points (center of the pixels)"""
- self.centroids.set_meta_to_lat_lon()
- self.centroids.meta = dict()
- self.check()
-
- def vector_to_raster(self, scheduler=None):
- """Change current point data to a raster with same resolution
-
- Parameters
- ----------
- scheduler : str, optional
- used for dask map_partitions. “threads”,
- “synchronous” or “processes”
"""
- points_df = gpd.GeoDataFrame()
- points_df['latitude'] = self.centroids.lat
- points_df['longitude'] = self.centroids.lon
- val_names = ['val' + str(i_ev) for i_ev in range(2 * self.size)]
- for i_ev, inten_name in enumerate(val_names):
- if i_ev < self.size:
- points_df[inten_name] = np.asarray(self.intensity[i_ev, :].toarray()).reshape(-1)
- else:
- points_df[inten_name] = np.asarray(self.fraction[i_ev - self.size, :].toarray()). \
- reshape(-1)
- raster, meta = u_coord.points_to_raster(points_df, val_names,
- crs=self.centroids.geometry.crs,
- scheduler=scheduler)
- self.intensity = sparse.csr_matrix(raster[:self.size, :, :].reshape(self.size, -1))
- self.fraction = sparse.csr_matrix(raster[self.size:, :, :].reshape(self.size, -1))
- self.centroids = Centroids(meta=meta)
+ self.centroids.gdf.to_crs(dst_crs, inplace=True)
self.check()
- def read_mat(self, *args, **kwargs):
- """This function is deprecated, use Hazard.from_mat."""
- LOGGER.warning("The use of Hazard.read_mat is deprecated."
- "Use Hazard.from_mat instead.")
- self.__dict__ = Hazard.from_mat(*args, **kwargs).__dict__
-
- @classmethod
- def from_mat(cls, file_name, var_names=None):
- """Read climada hazard generate with the MATLAB code in .mat format.
-
- Parameters
- ----------
- file_name : str
- absolute file name
- var_names : dict, optional
- name of the variables in the file,
- default: DEF_VAR_MAT constant
-
- Returns
- -------
- haz : climada.hazard.Hazard
- Hazard object from the provided MATLAB file
-
- Raises
- ------
- KeyError
- """
- # pylint: disable=protected-access
- if not var_names:
- var_names = DEF_VAR_MAT
- LOGGER.info('Reading %s', file_name)
- try:
- data = u_hdf5.read(file_name)
- try:
- data = data[var_names['field_name']]
- except KeyError:
- pass
-
- centroids = Centroids.from_mat(file_name, var_names=var_names['var_cent'])
- attrs = cls._read_att_mat(data, file_name, var_names, centroids)
- haz = cls(haz_type=u_hdf5.get_string(data[var_names['var_name']['per_id']]),
- centroids=centroids,
- **attrs
- )
- except KeyError as var_err:
- raise KeyError("Variable not in MAT file: " + str(var_err)) from var_err
- return haz
-
def read_excel(self, *args, **kwargs):
"""This function is deprecated, use Hazard.from_excel."""
LOGGER.warning("The use of Hazard.read_excel is deprecated."
@@ -1334,7 +1104,8 @@ def from_excel(cls, file_name, var_names=None, haz_type=None):
if haz_type is not None:
hazard_kwargs["haz_type"] = haz_type
try:
- centroids = Centroids.from_excel(file_name, var_names=var_names['col_centroids'])
+ centroids = Centroids._legacy_from_excel(
+ file_name, var_names=var_names['col_centroids'])
hazard_kwargs.update(cls._read_att_excel(file_name, var_names, centroids))
except KeyError as var_err:
raise KeyError("Variable not in Excel file: " + str(var_err)) from var_err
@@ -1426,7 +1197,6 @@ def select(self, event_names=None, event_id=None, date=None, orig=None,
LOGGER.info('No hazard centroids within extent and region')
return None
- sel_cen = sel_cen.nonzero()[0]
for (var_name, var_val) in self.__dict__.items():
if isinstance(var_val, np.ndarray) and var_val.ndim == 1 \
and var_val.size > 0:
@@ -1450,7 +1220,7 @@ def select(self, event_names=None, event_id=None, date=None, orig=None,
LOGGER.warning("Resetting the frequency is based on the calendar year of given"
" dates but the frequency unit here is %s. Consider setting the frequency"
" manually for the selection or changing the frequency unit to %s.",
- self.frequency_unit, DEF_FREQ_UNIT)
+ self.frequency_unit, u_const.DEF_FREQ_UNIT)
year_span_old = np.abs(dt.datetime.fromordinal(self.date.max()).year -
dt.datetime.fromordinal(self.date.min()).year) + 1
year_span_new = np.abs(dt.datetime.fromordinal(haz.date.max()).year -
@@ -1460,7 +1230,7 @@ def select(self, event_names=None, event_id=None, date=None, orig=None,
haz.sanitize_event_ids()
return haz
- def select_tight(self, buffer=NEAREST_NEIGHBOR_THRESHOLD/ONE_LAT_KM,
+ def select_tight(self, buffer=u_coord.NEAREST_NEIGHBOR_THRESHOLD / u_const.ONE_LAT_KM,
val='intensity'):
"""
Reduce hazard to those centroids spanning a minimal box which
@@ -1566,7 +1336,6 @@ def plot_rp_intensity(self, return_periods=(25, 50, 100, 250),
axis, inten_stats: matplotlib.axes._subplots.AxesSubplot, np.ndarray
intenstats is return_periods.size x num_centroids
"""
- self._set_coords_centroids()
inten_stats = self.local_exceedance_inten(np.array(return_periods))
colbar_name = 'Intensity (' + self.units + ')'
title = list()
@@ -1612,7 +1381,6 @@ def plot_intensity(self, event=None, centr=None, smooth=True, axis=None, adapt_f
------
ValueError
"""
- self._set_coords_centroids()
col_label = f'Intensity ({self.units})'
crs_epsg, _ = u_plot.get_transformation(self.centroids.geometry.crs)
if event is not None:
@@ -1662,7 +1430,6 @@ def plot_fraction(self, event=None, centr=None, smooth=True, axis=None,
------
ValueError
"""
- self._set_coords_centroids()
col_label = 'Fraction'
if event is not None:
if isinstance(event, str):
@@ -1792,7 +1559,7 @@ def set_frequency(self, yearrange=None):
if self.frequency_unit not in ['1/year', 'annual', '1/y', '1/a']:
LOGGER.warning("setting the frequency on a hazard object who's frequency unit"
"is %s and not %s will most likely lead to unexpected results",
- self.frequency_unit, DEF_FREQ_UNIT)
+ self.frequency_unit, u_const.DEF_FREQ_UNIT)
if not yearrange:
delta_time = dt.datetime.fromordinal(int(np.max(self.date))).year - \
dt.datetime.fromordinal(int(np.min(self.date))).year + 1
@@ -1810,35 +1577,68 @@ def size(self):
"""Return number of events."""
return self.event_id.size
- def write_raster(self, file_name, intensity=True):
- """Write intensity or fraction as GeoTIFF file. Each band is an event
+ def write_raster(self, file_name, variable='intensity', output_resolution=None):
+ """Write intensity or fraction as GeoTIFF file. Each band is an event.
+ Output raster is always a regular grid (same resolution for lat/lon).
+
+ Note that if output_resolution is not None, the data is rasterized to that
+ resolution. This is an expensive operation. For hazards that are already
+ a raster, output_resolution=None saves on this raster which is efficient.
+
+ If you want to save both fraction and intensity, create two separate files.
+ These can then be read together with the from_raster method.
Parameters
----------
file_name: str
file name to write in tif format
- intensity: bool
- if True, write intensity, otherwise write fraction
+ variable: str
+ if 'intensity', write intensity, if 'fraction' write fraction.
+ Default is 'intensity'
+ output_resolution: int
+ If not None, the data is rasterized to this resolution.
+ Default is None (resolution is estimated from the data).
+
+ See also
+ --------
+ from_raster:
+ method to read intensity and fraction raster files.
"""
- variable = self.intensity
- if not intensity:
- variable = self.fraction
- if self.centroids.meta:
- u_coord.write_raster(file_name, variable.toarray(), self.centroids.meta)
+
+ if variable == 'intensity':
+ var_to_write = self.intensity
+ elif variable =='fraction':
+ var_to_write = self.fraction
else:
- pixel_geom = self.centroids.calc_pixels_polygons()
- profile = self.centroids.meta
- profile.update(driver='GTiff', dtype=rasterio.float32, count=self.size)
- with rasterio.open(file_name, 'w', **profile) as dst:
+ raise ValueError(
+ f"The variable {variable} is not valid. Please use 'intensity' or 'fraction'."
+ )
+
+ meta = self.centroids.get_meta(resolution=output_resolution)
+ meta.update(driver='GTiff', dtype=rasterio.float32, count=self.size)
+ 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(file_name, var_to_write.toarray(), meta)
+ else:
+ geometry = self.centroids.get_pixel_shapes(res=res)
+ with rasterio.open(file_name, 'w', **meta) as dst:
LOGGER.info('Writing %s', file_name)
- for i_ev in range(variable.shape[0]):
- raster = rasterize(
- [(x, val) for (x, val) in
- zip(pixel_geom, np.array(variable[i_ev, :].toarray()).reshape(-1))],
- out_shape=(profile['height'], profile['width']),
- transform=profile['transform'], fill=0,
- all_touched=True, dtype=profile['dtype'], )
- dst.write(raster.astype(profile['dtype']), i_ev + 1)
+ for i_ev in range(self.size):
+ raster = rasterio.features.rasterize(
+ (
+ (geom, value)
+ for geom, value
+ in zip(geometry, var_to_write[i_ev].toarray().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)
def write_hdf5(self, file_name, todense=False):
"""Write hazard in hdf5 format.
@@ -1856,7 +1656,9 @@ def write_hdf5(self, file_name, todense=False):
str_dt = h5py.special_dtype(vlen=str)
for (var_name, var_val) in self.__dict__.items():
if var_name == 'centroids':
- self.centroids.write_hdf5(hf_data.create_group(var_name))
+ # Centroids have their own write_hdf5 method,
+ # which is invoked at the end of this method (s.b.)
+ pass
elif isinstance(var_val, sparse.csr_matrix):
if todense:
hf_data.create_dataset(var_name, data=var_val.toarray())
@@ -1884,6 +1686,7 @@ def write_hdf5(self, file_name, todense=False):
"%s being set to its default value.",
var_name, var_val.__class__.__name__, var_name
)
+ self.centroids.write_hdf5(file_name, mode='a')
def read_hdf5(self, *args, **kwargs):
"""This function is deprecated, use Hazard.from_hdf5."""
@@ -1916,9 +1719,8 @@ def from_hdf5(cls, file_name):
if var_name not in hf_data.keys():
continue
if var_name == 'centroids':
- hazard_kwargs["centroids"] = Centroids.from_hdf5(
- hf_data.get(var_name))
- elif isinstance(var_val, np.ndarray) and var_val.ndim == 1:
+ continue
+ if isinstance(var_val, np.ndarray) and var_val.ndim == 1:
hazard_kwargs[var_name] = np.array(hf_data.get(var_name))
elif isinstance(var_val, sparse.csr_matrix):
hf_csr = hf_data.get(var_name)
@@ -1936,15 +1738,10 @@ def from_hdf5(cls, file_name):
u_hdf5.to_string, np.array(hf_data.get(var_name)).tolist())]
else:
hazard_kwargs[var_name] = hf_data.get(var_name)
-
+ hazard_kwargs["centroids"] = Centroids.from_hdf5(file_name)
# Now create the actual object we want to return!
return cls(**hazard_kwargs)
- def _set_coords_centroids(self):
- """If centroids are raster, set lat and lon coordinates"""
- if self.centroids.meta and not self.centroids.coord.size:
- self.centroids.set_meta_to_lat_lon()
-
def _events_set(self):
"""Generate set of tuples with (event_name, event_date)"""
ev_set = set()
@@ -2042,15 +1839,20 @@ def _centr_plot(self, centr_idx, mat_var, col_name, axis=None, **kwargs):
except IndexError as err:
raise ValueError(f'Wrong centroid id: {centr_idx}.') from err
array_val = mat_var[:, centr_pos].toarray()
- title = f'Centroid {centr_idx}: ({coord[centr_pos, 0]}, {coord[centr_pos, 1]})'
+ title = (
+ f'Centroid {centr_idx}:'
+ f' ({np.around(coord[centr_pos, 0], 3)}, {np.around(coord[centr_pos, 1],3)})'
+ )
elif centr_idx < 0:
max_inten = np.asarray(np.sum(mat_var, axis=0)).reshape(-1)
centr_pos = np.argpartition(max_inten, centr_idx)[centr_idx:]
centr_pos = centr_pos[np.argsort(max_inten[centr_pos])][0]
array_val = mat_var[:, centr_pos].toarray()
- title = (f'{np.abs(centr_idx)}-largest Centroid. {centr_pos}:'
- f' ({coord[centr_pos, 0]}, {coord[centr_pos, 1]})')
+ title = (
+ f'{np.abs(centr_idx)}-largest Centroid. {centr_pos}:'
+ f' ({np.around(coord[centr_pos, 0], 3)}, {np.around(coord[centr_pos, 1], 3)})'
+ )
else:
array_val = np.max(mat_var, axis=1).toarray()
title = f'{self.haz_type} max intensity at each event'
@@ -2398,7 +2200,7 @@ def concat(cls, haz_list):
haz_concat.append(*haz_list)
return haz_concat
- def change_centroids(self, centroids, threshold=NEAREST_NEIGHBOR_THRESHOLD):
+ def change_centroids(self, centroids, threshold=u_coord.NEAREST_NEIGHBOR_THRESHOLD):
"""
Assign (new) centroids to hazard.
@@ -2438,31 +2240,23 @@ def change_centroids(self, centroids, threshold=NEAREST_NEIGHBOR_THRESHOLD):
haz_new_cent = copy.deepcopy(self)
haz_new_cent.centroids = centroids
- # indices for mapping matrices onto common centroids
- if centroids.meta:
- new_cent_idx = u_coord.match_grid_points(
- self.centroids.lon, self.centroids.lat,
- centroids.meta['width'], centroids.meta['height'],
- centroids.meta['transform'])
- if -1 in new_cent_idx:
- raise ValueError("At least one hazard centroid is out of"
- "the raster defined by centroids.meta."
- " Please choose a larger raster.")
- else:
- new_cent_idx = u_coord.match_coordinates(
- self.centroids.coord, centroids.coord, threshold=threshold
+
+ new_cent_idx = u_coord.match_coordinates(
+ self.centroids.coord, centroids.coord, threshold=threshold
+ )
+ if -1 in new_cent_idx:
+ raise ValueError(
+ "At least one hazard centroid is at a larger distance than the given threshold"
+ f" {threshold} from the given centroids. Please choose a larger threshold or"
+ " enlarge the centroids"
)
- if -1 in new_cent_idx:
- raise ValueError("At least one hazard centroid is at a larger "
- f"distance than the given threshold {threshold} "
- "from the given centroids. Please choose a "
- "larger threshold or enlarge the centroids")
if np.unique(new_cent_idx).size < new_cent_idx.size:
- raise ValueError("At least two hazard centroids are mapped to the same "
- "centroids. Please make sure that the given centroids "
- "cover the same area like the original centroids and "
- "are not of lower resolution.")
+ raise ValueError(
+ "At least two hazard centroids are mapped to the same centroids. Please make sure"
+ " that the given centroids cover the same area like the original centroids and are"
+ " not of lower resolution."
+ )
# re-assign attributes intensity and fraction
for attr_name in ["intensity", "fraction"]:
@@ -2579,3 +2373,70 @@ def _get_fraction(self, cent_idx=None):
if cent_idx is None:
return self.fraction
return self.fraction[:, cent_idx]
+
+
+def _values_from_raster_files(
+ file_names, meta, band=None, src_crs=None, window=None,
+ geometry=None, dst_crs=None, transform=None, width=None,
+ height=None, resampling=rasterio.warp.Resampling.nearest,
+):
+ """Read raster of bands and set 0 values to the masked ones.
+
+ Each band is an event. Select region using window or geometry. Reproject input by proving
+ dst_crs and/or (transform, width, height).
+
+ The main purpose of this function is to read intensity/fraction values from raster files for
+ use in Hazard.read_raster. It is implemented as a separate helper function (instead of a
+ class method) to allow for parallel computing.
+
+ Parameters
+ ----------
+ file_names : str
+ path of the file
+ meta : dict
+ description of the centroids raster
+ band : list(int), optional
+ band number to read. Default: [1]
+ src_crs : crs, optional
+ source CRS. Provide it if error without it.
+ window : rasterio.windows.Window, optional
+ window to read
+ geometry : list of shapely.geometry, optional
+ consider pixels only within these shapes
+ dst_crs : crs, optional
+ reproject to given crs
+ transform : rasterio.Affine
+ affine transformation to apply
+ wdith : float
+ number of lons for transform
+ height : float
+ number of lats for transform
+ resampling : rasterio.warp,.Resampling optional
+ resampling function used for reprojection to dst_crs
+
+ Raises
+ ------
+ ValueError
+
+ Returns
+ -------
+ inten : scipy.sparse.csr_matrix
+ Each row is an event.
+ """
+ if band is None:
+ band = [1]
+
+ values = []
+ for file_name in file_names:
+ tmp_meta, data = u_coord.read_raster(
+ file_name, band, src_crs, window, geometry, dst_crs,
+ transform, width, height, resampling,
+ )
+ if (tmp_meta['crs'] != meta['crs']
+ or tmp_meta['transform'] != meta['transform']
+ or tmp_meta['height'] != meta['height']
+ or tmp_meta['width'] != meta['width']):
+ raise ValueError('Raster data is inconsistent with contained raster.')
+ values.append(sparse.csr_matrix(data))
+
+ return sparse.vstack(values, format='csr')
diff --git a/climada/hazard/centroids/centr.py b/climada/hazard/centroids/centr.py
index ae124d151b..89ec5f0a34 100644
--- a/climada/hazard/centroids/centr.py
+++ b/climada/hazard/centroids/centr.py
@@ -20,1186 +20,891 @@
"""
import copy
+from deprecation import deprecated
import logging
from pathlib import Path
-from typing import Optional, Dict, Any
+from typing import Any, Literal, Union
+import warnings
+import h5py
import cartopy.crs as ccrs
import geopandas as gpd
-import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
-from pyproj.crs import CRS
+from pyproj.crs.crs import CRS
import rasterio
-from rasterio.warp import Resampling
-from scipy import sparse
from shapely.geometry.point import Point
-from climada.util.constants import (DEF_CRS,
- ONE_LAT_KM,
- NATEARTH_CENTROIDS)
+from climada.util.constants import DEF_CRS
import climada.util.coordinates as u_coord
-import climada.util.hdf5_handler as u_hdf5
import climada.util.plot as u_plot
__all__ = ['Centroids']
PROJ_CEA = CRS.from_user_input({'proj': 'cea'})
-DEF_VAR_MAT = {
- 'field_names': ['centroids', 'hazard'],
- 'var_name': {
- 'lat': 'lat',
- 'lon': 'lon',
- 'dist_coast': 'distance2coast_km',
- 'admin0_name': 'admin0_name',
- 'admin0_iso3': 'admin0_ISO3',
- 'comment': 'comment',
- 'region_id': 'NatId'
- }
-}
-"""MATLAB variable names"""
-
-DEF_VAR_EXCEL = {
- 'sheet_name': 'centroids',
- 'col_name': {
- 'region_id': 'region_id',
- 'lat': 'latitude',
- 'lon': 'longitude',
- }
-}
-"""Excel variable names"""
-
LOGGER = logging.getLogger(__name__)
+DEF_SHEET_NAME = 'centroids'
+EXP_SPECIFIC_COLS = [
+ 'value',
+ 'impf_',
+ 'centr_',
+ 'cover',
+ 'deductible'
+]
class Centroids():
- """Contains raster or vector centroids.
+ """Contains vector centroids as a GeoDataFrame
Attributes
----------
- meta : dict, optional
- rasterio meta dictionary containing raster properties: width, height, crs and transform
- must be present at least. The affine ransformation needs to be shearless (only stretching)
- and have positive x- and negative y-orientation.
- lat : np.array, optional
- latitude of size size
- lon : np.array, optional
- longitude of size size
- geometry : gpd.GeoSeries, optional
- contains lat and lon crs. Might contain geometry points for lat and lon
- area_pixel : np.array, optional
- area of size size
- dist_coast : np.array, optional
- distance to coast of size size
- on_land : np.array, optional
- on land (True) and on sea (False) of size size
+ lat : np.array
+ Latitudinal coordinates in the specified CRS (can be any unit).
+ lon : np.array
+ Longitudinal coordinates in the specified CRS (can be any unit).
+ crs : pyproj.CRS
+ Coordinate reference system. Default: EPSG:4326 (WGS84)
region_id : np.array, optional
- country region code of size size
- elevation : np.array, optional
- elevation of size size
+ Numeric country (or region) codes. Default: None
+ on_land : np.array, optional
+ Boolean array indicating on land (True) or off shore (False). Default: None
"""
- vars_check = {'lat', 'lon', 'geometry', 'area_pixel', 'dist_coast',
- 'on_land', 'region_id', 'elevation'}
- """Variables whose size will be checked"""
-
def __init__(
self,
- lat: Optional[np.ndarray] = None,
- lon: Optional[np.ndarray] = None,
- geometry: Optional[gpd.GeoSeries] = None,
- meta: Optional[Dict[Any, Any]] = None,
- area_pixel: Optional[np.ndarray] = None,
- on_land: Optional[np.ndarray] = None,
- region_id: Optional[np.ndarray] = None,
- elevation: Optional[np.ndarray] = None,
- dist_coast: Optional[np.ndarray] = None
+ *,
+ lat: Union[np.ndarray, list[float]],
+ lon: Union[np.ndarray, list[float]],
+ crs: Any = DEF_CRS,
+ region_id: Union[Literal["country"], None, np.ndarray, list[float]] = None,
+ on_land: Union[Literal["natural_earth"], None, np.ndarray, list[bool]] = None,
+ **kwargs,
):
"""Initialization
Parameters
----------
- lat : np.array, optional
- latitude of size size. Defaults to empty array
- lon : np.array, optional
- longitude of size size. Defaults to empty array
- geometry : gpd.GeoSeries, optional
- contains lat and lon crs. Might contain geometry points for lat and lon.
- Defaults to empty gpd.Geoseries with crs=DEF_CRS
- meta : dict, optional
- rasterio meta dictionary containing raster properties: width, height, crs and
- transform must be present at least. The affine ransformation needs to be
- shearless (only stretching) and have positive x- and negative y-orientation.
- Defaults to empty dict()
- area_pixel : np.array, optional
- area of size size. Defaults to empty array
- on_land : np.array, optional
- on land (True) and on sea (False) of size size. Defaults to empty array
- region_id : np.array, optional
- country region code of size size, Defaults to empty array
- elevation : np.array, optional
- elevation of size size. Defaults to empty array
- dist_coast : np.array, optional
- distance to coast of size size. Defaults to empty array
+ lat : np.array
+ Latitudinal coordinates in the specified CRS (can be any unit).
+ lon : np.array
+ Longitudinal coordinates in the specified CRS (can be any unit).
+ crs : str or anything accepted by pyproj.CRS.from_user_input()
+ Coordinate reference system. Default: EPSG:4326 (WGS84)
+ region_id : np.array or str, optional
+ Array of numeric country (or region) codes. If the special value "country" is given
+ admin-0 codes are automatically assigned. Default: None
+ on_land : np.array or str, optional
+ Boolean array indicating on land (True) or off shore (False). If the special value
+ "natural_earth" is given, the property is automatically determined from NaturalEarth
+ shapes. Default: None
+ kwargs : dict
+ Additional columns with data to store in the internal GeoDataFrame (gdf attribute).
"""
- self.lat = lat if lat is not None else np.array([])
- self.lon = lon if lon is not None else np.array([])
- self.geometry = geometry if geometry is not None else gpd.GeoSeries(crs=DEF_CRS)
- self.meta = meta if meta is not None else dict()
- self.area_pixel = area_pixel if area_pixel is not None else np.array([])
- self.on_land = on_land if on_land is not None else np.array([])
- self.region_id = region_id if region_id is not None else np.array([])
- self.elevation = elevation if elevation is not None else np.array([])
- self.dist_coast = dist_coast if dist_coast is not None else np.array([])
-
- def check(self):
- """Check integrity of stored information.
-
- Checks that either `meta` attribute is set, or `lat`, `lon` and `geometry.crs`.
- Checks sizes of (optional) data attributes."""
- n_centr = self.size
- for var_name, var_val in self.__dict__.items():
- if var_name in self.vars_check:
- if var_val.size > 0 and var_val.size != n_centr:
- raise ValueError(f'Wrong {var_name} size: {n_centr} != {var_val.size}.')
- if self.meta:
- for name in ['width', 'height', 'crs', 'transform']:
- if name not in self.meta.keys():
- raise ValueError('Missing meta information: %s' % name)
- xres, xshear, _xoff, yshear, yres, _yoff = self.meta['transform'][:6]
- if xshear != 0 or yshear != 0:
- raise ValueError('Affine transformations with shearing components are not '
- 'supported.')
- if yres > 0 or xres < 0:
- raise ValueError('Affine transformations with positive y-orientation '
- 'or negative x-orientation are not supported.')
-
- def equal(self, centr):
- """Return True if two centroids equal, False otherwise
+ self.gdf = gpd.GeoDataFrame(
+ data={
+ 'geometry': gpd.points_from_xy(lon, lat, crs=crs),
+ 'region_id': region_id,
+ 'on_land': on_land,
+ **kwargs,
+ }
+ )
- Parameters
- ----------
- centr : Centroids
- centroids to compare
+ if isinstance(region_id, str):
+ LOGGER.info('Setting region id to %s level.', region_id)
+ self.set_region_id(level=region_id, overwrite=True)
+ if isinstance(on_land, str):
+ LOGGER.info('Setting on land from %s source.', on_land)
+ self.set_on_land(source=on_land, overwrite=True)
- Returns
- -------
- eq : bool
- """
- if self.meta and centr.meta:
- return (u_coord.equal_crs(self.meta['crs'], centr.meta['crs'])
- and self.meta['height'] == centr.meta['height']
- and self.meta['width'] == centr.meta['width']
- and self.meta['transform'] == centr.meta['transform'])
- return (u_coord.equal_crs(self.crs, centr.crs)
- and self.lat.shape == centr.lat.shape
- and self.lon.shape == centr.lon.shape
- and np.allclose(self.lat, centr.lat)
- and np.allclose(self.lon, centr.lon))
+ @property
+ def lat(self):
+ """ Return latitudes """
+ return self.gdf.geometry.y.values
- @staticmethod
- def from_base_grid(land=False, res_as=360, base_file=None):
- """Initialize from base grid data provided with CLIMADA
+ @property
+ def lon(self):
+ """ Return longitudes """
+ return self.gdf.geometry.x.values
- Parameters
- ----------
- land : bool, optional
- If True, restrict to grid points on land. Default: False.
- res_as : int, optional
- Base grid resolution in arc-seconds (one of 150, 360). Default: 360.
- base_file : str, optional
- If set, read this file instead of one provided with climada.
- """
- if base_file is None:
- base_file = NATEARTH_CENTROIDS[res_as]
-
- centroids = Centroids.from_hdf5(base_file)
- if centroids.meta:
- xres, xshear, xoff, yshear, yres, yoff = centroids.meta['transform'][:6]
- shape = (centroids.meta['height'], centroids.meta['width'])
- if yres > 0:
- # make sure y-orientation is negative
- centroids.meta['transform'] = rasterio.Affine(xres, xshear, xoff, yshear,
- -yres, yoff + (shape[0] - 1) * yres)
- # flip y-axis in data arrays
- for name in ["region_id", "dist_coast"]:
- if not hasattr(centroids, name):
- continue
- data = getattr(centroids, name)
- if data.size == 0:
- continue
- setattr(centroids, name, np.flipud(data.reshape(shape)).reshape(-1))
- if land:
- land_reg_ids = list(range(1, 1000))
- land_reg_ids.remove(10) # Antarctica
- centroids = centroids.select(reg_id=land_reg_ids)
-
- centroids.check()
- return centroids
+ @property
+ def geometry(self):
+ """ Return the geometry """
+ return self.gdf['geometry']
- @classmethod
- def from_geodataframe(cls, gdf, geometry_alias='geom'):
- """Create Centroids instance from GeoDataFrame.
+ @property
+ def on_land(self):
+ """ Get the on_land property """
+ if self.gdf.on_land.isna().all():
+ return None
+ return self.gdf['on_land'].values
+
+ @property
+ def region_id(self):
+ """ Get the assigned region_id """
+ if self.gdf.region_id.isna().all():
+ return None
+ return self.gdf['region_id'].values
+
+ @property
+ def crs(self):
+ """ Get the crs"""
+ return self.gdf.crs
- .. deprecated:: 3.3
- This method will be removed in a future version. Pass the data you want to
- construct the Centroids with to the constructor instead.
+ @property
+ def size(self):
+ """Get size (number of lat/lon pairs)"""
+ return self.gdf.shape[0]
- The geometry, lat, and lon attributes are set from the GeoDataFrame.geometry attribute,
- while the columns are copied as attributes to the Centroids object in the form of
- numpy.ndarrays using pandas.Series.to_numpy. The Series dtype will thus be respected.
+ @property
+ def shape(self):
+ """Get shape [lat, lon] assuming rastered data."""
+ return (np.unique(self.lat).size, np.unique(self.lon).size)
- Columns named lat or lon are ignored, as they would overwrite the coordinates extracted
- from the point features. If the geometry attribute bears an alias, it can be dropped by
- setting the geometry_alias parameter.
+ @property
+ def total_bounds(self):
+ """Get total bounds (minx, miny, maxx, maxy)."""
+ return self.gdf.total_bounds
- If the GDF includes a region_id column, but no on_land column, then on_land=True is
- inferred for those centroids that have a set region_id.
+ @property
+ def coord(self):
+ """Get [lat, lon] array."""
+ return np.stack([self.lat, self.lon], axis=1)
- Example
- -------
- >>> gdf = geopandas.read_file('centroids.shp')
- >>> gdf.region_id = gdf.region_id.astype(int) # type coercion
- >>> centroids = Centroids.from_geodataframe(gdf)
+ def __eq__(self, other):
+ """ dunder method for Centroids comparison.
+ returns True if two centroids equal, False otherwise
Parameters
----------
- gdf : GeoDataFrame
- Where the geometry column needs to consist of point features. See above for details on
- processing.
- geometry_alias : str, opt
- Alternate name for the geometry column; dropped to avoid duplicate assignment.
+ other : Centroids
+ object to compare with
Returns
-------
- centr : Centroids
- Centroids with data from given GeoDataFrame
+ bool
"""
- LOGGER.warning(
- "Centroids.from_geodataframe has been deprecated and will be removed in a "
- "future version. Use ther default constructor instead."
- )
+ if not u_coord.equal_crs(self.crs, other.crs):
+ return False
- geometry = gdf.geometry
- lat = gdf.geometry.y.to_numpy(copy=True)
- lon = gdf.geometry.x.to_numpy(copy=True)
- centroids = cls(lat=lat, lon=lon, geometry=geometry)
+ try:
+ pd.testing.assert_frame_equal(self.gdf, other.gdf, check_like=True)
+ return True
+ except AssertionError:
+ return False
- for col in gdf.columns:
- if col in [geometry_alias, 'geometry', 'lat', 'lon']:
- continue # skip these, because they're already set above
- val = gdf[col].to_numpy(copy=True)
- setattr(centroids, col, val)
+ def to_default_crs(self, inplace=True):
+ """Project the current centroids to the default CRS (epsg4326)
- if centroids.on_land.size == 0:
- try:
- centroids.on_land = ~np.isnan(centroids.region_id)
- except KeyError:
- pass
+ Parameters
+ ----------
+ inplace: bool
+ if True, modifies the centroids in place.
+ if False, return projected centroids object.
+ Default is True.
- return centroids
+ Returns
+ -------
+ Centroids or None (if inplace is True)
- @classmethod
- def from_pix_bounds(cls, xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon, crs=DEF_CRS):
- """Create Centroids object with meta attribute according to pixel border data.
+ """
+ return self.to_crs(DEF_CRS, inplace=inplace)
- .. deprecated:: 3.3
- This method will be removed in a future version. CLIMADA will only support
- regular grids with a constant lat/lon resolution then. Use
- :py:meth:`from_pnt_bounds` instead.
+ def to_crs(self, crs, inplace=False):
+ """ Project the current centroids to the desired crs
Parameters
----------
- xf_lat : float
- upper latitude (top)
- xo_lon : float
- left longitude
- d_lat : float
- latitude step (negative)
- d_lon : float
- longitude step (positive)
- n_lat : int
- number of latitude points
- n_lon : int
- number of longitude points
- crs : dict() or rasterio.crs.CRS, optional
- CRS. Default: DEF_CRS
+ crs : str
+ coordinate reference system
+ inplace: bool, default False
+ if True, modifies the centroids in place.
+ if False, returns a copy.
Returns
-------
- centr : Centroids
- Centroids with meta according to given pixel border data.
+ Centroids or None (if inplace is True)
"""
- LOGGER.warning(
- "Centroids.from_pix_bounds has been deprecated and will be removed in a "
- "future version. Use Centroids.from_pnt_bounds instead."
- )
-
- meta = {
- 'dtype': 'float32',
- 'width': n_lon,
- 'height': n_lat,
- 'crs': crs,
- 'transform': rasterio.Affine(d_lon, 0.0, xo_lon, 0.0, d_lat, xf_lat),
- }
-
- return cls(meta=meta)
-
- def set_raster_from_pnt_bounds(self, *args, **kwargs):
- """This function is deprecated, use Centroids.from_pnt_bounds instead."""
- LOGGER.warning("The use of Centroids.set_raster_from_pnt_bounds is deprecated. "
- "Use Centroids.from_pnt_bounds instead.")
- self.__dict__ = Centroids.from_pnt_bounds(*args, **kwargs).__dict__
+ if inplace:
+ self.gdf.to_crs(crs, inplace=True)
+ return None
+ return Centroids.from_geodataframe(self.gdf.to_crs(crs))
@classmethod
- def from_pnt_bounds(cls, points_bounds, res, crs=DEF_CRS):
- """Create Centroids object with meta attribute according to points border data.
-
- raster border = point border + res/2
+ def from_geodataframe(cls, gdf):
+ """Initialize centroids from a geodataframe
Parameters
----------
- points_bounds : tuple
- points' lon_min, lat_min, lon_max, lat_max
- res : float
- desired resolution in same units as points_bounds
- crs : dict() or rasterio.crs.CRS, optional
- CRS. Default: DEF_CRS
+ gdf : GeoDataFrame
+ Input geodataframe with centroids as points
+ in the geometry column. All other columns are
+ attached to the centroids geodataframe.
Returns
-------
- centr : Centroids
- Centroids with meta according to given points border data.
+ Centroids
+ Centroids built from the geodataframe.
+
+ Raises
+ ------
+ ValueError
"""
- rows, cols, ras_trans = u_coord.pts_to_raster_meta(points_bounds, (res, -res))
- meta = {
- 'width': cols,
- 'height': rows,
- 'crs': crs,
- 'transform': ras_trans,
- }
- return cls(meta=meta)
+ if (gdf.geom_type != 'Point').any():
+ raise ValueError(
+ 'The inpute geodataframe contains geometries that are not points.'
+ )
- def set_lat_lon(self, *args, **kwargs):
- """This function is deprecated, use Centroids.from_lat_lon instead."""
- LOGGER.warning("The use of Centroids.set_lat_lon is deprecated. "
- "Use Centroids.from_lat_lon instead.")
- self.__dict__ = Centroids.from_lat_lon(*args, **kwargs).__dict__
+ # Don't forget to make a copy!!
+ # This is a bit ugly, but avoids to recompute the geometries
+ # in the init. For large datasets this saves computation time
+ centroids = cls(lat=[1], lon=[1]) #make "empty" centroids
+ centroids.gdf = gdf.copy(deep=True)
+ if gdf.crs is None:
+ centroids.gdf.set_crs(DEF_CRS, inplace=True)
+ return centroids
@classmethod
- def from_lat_lon(cls, lat, lon, crs=DEF_CRS):
- """Create Centroids object from given latitude, longitude and CRS.
+ def from_exposures(cls, exposures):
+ """Generate centroids from the locations of exposures.
+
+ The properties "region_id" and "on_land" are also extracted from the Exposures object if
+ available. The columns "value", "impf_*", "centr_*", "cover", and "deductible" are not
+ used.
Parameters
----------
- lat : np.array
- latitude
- lon : np.array
- longitude
- crs : dict() or rasterio.crs.CRS, optional
- CRS. Default: DEF_CRS
+ exposures : Exposures
+ Exposures from which to take the centroids locations (as well as region_id and on_land
+ if available).
Returns
-------
- centr : Centroids
- Centroids with points according to given coordinates
+ Centroids
+
+ Raises
+ ------
+ ValueError
"""
- lat = np.asarray(lat)
- lon = np.asarray(lon)
- geometry = gpd.GeoSeries(crs=crs)
- return cls(lat=lat, lon=lon, geometry=geometry)
-
- def set_raster_file(self, file_name, band=None, **kwargs):
- """This function is deprecated, use Centroids.from_raster_file
- and Centroids.values_from_raster_files instead."""
- LOGGER.warning("The use of Centroids.set_raster_file is deprecated. "
- "Use Centroids.from_raster_file and "
- "Centroids.values_from_raster_files instead.")
- if not self.meta:
- self.__dict__ = Centroids.from_raster_file(file_name, **kwargs).__dict__
- return self.values_from_raster_files([file_name], band=band, **kwargs)
+ col_names = [
+ column for column in exposures.gdf.columns
+ if not any(pattern in column for pattern in EXP_SPECIFIC_COLS)
+ ]
+
+ # Legacy behaviour
+ # Exposures can be without geometry column
+ #TODO: remove once exposures is real geodataframe with geometry.
+ if 'geometry' in exposures.gdf.columns:
+ gdf = exposures.gdf[col_names]
+ return cls.from_geodataframe(gdf)
+
+ if 'latitude' in exposures.gdf.columns and 'longitude' in exposures.gdf.columns:
+ gdf = exposures.gdf[col_names]
+ return cls(
+ lat=exposures.gdf['latitude'],
+ lon=exposures.gdf['longitude'],
+ crs=exposures.crs,
+ **dict(gdf.items()),
+ )
+
+ raise ValueError(
+ "The given exposures object has no coordinates information."
+ "The exposures' GeoDataFrame must have either point geometries"
+ " or latitude and longitude values."
+ )
@classmethod
- def from_raster_file(cls, file_name, src_crs=None, window=None,
- geometry=None, dst_crs=None, transform=None, width=None,
- height=None, resampling=Resampling.nearest):
- """Create a new Centroids object from a raster file
+ def from_pnt_bounds(cls, points_bounds, res, crs=DEF_CRS):
+ """Create Centroids object from coordinate bounds and resolution.
- Select region using window or geometry. Reproject input by providing
- dst_crs and/or (transform, width, height).
+ The result contains all points from a regular raster with the given resolution and CRS,
+ covering the given bounds. Note that the raster bounds are larger than the points' bounds
+ by res/2.
Parameters
----------
- file_name : str
- path of the file
- src_crs : crs, optional
- source CRS. Provide it if error without it.
- window : rasterio.windows.Window, optional
- window to read
- geometry : list of shapely.geometry, optional
- consider pixels only within these shapes
- dst_crs : crs, optional
- reproject to given crs
- transform : rasterio.Affine
- affine transformation to apply
- wdith : float
- number of lons for transform
- height : float
- number of lats for transform
- resampling : rasterio.warp,.Resampling optional
- resampling function used for reprojection to dst_crs
+ points_bounds : tuple
+ The bounds (lon_min, lat_min, lon_max, lat_max) of the point coordinates.
+ res : float
+ The desired resolution in same units as `points_bounds`.
+ crs : dict() or rasterio.crs.CRS, optional
+ Coordinate reference system. Default: DEF_CRS
Returns
-------
- centr : Centroids
- Centroids with meta attribute according to the given raster file
+ Centroids
"""
- meta, _ = u_coord.read_raster(
- file_name, [1], src_crs, window, geometry, dst_crs,
- transform, width, height, resampling)
- return cls(meta=meta)
+ height, width, transform = u_coord.pts_to_raster_meta(points_bounds, (res, -res))
+ return cls.from_meta({
+ "crs": crs,
+ "width": width,
+ "height": height,
+ "transform": transform,
+ })
- def values_from_raster_files(self, file_names, band=None, src_crs=None, window=None,
- geometry=None, dst_crs=None, transform=None, width=None,
- height=None, resampling=Resampling.nearest):
- """Read raster of bands and set 0 values to the masked ones.
+ def append(self, centr):
+ """Append Centroids
- Each band is an event. Select region using window or geometry. Reproject input by proving
- dst_crs and/or (transform, width, height).
+ Note that the result might contain duplicate points if the object to append has an overlap
+ with the current object.
Parameters
----------
- file_names : str
- path of the file
- band : list(int), optional
- band number to read. Default: [1]
- src_crs : crs, optional
- source CRS. Provide it if error without it.
- window : rasterio.windows.Window, optional
- window to read
- geometry : list of shapely.geometry, optional
- consider pixels only within these shapes
- dst_crs : crs, optional
- reproject to given crs
- transform : rasterio.Affine
- affine transformation to apply
- wdith : float
- number of lons for transform
- height : float
- number of lats for transform
- resampling : rasterio.warp,.Resampling optional
- resampling function used for reprojection to dst_crs
+ centr : Centroids
+ Centroids to append. The centroids need to have the same CRS.
Raises
------
ValueError
- Returns
- -------
- inten : scipy.sparse.csr_matrix
- Each row is an event.
+ See Also
+ --------
+ union : Union of Centroid objects.
+ remove_duplicate_points : Remove duplicate points in a Centroids object.
"""
- if band is None:
- band = [1]
-
- values = []
- for file_name in file_names:
- tmp_meta, data = u_coord.read_raster(
- file_name, band, src_crs, window, geometry, dst_crs,
- transform, width, height, resampling)
- if (tmp_meta['crs'] != self.meta['crs']
- or tmp_meta['transform'] != self.meta['transform']
- or tmp_meta['height'] != self.meta['height']
- or tmp_meta['width'] != self.meta['width']):
- raise ValueError('Raster data is inconsistent with contained raster.')
- values.append(sparse.csr_matrix(data))
-
- return sparse.vstack(values, format='csr')
-
-
- def set_vector_file(self, file_name, inten_name=None, **kwargs):
- """This function is deprecated, use Centroids.from_vector_file
- and Centroids.values_from_vector_files instead."""
- LOGGER.warning("The use of Centroids.set_vector_file is deprecated. "
- "Use Centroids.from_vector_file and "
- "Centroids.values_from_vector_files instead.")
- if not self.geometry.any():
- self.__dict__ = Centroids.from_vector_file(file_name, **kwargs).__dict__
- return self.values_from_vector_files([file_name], val_names=inten_name, **kwargs)
+ if not u_coord.equal_crs(self.crs, centr.crs):
+ raise ValueError(
+ "The centroids have different Coordinate-Reference-Systems (CRS)"
+ )
+ self.gdf = pd.concat([self.gdf, centr.gdf])
- @classmethod
- def from_vector_file(cls, file_name, dst_crs=None):
- """Create Centroids object from vector file (any format supported by fiona).
+ def union(self, *others):
+ """Create the union of Centroids objects
+
+ All centroids must have the same CRS. Points that are contained in more than one of the
+ Centroids objects will only be contained once (i.e. duplicates are removed).
Parameters
----------
- file_name : str
- vector file with format supported by fiona and 'geometry' field.
- dst_crs : crs, optional
- reproject to given crs
+ others : list of Centroids
+ Centroids contributing to the union.
Returns
-------
- centr : Centroids
- Centroids with points according to the given vector file
+ centroids : Centroids
+ Centroids object containing the union of all Centroids.
"""
- lat, lon, geometry, _ = u_coord.read_vector(
- file_name, [], dst_crs=dst_crs)
- return cls(lat=lat, lon=lon, geometry=geometry)
-
- def values_from_vector_files(self, file_names, val_names=None, dst_crs=None):
- """Read intensity or other data from vector files, making sure that geometry is compatible.
+ centroids = copy.deepcopy(self)
+ for cent in others:
+ centroids.append(cent)
+ return centroids.remove_duplicate_points()
- If the geometry of the shapes in any of the given files does not agree with the
- geometry of this Centroids instance, a ValueError is raised.
+ def remove_duplicate_points(self):
+ """Return a copy of centroids with duplicate points removed
Parameters
----------
- file_names : list(str)
- vector files with format supported by fiona and 'geometry' field.
- val_names : list(str), optional
- list of names of the columns of the values. Default: ['intensity']
- dst_crs : crs, optional
- reproject to given crs
-
- Raises
- ------
- ValueError
+ centr : Centroids
+ Centroids with or without duplicate points
Returns
-------
- values : scipy.sparse.csr_matrix
- Sparse array of shape (len(val_name), len(geometry)).
+ centroids : Centroids
+ A new Centroids object that contains a subselection of the original centroids without
+ duplicates. Note that a copy is returned even if there were no duplicates.
"""
- if val_names is None:
- val_names = ["intensity"]
+ return self.from_geodataframe(self.gdf.drop_duplicates(subset=["geometry"]))
- values = []
- for file_name in file_names:
- tmp_lat, tmp_lon, tmp_geometry, data = u_coord.read_vector(
- file_name, val_names, dst_crs=dst_crs
- )
- try:
- assert u_coord.equal_crs(tmp_geometry.crs, self.geometry.crs)
- np.testing.assert_allclose(tmp_lat, self.lat)
- np.testing.assert_allclose(tmp_lon, self.lon)
- except AssertionError as exc:
- raise ValueError(
- "Vector data inconsistent with contained vector"
- ) from exc
- values.append(sparse.csr_matrix(data))
-
- return sparse.vstack(values, format="csr")
-
- def read_mat(self, *args, **kwargs):
- """This function is deprecated, use Centroids.from_mat instead."""
- LOGGER.warning("The use of Centroids.read_mat is deprecated."
- "Use Centroids.from_mat instead.")
- self.__dict__ = Centroids.from_mat(*args, **kwargs).__dict__
+ def select(self, reg_id=None, extent=None, sel_cen=None):
+ """Return new Centroids object containing points following certain criteria
- @classmethod
- def from_mat(cls, file_name, var_names=None):
- """Read centroids from CLIMADA's MATLAB version.
+ It is currently possible to filter by region (reg_id), by geographical extent (extent), or
+ by an explicit list of indices/a mask (sel_cen). If more than one criterion is given, all
+ of them must be satisfied for a point to be included in the selection.
Parameters
----------
- file_name : str
- absolute or relative file name
- var_names : dict, optional
- name of the variables
-
- Raises
- ------
- KeyError
+ reg_id : int or list of int, optional
+ Numeric ID (or IDs) of the region (or regions) to restrict to, according to the values
+ in the region_id property. Default: None
+ extent : tuple, optional
+ The geographical extent (min_lon, max_lon, min_lat, max_lat) to restrict to, including
+ the boundary. If the value for min_lon is greater than max_lon, the extent is
+ interpreted to cross the antimeridian ([max_lon, 180] and [-180, min_lon]).
+ Default: None
+ sel_cen : np.ndarray of int or bool, optional
+ Boolean mask, or list of indices to restrict to. Default: None
Returns
-------
- centr : Centroids
- Centroids with data from the given file
+ centroids : Centroids
+ Sub-selection of this object
"""
- LOGGER.info('Reading %s', file_name)
- if var_names is None:
- var_names = DEF_VAR_MAT
-
- cent = u_hdf5.read(file_name)
- # Try open encapsulating variable FIELD_NAMES
- num_try = 0
- for field in var_names['field_names']:
- try:
- cent = cent[field]
- break
- except KeyError:
- num_try += 1
- if num_try == len(var_names['field_names']):
- LOGGER.warning("Variables are not under: %s.", var_names['field_names'])
+ sel_cen_bool = sel_cen
+ if sel_cen is not None and sel_cen.dtype.kind == 'i':
+ # if needed, convert indices to bool
+ sel_cen_bool = np.zeros(self.size, dtype=bool)
+ sel_cen_bool[np.unique(sel_cen)] = True
- try:
- cen_lat = np.squeeze(cent[var_names['var_name']['lat']])
- cen_lon = np.squeeze(cent[var_names['var_name']['lon']])
- centr = cls.from_lat_lon(cen_lat, cen_lon)
-
- try:
- centr.dist_coast = np.squeeze(cent[var_names['var_name']['dist_coast']])
- except KeyError:
- pass
- try:
- centr.region_id = np.squeeze(cent[var_names['var_name']['region_id']])
- except KeyError:
- pass
- except KeyError as err:
- raise KeyError("Not existing variable: %s" % str(err)) from err
-
- return centr
+ sel_cen_mask = self.select_mask(sel_cen=sel_cen_bool, reg_id=reg_id, extent=extent)
+ return Centroids.from_geodataframe(self.gdf.iloc[sel_cen_mask])
- def read_excel(self, *args, **kwargs):
- """This function is deprecated, use Centroids.from_excel instead."""
- LOGGER.warning("The use of Centroids.read_excel is deprecated."
- "Use Centroids.from_excel instead.")
- self.__dict__ = Centroids.from_excel(*args, **kwargs).__dict__
- @classmethod
- def from_excel(cls, file_name, var_names=None):
- """Generate a new centroids object from an excel file with column names in var_names.
+ def select_mask(self, sel_cen=None, reg_id=None, extent=None):
+ """Create mask of selected centroids
Parameters
----------
- file_name : str
- absolute or relative file name
- var_names : dict, default
- name of the variables
-
- Raises
- ------
- KeyError
+ sel_cen: np.ndarray of bool, optional
+ Boolean array, with size matching the number of centroids. Default: None
+ reg_id : int or list of int, optional
+ Numeric ID (or IDs) of the region (or regions) to restrict to, according to the values
+ in the region_id property. Default: None
+ extent : tuple, optional
+ The geographical extent (min_lon, max_lon, min_lat, max_lat) to restrict to, including
+ the boundary. If the value for min_lon is greater than lon_max, the extent is
+ interpreted to cross the antimeridian ([lon_max, 180] and [-180, lon_min]).
+ Default: None
Returns
-------
- centr : Centroids
- Centroids with data from the given file
+ sel_cen : np.ndarray of bool
+ Boolean array (mask) with value True for centroids in selection.
"""
- LOGGER.info('Reading %s', file_name)
- if var_names is None:
- var_names = DEF_VAR_EXCEL
+ if sel_cen is None:
+ sel_cen = np.ones(self.size, dtype=bool)
+ if reg_id is not None:
+ sel_cen &= np.isin(self.region_id, reg_id)
+ if extent is not None:
+ lon_min, lon_max, lat_min, lat_max = extent
+ lon_max += 360 if lon_min > lon_max else 0
+ lon_normalized = u_coord.lon_normalize(
+ self.lon.copy(), center=0.5 * (lon_min + lon_max))
+ sel_cen &= (
+ (lon_normalized >= lon_min) & (lon_normalized <= lon_max) &
+ (self.lat >= lat_min) & (self.lat <= lat_max)
+ )
+ return sel_cen
- try:
- dfr = pd.read_excel(file_name, var_names['sheet_name'])
- centr = cls.from_lat_lon(dfr[var_names['col_name']['lat']],
- dfr[var_names['col_name']['lon']])
- try:
- centr.region_id = dfr[var_names['col_name']['region_id']]
- except KeyError:
- pass
+ #TODO replace with nice GeoDataFrame util plot method.
+ def plot(self, axis=None, figsize=(9, 13), **kwargs):
+ """Plot centroids scatter points over earth
- except KeyError as err:
- raise KeyError("Not existing variable: %s" % str(err)) from err
+ Parameters
+ ----------
+ axis : matplotlib.axes._subplots.AxesSubplot, optional
+ axis to use
+ figsize: (float, float), optional
+ figure size for plt.subplots
+ The default is (9, 13)
+ kwargs : optional
+ arguments for scatter matplotlib function
- return centr
+ Returns
+ -------
+ axis : matplotlib.axes._subplots.AxesSubplot
+ """
+ pad = np.abs(u_coord.get_resolution(self.lat, self.lon)).min()
- def clear(self):
- """Clear vector and raster data."""
- self.__init__()
+ proj_data, _ = u_plot.get_transformation(self.crs)
+ proj_plot = proj_data
+ if isinstance(proj_data, ccrs.PlateCarree):
+ # use different projections for plot and data to shift the central lon in the plot
+ xmin, ymin, xmax, ymax = u_coord.latlon_bounds(self.lat, self.lon, buffer=pad)
+ proj_plot = ccrs.PlateCarree(central_longitude=0.5 * (xmin + xmax))
+ else:
+ xmin, ymin, xmax, ymax = (self.lon.min() - pad, self.lat.min() - pad,
+ self.lon.max() + pad, self.lat.max() + pad)
- def append(self, centr):
- """Append centroids points.
+ if not axis:
+ _fig, axis, _fontsize = u_plot.make_map(proj=proj_plot, figsize=figsize)
- If centr or self are rasters they are converted to points first using
- Centroids.set_meta_to_lat_lon. Note that self is modified in-place,
- and meta is set to {}. Thus, raster information in self is lost.
+ axis.set_extent((xmin, xmax, ymin, ymax), crs=proj_data)
+ u_plot.add_shapes(axis)
+ axis.scatter(self.lon, self.lat, transform=proj_data, **kwargs)
+ plt.tight_layout()
+ return axis
- Note: this is a wrapper for centroids.union.
+ def set_region_id(self, level='country', overwrite=False):
+ """Set region_id as country ISO numeric code attribute for every pixel or point.
Parameters
----------
- centr : Centroids
- Centroids to append. The centroids need to have the same CRS.
-
- See Also
- --------
- union : Union of Centroid objects.
+ level: str, optional
+ The admin level on which to assign centroids. Currently only 'country' (admin0) is
+ implemented. Default: 'country'
+ overwrite : bool, optional
+ If True, overwrite the existing region_id information. If False, region_id is set
+ only if region_id is missing (None). Default: False
"""
- self.__dict__.update(self.union(centr).__dict__)
+ if overwrite or self.region_id is None:
+ LOGGER.debug('Setting region_id %s points.', str(self.size))
+ if level == 'country':
+ ne_geom = self._ne_crs_geom()
+ self.gdf['region_id'] = u_coord.get_country_code(
+ ne_geom.y.values, ne_geom.x.values,
+ )
+ else:
+ raise NotImplementedError(
+ 'The region id can only be assigned for countries so far'
+ )
+ def set_on_land(self, source='natural_earth', overwrite=False):
+ """Set on_land attribute for every pixel or point.
- def union(self, *others):
+ Parameters
+ ----------
+ source: str, optional
+ The source of the on-land information. Currently, only 'natural_earth' (based on shapes
+ from NaturalEarth, https://www.naturalearthdata.com/) is implemented.
+ Default: 'natural_earth'.
+ overwrite : bool, optional
+ If True, overwrite the existing on_land information. If False, on_land is set
+ only if on_land is missing (None). Default: False
"""
- Create the union of centroids from the inputs.
+ if overwrite or self.on_land is None:
+ LOGGER.debug('Setting on_land %s points.', str(self.lat.size))
+ if source=='natural_earth':
+ ne_geom = self._ne_crs_geom()
+ self.gdf['on_land'] = u_coord.coord_on_land(
+ ne_geom.y.values, ne_geom.x.values
+ )
+ else:
+ raise NotImplementedError(
+ 'The on land variables can only be automatically assigned using natural earth.'
+ )
+
+ def get_pixel_shapes(self, res=None, **kwargs):
+ """Create a GeoSeries of the quadratic pixel shapes at the centroid locations
- The centroids are combined together point by point.
- Rasters are converted to points and raster information is lost
- in the output. All centroids must have the same CRS.
+ Note that this assumes that the centroids define a regular grid of pixels.
- In any case, the attribute .geometry is computed for all centroids.
- This requires a CRS to be defined. If Centroids.crs is None, the
- default DEF_CRS is set for all centroids (self and others).
+ Parameters
+ ----------
+ res : float, optional
+ The resolution of the regular grid the pixels are taken from. If not given, it is
+ estimated using climada.util.coordinates.get_resolution. Default: None
+ kwargs : optional
+ Additional keyword arguments are passed to climada.util.coordinates.get_resolution.
- When at least one centroids has one of the following property
- defined, it is also computed for all others.
- .area_pixel, .dist_coast, .on_land, .region_id, .elevetaion'
+ Returns
+ -------
+ GeoSeries
- !Caution!: the input objects (self and others) are modified in place.
- Missing properties are added, existing ones are not overwritten.
+ See also
+ --------
+ climada.util.coordinates.get_resolution
+ """
+ if res is None:
+ res = np.abs(u_coord.get_resolution(self.lat, self.lon, **kwargs)).min()
+ geom = self.geometry.copy()
+ # unset CRS to avoid warnings about geographic CRS when using `GeoSeries.buffer`
+ geom.crs = None
+ return geom.buffer(
+ # resolution=1, cap_style=3: squared buffers
+ # https://shapely.readthedocs.io/en/latest/manual.html#object.buffer
+ distance=res / 2, resolution=1, cap_style=3,
+ # reset CRS (see above)
+ ).set_crs(self.crs)
+
+ def get_area_pixel(self, min_resol=1.0e-8):
+ """Compute the area per centroid in the CEA projection
+
+ Note that this assumes that the centroids define a regular grid of pixels (area in m²).
Parameters
----------
- others : any number of climada.hazard.Centroids()
- Centroids to form the union with
+ min_resol : float, optional
+ When estimating the grid resolution, use this as the minimum resolution in lat and lon.
+ It is passed to climada.util.coordinates.get_resolution. Default: 1.0e-8
Returns
-------
- centroids : Centroids
- Centroids containing the union of the centroids in others.
-
- Raises
- ------
- ValueError
+ areapixels : np.array
+ Area of each pixel in square meters.
"""
- # restrict to non-empty centroids
- cent_list = [c for c in (self,) + others if c.size > 0 or c.meta] # pylint: disable=no-member
- if len(cent_list) == 0 or len(others) == 0:
- return copy.deepcopy(self)
-
- # check if all centroids are identical
- if all([cent_list[0].equal(cent) for cent in cent_list[1:]]):
- return copy.deepcopy(cent_list[0])
-
- # convert all raster centroids to point centroids
- for cent in cent_list:
- if cent.meta and not cent.lat.any():
- cent.set_meta_to_lat_lon()
-
- # make sure that all Centroids have the same CRS
- for cent in cent_list:
- if cent.crs is None:
- cent.geometry = cent.geometry.set_crs(DEF_CRS)
- if not u_coord.equal_crs(cent.crs, cent_list[0].crs):
- raise ValueError('In a union, all Centroids need to have the same CRS: '
- f'{cent.crs} != {cent_list[0].crs}')
-
- # set attributes that are missing in some but defined in others
- for attr in ["geometry", "area_pixel", "dist_coast", "on_land", "region_id", "elevation"]:
- if np.any([getattr(cent, attr).size > 0 for cent in cent_list]):
- for cent in cent_list:
- if not getattr(cent, attr).size > 0:
- fun_name = f"set_{attr}{'_points' if attr == 'geometry' else ''}"
- getattr(Centroids, fun_name)(cent)
-
- # create new Centroids object and set concatenated attributes
- centroids = Centroids()
- for attr_name, attr_val in vars(cent_list[0]).items():
- if isinstance(attr_val, np.ndarray) and attr_val.ndim == 1:
- attr_val_list = [getattr(cent, attr_name) for cent in cent_list]
- setattr(centroids, attr_name, np.hstack(attr_val_list))
- elif isinstance(attr_val, gpd.GeoSeries):
- attr_val_list = [getattr(cent, attr_name) for cent in cent_list]
- setattr(centroids, attr_name, pd.concat(attr_val_list, ignore_index=True))
-
- # finally, remove duplicate points
- return centroids.remove_duplicate_points()
+ LOGGER.debug('Computing pixel area for %d centroids.', self.size)
+ xy_pixels = self.get_pixel_shapes(min_resol=min_resol)
+ if PROJ_CEA != xy_pixels.crs:
+ xy_pixels = xy_pixels.to_crs(crs={'proj': 'cea'})
+ return xy_pixels.area.values
- def get_closest_point(self, x_lon, y_lat, scheduler=None):
+ def get_closest_point(self, x_lon, y_lat):
"""Returns closest centroid and its index to a given point.
Parameters
----------
x_lon : float
- x coord (lon)
+ Longitudinal (x) coordinate.
y_lat : float
- y coord (lat)
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
+ Latitudinal (y) coordinate.
Returns
-------
x_close : float
x-coordinate (longitude) of closest centroid.
y_close : float
- y-coordinate (latitude) of closest centroids.
+ y-coordinate (latitude) of closest centroid.
idx_close : int
Index of centroid in internal ordering of centroids.
"""
- if self.meta:
- if not self.lat.size or not self.lon.size:
- self.set_meta_to_lat_lon()
- i_lat, i_lon = rasterio.transform.rowcol(self.meta['transform'], x_lon, y_lat)
- i_lat = np.clip(i_lat, 0, self.meta['height'] - 1)
- i_lon = np.clip(i_lon, 0, self.meta['width'] - 1)
- close_idx = int(i_lat * self.meta['width'] + i_lon)
- else:
- self.set_geometry_points(scheduler)
- close_idx = self.geometry.distance(Point(x_lon, y_lat)).values.argmin()
+ close_idx = self.geometry.distance(Point(x_lon, y_lat)).values.argmin()
return self.lon[close_idx], self.lat[close_idx], close_idx
- def set_region_id(self, scheduler=None):
- """Set region_id as country ISO numeric code attribute for every pixel or point.
+ # NOT REALLY AN ELEVATION FUNCTION, JUST READ RASTER
+ def get_elevation(self, topo_path):
+ """Return elevation attribute for every pixel or point in meters.
Parameters
----------
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
+ topo_path : str
+ Path to a raster file containing gridded elevation data.
+
+ Returns
+ -------
+ values : np.array of shape (npoints,)
+ Interpolated elevation values from raster file for each given coordinate point.
"""
- ne_geom = self._ne_crs_geom(scheduler)
- LOGGER.debug('Setting region_id %s points.', str(self.lat.size))
- self.region_id = u_coord.get_country_code(
- ne_geom.geometry[:].y.values, ne_geom.geometry[:].x.values)
+ return u_coord.read_raster_sample(topo_path, self.lat, self.lon)
- def set_area_pixel(self, min_resol=1.0e-8, scheduler=None):
- """Set `area_pixel` attribute for every pixel or point (area in m*m).
+ def get_dist_coast(self, signed=False, precomputed=False):
+ """Get dist_coast attribute for every pixel or point in meters.
Parameters
----------
- min_resol : float, optional
- if centroids are points, use this minimum resolution in lat and lon. Default: 1.0e-8
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
+ signed : bool, optional
+ If True, use signed distances (positive off shore and negative on land). Default: False
+ precomputed : bool, optional
+ If True, use precomputed distances (from NASA). Default: False
+
+ Returns
+ -------
+ dist : np.array
+ (Signed) distance to coast in meters.
"""
- if self.meta:
- if hasattr(self.meta['crs'], 'linear_units') and \
- str.lower(self.meta['crs'].linear_units) in ['m', 'metre', 'meter']:
- self.area_pixel = np.zeros((self.meta['height'], self.meta['width']))
- self.area_pixel *= abs(self.meta['transform'].a) * abs(self.meta['transform'].e)
- return
- if abs(abs(self.meta['transform'].a) -
- abs(self.meta['transform'].e)) > 1.0e-5:
- raise ValueError('Area can not be computed for not squared pixels.')
- res = self.meta['transform'].a
- else:
- res = u_coord.get_resolution(self.lat, self.lon, min_resol=min_resol)
- res = np.abs(res).min()
- self.set_geometry_points(scheduler)
- LOGGER.debug('Setting area_pixel %s points.', str(self.lat.size))
- xy_pixels = self.geometry.buffer(res / 2).envelope
- if PROJ_CEA == self.geometry.crs:
- self.area_pixel = xy_pixels.area.values
- else:
- self.area_pixel = xy_pixels.to_crs(crs={'proj': 'cea'}).area.values
+ ne_geom = self._ne_crs_geom()
+ if precomputed:
+ return u_coord.dist_to_coast_nasa(
+ ne_geom.y.values, ne_geom.x.values, highres=True, signed=signed)
+ LOGGER.debug('Computing distance to coast for %s centroids.', str(self.size))
+ return u_coord.dist_to_coast(ne_geom, signed=signed)
- def set_area_approx(self, min_resol=1.0e-8):
- """Set `area_pixel` attribute for every pixel or point (approximate area in m*m).
+ def get_meta(self, resolution=None):
+ """Returns a meta raster based on the centroids bounds.
- Values are differentiated per latitude. Faster than `set_area_pixel`.
+ Note that this function is not perfectly inverse with `from_meta` since `get_meta` enforces
+ a grid with equal resolution in x- and y-direction with coordinates increasing in
+ x-direction and decreasing in y-direction.
Parameters
----------
- min_resol : float, optional
- if centroids are points, use this minimum resolution in lat and lon. Default: 1.0e-8
+ resolution : float, optional
+ Resolution of the raster. If not given, the resolution is estimated from the centroids
+ by assuming that they form a regular raster. Default: None
+
+ Returns
+ -------
+ meta: dict
+ meta raster representation of the centroids
"""
- if self.meta:
- if hasattr(self.meta['crs'], 'linear_units') and \
- str.lower(self.meta['crs'].linear_units) in ['m', 'metre', 'meter']:
- self.area_pixel = np.zeros((self.meta['height'], self.meta['width']))
- self.area_pixel *= abs(self.meta['transform'].a) * abs(self.meta['transform'].e)
- return
- res_lat, res_lon = self.meta['transform'].e, self.meta['transform'].a
- lat_unique = np.arange(self.meta['transform'].f + res_lat / 2,
- self.meta['transform'].f + self.meta['height'] * res_lat,
- res_lat)
- lon_unique_len = self.meta['width']
- res_lat = abs(res_lat)
- else:
- res_lat, res_lon = np.abs(
- u_coord.get_resolution(self.lat, self.lon, min_resol=min_resol))
- lat_unique = np.array(np.unique(self.lat))
- lon_unique_len = len(np.unique(self.lon))
- if PROJ_CEA == self.geometry.crs:
- self.area_pixel = np.repeat(res_lat * res_lon, lon_unique_len)
- return
-
- LOGGER.debug('Setting area_pixel approx %s points.', str(self.lat.size))
- res_lat = res_lat * ONE_LAT_KM * 1000
- res_lon = res_lon * ONE_LAT_KM * 1000 * np.cos(np.radians(lat_unique))
- area_approx = np.repeat(res_lat * res_lon, lon_unique_len)
- if area_approx.size == self.size:
- self.area_pixel = area_approx
- else:
- raise ValueError('Pixel area of points can not be computed.')
+ if resolution is None:
+ resolution = np.abs(u_coord.get_resolution(self.lat, self.lon)).min()
+ xmin, ymin, xmax, ymax = self.gdf.total_bounds
+ rows, cols, ras_trans = u_coord.pts_to_raster_meta(
+ (xmin, ymin, xmax, ymax),
+ (resolution, -resolution),
+ )
+ meta = {
+ 'crs': self.crs,
+ 'height': rows,
+ 'width': cols,
+ 'transform': ras_trans,
+ }
+ return meta
- def set_elevation(self, topo_path):
- """Set elevation attribute for every pixel or point in meters.
+ ##
+ # I/O methods
+ ##
- Parameters
- ----------
- topo_path : str
- Path to a raster file containing gridded elevation data.
- """
- if not self.coord.size:
- self.set_meta_to_lat_lon()
- self.elevation = u_coord.read_raster_sample(topo_path, self.lat, self.lon)
+ @classmethod
+ def from_raster_file(cls, file_name, src_crs=None, window=None, geometry=None,
+ dst_crs=None, transform=None, width=None, height=None,
+ resampling=rasterio.warp.Resampling.nearest, return_meta=False):
+ """Create a new Centroids object from a raster file
- def set_dist_coast(self, signed=False, precomputed=False, scheduler=None):
- """Set dist_coast attribute for every pixel or point in meters.
+ Select region using window or geometry. Reproject input by providing
+ dst_crs and/or (transform, width, height).
Parameters
----------
- signed : bool
- If True, use signed distances (positive off shore and negative on land). Default: False.
- precomputed : bool
- If True, use precomputed distances (from NASA). Default: False.
- scheduler : str
- Used for dask map_partitions. "threads", "synchronous" or "processes"
+ file_name : str
+ path of the file
+ src_crs : crs, optional
+ source CRS. Provide it if error without it.
+ window : rasterio.windows.Window, optional
+ window to read
+ geometry : list of shapely.geometry, optional
+ consider pixels only within these shapes
+ dst_crs : crs, optional
+ reproject to given crs
+ transform : rasterio.Affine
+ affine transformation to apply
+ wdith : float
+ number of lons for transform
+ height : float
+ number of lats for transform
+ resampling : rasterio.warp.Resampling optional
+ resampling function used for reprojection to dst_crs,
+ default: nearest
+ return_meta : bool, optional
+ default: False
+
+ Returns
+ -------
+ centr : Centroids
+ Centroids according to the given raster file
+ meta : dict, optional if return_meta is True
+ Raster meta (height, width, transform, crs).
"""
- if (not self.lat.size or not self.lon.size) and not self.meta:
- LOGGER.warning('No lat/lon, no meta, nothing to do!')
- return
- if precomputed:
- if not self.lat.size or not self.lon.size:
- self.set_meta_to_lat_lon()
- self.dist_coast = u_coord.dist_to_coast_nasa(
- self.lat, self.lon, highres=True, signed=signed)
- else:
- ne_geom = self._ne_crs_geom(scheduler)
- LOGGER.debug('Computing distance to coast for %s centroids.', str(self.lat.size))
- self.dist_coast = u_coord.dist_to_coast(ne_geom, signed=signed)
+ meta, _ = u_coord.read_raster(
+ file_name, [1], src_crs, window, geometry, dst_crs,
+ transform, width, height, resampling,
+ )
+ centr = cls.from_meta(meta)
+ return (centr, meta) if return_meta else centr
- def set_on_land(self, scheduler=None):
- """Set on_land attribute for every pixel or point.
+ @classmethod
+ def from_meta(cls, meta):
+ """initiate centroids from meta raster definition
Parameters
----------
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
- """
- ne_geom = self._ne_crs_geom(scheduler)
- LOGGER.debug('Setting on_land %s points.', str(self.lat.size))
- self.on_land = u_coord.coord_on_land(
- ne_geom.geometry[:].y.values, ne_geom.geometry[:].x.values)
-
- def remove_duplicate_points(self):
- """Return Centroids with removed duplicated points
+ meta : dict
+ meta description of raster
Returns
-------
- cen : Centroids
- Sub-selection of this object.
+ Centroid
+ Centroids initialized for raster described by meta.
"""
- if not self.lat.any() and not self.meta:
- return self
- if self.lat.size > 0:
- coords_view = self.coord.astype(np.float64).view(dtype='float64,float64')
- sel_cen = np.sort(np.unique(coords_view, return_index=True)[1])
- else:
- geom_wkb = self.geometry.apply(lambda geom: geom.wkb)
- sel_cen = geom_wkb.drop_duplicates().index
- return self.select(sel_cen=sel_cen)
+ crs = meta['crs']
+ lat, lon = _meta_to_lat_lon(meta)
+ return cls(lon=lon, lat=lat, crs=crs)
- def select(self, reg_id=None, extent=None, sel_cen=None):
- """Return Centroids with points in the given reg_id or within mask
+ @classmethod
+ def from_vector_file(cls, file_name, dst_crs=None):
+ """Create Centroids object from vector file (any format supported by fiona).
Parameters
----------
- reg_id : int
- region to filter according to region_id values
- extent : tuple
- Format (min_lon, max_lon, min_lat, max_lat) tuple.
- If min_lon > lon_max, the extend crosses the antimeridian and is
- [lon_max, 180] + [-180, lon_min]
- Borders are inclusive.
- sel_cen : np.array
- 1-dim mask, overrides reg_id and extent
+ file_name : str
+ vector file with format supported by fiona and 'geometry' field.
+ dst_crs : crs, optional
+ reproject to given crs
+ If no crs is given in the file, simply sets the crs.
Returns
-------
- cen : Centroids
- Sub-selection of this object
- """
- if sel_cen is None:
- sel_cen = self.select_mask(reg_id=reg_id, extent=extent)
-
- if not self.lat.size or not self.lon.size:
- self.set_meta_to_lat_lon()
-
- centr = Centroids.from_lat_lon(self.lat[sel_cen], self.lon[sel_cen], self.geometry.crs)
- if self.area_pixel.size:
- centr.area_pixel = self.area_pixel[sel_cen]
- if self.region_id.size:
- centr.region_id = self.region_id[sel_cen]
- if self.on_land.size:
- centr.on_land = self.on_land[sel_cen]
- if self.dist_coast.size:
- centr.dist_coast = self.dist_coast[sel_cen]
- return centr
-
- def select_mask(self, reg_id=None, extent=None):
+ centr : Centroids
+ Centroids with points according to the given vector file
"""
- Make mask of selected centroids
+
+ centroids = cls.from_geodataframe(gpd.read_file(file_name))
+ if dst_crs is not None:
+ if centroids.crs:
+ centroids.to_crs(dst_crs, inplace=True)
+ else:
+ centroids.gdf.set_crs(dst_crs, inplace=True)
+ return centroids
+
+ @classmethod
+ def from_csv(cls, file_path, **kwargs):
+ """Generate centroids from a CSV file with column names in var_names.
Parameters
----------
- reg_id : int
- region to filter according to region_id values
- extent : tuple
- Format (min_lon, max_lon, min_lat, max_lat) tuple.
- If min_lon > lon_max, the extend crosses the antimeridian and is
- [lon_max, 180] + [-180, lon_min]
- Borders are inclusive.
+ file_path : str
+ path to CSV file to be read
+ kwargs : dict
+ Additional keyword arguments to pass on to pandas.read_csv.
Returns
-------
- sel_cen : 1d array of booleans
- 1d mask of selected centroids
-
+ Centroids
"""
- sel_cen = np.ones(self.size, dtype=bool)
- if reg_id is not None:
- sel_cen &= np.isin(self.region_id, reg_id)
- if extent is not None:
- lon_min, lon_max, lat_min, lat_max = extent
- lon_max += 360 if lon_min > lon_max else 0
- lon_normalized = u_coord.lon_normalize(
- self.lon.copy(), center=0.5 * (lon_min + lon_max))
- sel_cen &= (
- (lon_normalized >= lon_min) & (lon_normalized <= lon_max) &
- (self.lat >= lat_min) & (self.lat <= lat_max)
- )
- return sel_cen
+ return cls._from_dataframe(pd.read_csv(file_path, **kwargs))
- def set_lat_lon_to_meta(self, min_resol=1.0e-8):
- """Compute meta from lat and lon values.
+ def write_csv(self, file_path):
+ """Save centroids as CSV file
Parameters
----------
- min_resol : float, optional
- Minimum centroids resolution to use in the raster. Default: 1.0e-8.
+ file_path : str, Path
+ absolute or relative file path and name to write to
"""
- res = u_coord.get_resolution(self.lon, self.lat, min_resol=min_resol)
- rows, cols, ras_trans = u_coord.pts_to_raster_meta(self.total_bounds, res)
- LOGGER.debug('Resolution points: %s', str(res))
- self.meta = {
- 'width': cols,
- 'height': rows,
- 'crs': self.crs,
- 'transform': ras_trans,
- }
+ file_path = Path(file_path).with_suffix('.csv')
+ LOGGER.info('Writing %s', file_path)
+ self._centroids_to_dataframe().to_csv(file_path, index=False)
- def set_meta_to_lat_lon(self):
- """Compute lat and lon of every pixel center from meta raster."""
- if self.meta:
- xgrid, ygrid = u_coord.raster_to_meshgrid(
- self.meta['transform'], self.meta['width'], self.meta['height'])
- self.lon = xgrid.flatten()
- self.lat = ygrid.flatten()
- self.geometry = gpd.GeoSeries(crs=self.meta['crs'])
- def plot(self, axis=None, figsize=(9, 13), **kwargs):
- """Plot centroids scatter points over earth.
+ @classmethod
+ def from_excel(cls, file_path, sheet_name=None):
+ """Generate a new centroids object from an excel file with column names in var_names.
Parameters
----------
- axis : matplotlib.axes._subplots.AxesSubplot, optional
- axis to use
- figsize: (float, float), optional
- figure size for plt.subplots
- The default is (9, 13)
- kwargs : optional
- arguments for scatter matplotlib function
+ file_path : str
+ absolute or relative file path
+ sheet_name : str, optional
+ name of sheet in excel file containing centroid information
+ Default: "centroids"
Returns
-------
- axis : matplotlib.axes._subplots.AxesSubplot
+ centr : Centroids
+ Centroids with data from the given excel file
"""
- if self.meta and not self.coord.size:
- self.set_meta_to_lat_lon()
- pad = np.abs(u_coord.get_resolution(self.lat, self.lon)).min()
+ if sheet_name is None:
+ sheet_name = 'centroids'
+ df = pd.read_excel(file_path, sheet_name)
+ return cls._from_dataframe(df)
- proj_data, _ = u_plot.get_transformation(self.crs)
- proj_plot = proj_data
- if isinstance(proj_data, ccrs.PlateCarree):
- # use different projections for plot and data to shift the central lon in the plot
- xmin, ymin, xmax, ymax = u_coord.latlon_bounds(self.lat, self.lon, buffer=pad)
- proj_plot = ccrs.PlateCarree(central_longitude=0.5 * (xmin + xmax))
- else:
- xmin, ymin, xmax, ymax = (self.lon.min() - pad, self.lat.min() - pad,
- self.lon.max() + pad, self.lat.max() + pad)
-
- if not axis:
- _, axis, _fontsize = u_plot.make_map(proj=proj_plot, figsize=figsize)
-
- axis.set_extent((xmin, xmax, ymin, ymax), crs=proj_data)
- u_plot.add_shapes(axis)
- axis.scatter(self.lon, self.lat, transform=proj_data, **kwargs)
- plt.tight_layout()
- return axis
-
- def calc_pixels_polygons(self, scheduler=None):
- """Return a gpd.GeoSeries with a polygon for every pixel
+ def write_excel(self, file_path):
+ """Save centroids as excel file
Parameters
----------
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
-
- Returns
- -------
- geo : gpd.GeoSeries
+ file_path : str, Path
+ absolute or relative file path and name to write to
"""
- if not self.meta:
- self.set_lat_lon_to_meta()
- if abs(abs(self.meta['transform'].a) -
- abs(self.meta['transform'].e)) > 1.0e-5:
- raise ValueError('Area can not be computed for not squared pixels.')
- self.set_geometry_points(scheduler)
- return self.geometry.buffer(self.meta['transform'].a / 2).envelope
-
- def empty_geometry_points(self):
- """Removes all points in geometry.
-
- Useful when centroids is used in multiprocessing function."""
- self.geometry = gpd.GeoSeries(crs=self.geometry.crs)
+ file_path = Path(file_path).with_suffix('.xlsx')
+ LOGGER.info('Writing %s', file_path)
+ self._centroids_to_dataframe().to_excel(
+ file_path, sheet_name=DEF_SHEET_NAME, index=False,
+ )
- def write_hdf5(self, file_data):
- """Write centroids attributes into hdf5 format.
+ def write_hdf5(self, file_name, mode='w'):
+ """Write data frame and metadata in hdf5 format
Parameters
----------
- file_data : str or h5
- If string, path to write data. If h5 object, the datasets will be generated there.
+ file_name : str
+ (path and) file name to write to.
"""
- if isinstance(file_data, str):
- LOGGER.info('Writing %s', file_data)
- with h5py.File(file_data, 'w') as data:
- self._write_hdf5(data)
- else:
- self._write_hdf5(file_data)
-
- def _write_hdf5(self, data):
- str_dt = h5py.special_dtype(vlen=str)
- for centr_name, centr_val in self.__dict__.items():
- if isinstance(centr_val, np.ndarray):
- data.create_dataset(centr_name, data=centr_val, compression="gzip")
- elif centr_name == 'meta' and centr_val:
- centr_meta = data.create_group(centr_name)
- for key, value in centr_val.items():
- if value is None:
- LOGGER.info("Skip writing Centroids.meta['%s'] for it is None.", key)
- elif key not in ('crs', 'transform'):
- if not isinstance(value, str):
- centr_meta.create_dataset(key, (1,), data=value, dtype=type(value))
- else:
- hf_str = centr_meta.create_dataset(key, (1,), dtype=str_dt)
- hf_str[0] = value
- elif key == 'transform':
- centr_meta.create_dataset(
- key, (6,),
- data=[value.a, value.b, value.c, value.d, value.e, value.f],
- dtype=float)
- elif centr_name == 'geometry':
- LOGGER.debug("Skip writing Centroids.geometry")
- else:
- LOGGER.info("Skip writing Centroids.%s:%s, it's neither an array nor a non-empty"
- " meta object", centr_name, centr_val.__class__.__name__)
- hf_str = data.create_dataset('crs', (1,), dtype=str_dt)
- hf_str[0] = CRS.from_user_input(self.crs).to_wkt()
+ LOGGER.info('Writing %s', file_name)
+ store = pd.HDFStore(file_name, mode=mode)
+ pandas_df = pd.DataFrame(self.gdf)
+ for col in pandas_df.columns:
+ if str(pandas_df[col].dtype) == "geometry":
+ pandas_df[col] = np.asarray(self.gdf[col])
+
+ # Avoid pandas PerformanceWarning when writing HDF5 data
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore", category=pd.errors.PerformanceWarning)
+ # Write dataframe
+ store.put('centroids', pandas_df)
+
+ store.get_storer('centroids').attrs.metadata = {
+ 'crs': CRS.from_user_input(self.crs).to_wkt()
+ }
+
+ store.close()
- def read_hdf5(self, *args, **kwargs):
- """This function is deprecated, use Centroids.from_hdf5 instead."""
- LOGGER.warning("The use of Centroids.read_hdf5 is deprecated."
- "Use Centroids.from_hdf5 instead.")
- self.__dict__ = Centroids.from_hdf5(*args, **kwargs).__dict__
@classmethod
- def from_hdf5(cls, file_data):
+ def from_hdf5(cls, file_name):
"""Create a centroids object from a HDF5 file.
Parameters
@@ -1211,30 +916,60 @@ def from_hdf5(cls, file_data):
-------
centr : Centroids
Centroids with data from the given file
+
+ Raises
+ ------
+ FileNotFoundError
"""
- if isinstance(file_data, (str, Path)):
- LOGGER.info('Reading %s', file_data)
- with h5py.File(file_data, 'r') as data:
- return cls._from_hdf5(data)
- else:
- return cls._from_hdf5(file_data)
+ if not Path(file_name).is_file():
+ raise FileNotFoundError(str(file_name))
+ try:
+ with pd.HDFStore(file_name, mode='r') as store:
+ metadata = store.get_storer('centroids').attrs.metadata
+ # in previous versions of CLIMADA and/or geopandas,
+ # the CRS was stored in '_crs'/'crs'
+ crs = metadata.get('crs')
+ gdf = gpd.GeoDataFrame(store['centroids'], crs=crs)
+ except TypeError:
+ with h5py.File(file_name, 'r') as data:
+ gdf = cls._gdf_from_legacy_hdf5(data.get('centroids'))
+ except KeyError:
+ with h5py.File(file_name, 'r') as data:
+ gdf = cls._gdf_from_legacy_hdf5(data)
+
+ return cls.from_geodataframe(gdf)
+
+ ##
+ # Private methods
+ ##
@classmethod
- def _from_hdf5(cls, data):
- centr = None
+ def _from_dataframe(cls, df):
+ if 'crs' in df.columns:
+ crs = df['crs'].iloc[0]
+ else:
+ LOGGER.info("No 'crs' column provided in file, setting CRS to WGS84 default.")
+ crs = DEF_CRS
+
+ extra_values = {
+ col: df[col]
+ for col in df.columns
+ if col not in ['lat', 'lon', 'crs']
+ }
+
+ return cls(lat=df['lat'], lon=df['lon'], **extra_values, crs=crs)
+
+ @staticmethod
+ def _gdf_from_legacy_hdf5(data):
crs = DEF_CRS
if data.get('crs'):
crs = u_coord.to_crs_user_input(data.get('crs')[0])
if data.get('lat') and data.get('lat').size:
- centr = cls.from_lat_lon(
- np.array(data.get('lat')),
- np.array(data.get('lon')),
- crs=crs)
+ latitude = np.array(data.get('lat'))
+ longitude = np.array(data.get('lon'))
elif data.get('latitude') and data.get('latitude').size:
- centr = cls.from_lat_lon(
- np.array(data.get('latitude')),
- np.array(data.get('longitude')),
- crs=crs)
+ latitude = np.array(data.get('latitude'))
+ longitude = np.array(data.get('longitude'))
else:
centr_meta = data.get('meta')
meta = dict()
@@ -1244,152 +979,130 @@ def _from_hdf5(cls, data):
meta[key] = value[0]
else:
meta[key] = rasterio.Affine(*value)
- centr = cls(meta=meta)
+ latitude, longitude = _meta_to_lat_lon(meta)
+ extra_values = {}
for centr_name in data.keys():
- if centr_name not in ('crs', 'lat', 'lon', 'meta'):
- setattr(centr, centr_name, np.array(data.get(centr_name)))
- return centr
-
- @property
- def crs(self):
- """Get CRS of raster or vector."""
- if self.meta:
- return self.meta['crs']
- if self.geometry.crs:
- return self.geometry.crs
- return DEF_CRS
-
- @property
- def size(self):
- """Get number of pixels or points."""
- if self.meta:
- return int(self.meta['height'] * self.meta['width'])
- return self.lat.size
+ if centr_name not in ('crs', 'lat', 'lon', 'meta', 'latitude', 'longitude'):
+ values = np.array(data.get(centr_name))
+ if latitude.size != 0 and values.size != 0 :
+ extra_values[centr_name] = values
+
+ return gpd.GeoDataFrame(
+ extra_values,
+ geometry=gpd.points_from_xy(x=longitude, y=latitude, crs=crs),
+ )
- @property
- def shape(self):
- """Get shape of rastered data."""
+ @classmethod
+ def _legacy_from_excel(cls, file_name, var_names):
+ LOGGER.info('Reading %s', file_name)
try:
- if self.meta:
- return (self.meta['height'], self.meta['width'])
- return (np.unique(self.lat).size, np.unique(self.lon).size)
- except AttributeError:
- return ()
-
- @property
- def total_bounds(self):
- """Get total bounds (left, bottom, right, top)."""
- if self.meta:
- left = self.meta['transform'].xoff
- right = left + self.meta['transform'][0] * self.meta['width']
- if left > right:
- left, right = right, left
- top = self.meta['transform'].yoff
- bottom = top + self.meta['transform'][4] * self.meta['height']
- if bottom > top:
- bottom, top = top, bottom
- return left, bottom, right, top
- return self.lon.min(), self.lat.min(), self.lon.max(), self.lat.max()
-
- @property
- def coord(self):
- """Get [lat, lon] array."""
- return np.stack([self.lat, self.lon], axis=1)
+ df = pd.read_excel(file_name, var_names['sheet_name'])
+ df = df.rename(columns=var_names['col_name'])
+ except KeyError as err:
+ raise KeyError("Not existing variable: %s" % str(err)) from err
+ return cls._from_dataframe(df)
- def set_geometry_points(self, scheduler=None):
- """Set `geometry` attribute with Points from `lat`/`lon` attributes.
+ def _centroids_to_dataframe(self):
+ """Create dataframe from Centroids object to facilitate
+ saving in different file formats.
- Parameters
- ----------
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
+ Returns
+ -------
+ df : DataFrame
"""
- def apply_point(df_exp):
- return df_exp.apply((lambda row: Point(row.longitude, row.latitude)), axis=1)
- if not self.geometry.size:
- LOGGER.info('Convert centroids to GeoSeries of Point shapes.')
- if (not self.lat.any() or not self.lon.any()) and self.meta:
- self.set_meta_to_lat_lon()
- if not scheduler:
- self.geometry = gpd.GeoSeries(
- gpd.points_from_xy(self.lon, self.lat), crs=self.geometry.crs)
- else:
- import dask.dataframe as dd
- from multiprocessing import cpu_count
- ddata = dd.from_pandas(self, npartitions=cpu_count())
- self.geometry = (ddata
- .map_partitions(apply_point, meta=Point)
- .compute(scheduler=scheduler))
-
- def _ne_crs_geom(self, scheduler=None):
+ df = pd.DataFrame(self.gdf)
+ df['lon'] = self.gdf['geometry'].x
+ df['lat'] = self.gdf['geometry'].y
+ df['crs'] = CRS.from_user_input(self.crs).to_wkt()
+ df = df.drop(['geometry'], axis=1)
+ return df
+
+ def _ne_crs_geom(self):
"""Return `geometry` attribute in the CRS of Natural Earth.
- Parameters
- ----------
- scheduler : str
- used for dask map_partitions. “threads”, “synchronous” or “processes”
-
Returns
-------
geo : gpd.GeoSeries
"""
- if not self.lat.size or not self.lon.size:
- self.set_meta_to_lat_lon()
- if u_coord.equal_crs(self.geometry.crs, u_coord.NE_CRS) and self.geometry.size:
- return self.geometry
- self.set_geometry_points(scheduler)
- return self.geometry.to_crs(u_coord.NE_CRS)
-
- def __deepcopy__(self, memo):
- """Avoid error deep copy in gpd.GeoSeries by setting only the crs."""
- cls = self.__class__
- result = cls.__new__(cls)
- memo[id(self)] = result
- for key, value in self.__dict__.items():
- if key == 'geometry':
- setattr(result, key, gpd.GeoSeries(crs=self.geometry.crs))
- else:
- setattr(result, key, copy.deepcopy(value, memo))
- return result
+ if u_coord.equal_crs(self.gdf.crs, u_coord.NE_CRS):
+ return self.gdf.geometry
+ return self.to_crs(u_coord.NE_CRS, inplace=False).geometry
+
+ ##
+ # Deprecated methods
+ ##
+
+ @classmethod
+ @deprecated(details="Reading Centroids data from matlab files is not supported anymore."
+ "This method has been removed with climada 5.0")
+ def from_mat(cls, file_name, var_names=None):
+ """Reading Centroids data from matlab files is not supported anymore.
+ This method has been removed with climada 5.0"""
+ raise NotImplementedError("You are suggested to use an old version of climada (<=4.*) and"
+ " convert the file to hdf5 format.")
+
+ @staticmethod
+ @deprecated(details="This method has been removed with climada 5.0")
+ def from_base_grid(land=False, res_as=360, base_file=None):
+ """This method has been removed with climada 5.0"""
+ raise NotImplementedError("Create the Centroids from a custom base file or from Natural"
+ " Earth (files are available in Climada, look up ``climada.util"
+ ".constants.NATEARTH_CENTROIDS`` for their location)")
+
+ @classmethod
+ @deprecated(details="This method will be removed in a future version."
+ " Simply use the constructor instead.")
+ def from_lat_lon(cls, lat, lon, crs="EPSG:4326"):
+ """deprecated, use the constructor instead"""
+ return cls(lat=lat, lon=lon, crs=crs)
+
+ @deprecated(details="This method is futile and will be removed in a future version."
+ " `Centroids.get_area_pixel` can be run without initialization.")
+ def set_area_pixel(self, min_resol=1e-08, scheduler=None):
+ """deprecated, obsolete"""
+
+ @deprecated(details="This method is futile and will be removed in a future version."
+ " `Centroids.get_area_pixel` can be run without initialization.")
+ def set_area_approx(self, min_resol=1e-08):
+ """deprecated, obsolete"""
+
+ @deprecated(details="This method is futile and will be removed in a future version."
+ " `Centroids.get_dist_coast` can be run without initialization.")
+ def set_dist_coast(self, signed=False, precomputed=False, scheduler=None):
+ """deprecated, obsolete"""
+
+ @deprecated(details="This method has no effect and will be removed in a future version."
+ " In the current version of climada the geometry points of a `Centroids` object"
+ " cannot be removed as they are the backbone of the Centroids' GeoDataFrame.")
+ def empty_geometry_points(self):
+ """"deprecated, has no effect, which may be unexpected: no geometry points will be removed,
+ the centroids' GeoDataFrame is built on them!
+ """
+
+ @deprecated(details="This method has no effect and will be removed in a future version.")
+ def set_meta_to_lat_lon(self):
+ """deprecated, has no effect"""
+ @deprecated(details="This method has no effect and will be removed in a future version.")
+ def set_lat_lon_to_meta(self, min_resol=1e-08):
+ """deprecated, has no effect"""
-def generate_nat_earth_centroids(res_as=360, path=None, dist_coast=False):
- """Generate hdf5 file containing Centroids of given resolution.
- For reproducibility, this is the function that generates the centroids files in
- `NATEARTH_CENTROIDS`. These files are provided with CLIMADA so that this function should never
- be called!
+def _meta_to_lat_lon(meta):
+ """Compute lat and lon of every pixel center from meta raster.
Parameters
----------
- res_as : int
- Resolution of file in arc-seconds. Default: 360.
- path : str, optional
- If set, write resulting hdf5 file here instead of the default location.
- dist_coast : bool, optional
- If True, read distance from a NASA dataset (see util.coordinates.dist_to_coast_nasa).
- Default: False.
+ meta : dict
+ meta description of raster
+
+ Returns
+ -------
+ latitudes : np.ndarray
+ Latitudinal coordinates of pixel centers.
+ longitudes : np.ndarray
+ Longitudinal coordinates of pixel centers.
"""
- if path is None and res_as not in [150, 360]:
- raise ValueError("Only 150 and 360 arc-seconds are supported!")
-
- res_deg = res_as / 3600
- lat_dim = np.arange(-90 + res_deg, 90, res_deg)
- lon_dim = np.arange(-180 + res_deg, 180 + res_deg, res_deg)
- lon, lat = [ar.ravel() for ar in np.meshgrid(lon_dim, lat_dim)]
- natids = np.uint16(u_coord.get_country_code(lat, lon, gridded=False))
-
- cen = Centroids.from_lat_lon(lat, lon)
- cen.region_id = natids
- cen.set_lat_lon_to_meta()
- cen.lat = np.array([])
- cen.lon = np.array([])
-
- if path is None:
- path = NATEARTH_CENTROIDS[res_as]
-
- if dist_coast:
- cen.set_dist_coast(precomputed=True, signed=False)
- cen.dist_coast = np.float16(cen.dist_coast)
- cen.write_hdf5(path)
+ xgrid, ygrid = u_coord.raster_to_meshgrid(meta['transform'], meta['width'], meta['height'])
+ return ygrid.ravel(), xgrid.ravel()
diff --git a/climada/hazard/centroids/test/test_centr.py b/climada/hazard/centroids/test/test_centr.py
index 5397b01a73..cb0f47aedd 100644
--- a/climada/hazard/centroids/test/test_centr.py
+++ b/climada/hazard/centroids/test/test_centr.py
@@ -19,115 +19,767 @@
Test CentroidsVector and CentroidsRaster classes.
"""
import unittest
+from unittest.mock import patch
from pathlib import Path
import numpy as np
import pandas as pd
import geopandas as gpd
+import shapely
+import itertools
+import rasterio
+from pyproj.crs.crs import CRS
+from rasterio.windows import Window
+from shapely.geometry.point import Point
+from cartopy.io import shapereader
+
+from climada import CONFIG
from climada.hazard.centroids.centr import Centroids
-from climada.util.constants import GLB_CENTROIDS_MAT, HAZ_TEMPLATE_XLS
-from climada.test import get_test_file
-
-HAZ_TEST_MAT = get_test_file('atl_prob_no_name', file_format='matlab')
-
-
-class TestCentroidsReader(unittest.TestCase):
- """Test read functions Centroids"""
-
- def test_mat_pass(self):
- """Read a centroid mat file correctly."""
- centroids = Centroids.from_mat(HAZ_TEST_MAT)
-
- n_centroids = 100
- self.assertEqual(centroids.coord.shape, (n_centroids, 2))
- self.assertEqual(centroids.coord[0][0], 21)
- self.assertEqual(centroids.coord[0][1], -84)
- self.assertEqual(centroids.coord[n_centroids - 1][0], 30)
- self.assertEqual(centroids.coord[n_centroids - 1][1], -75)
-
- def test_mat_global_pass(self):
- """Test read GLB_CENTROIDS_MAT"""
- centroids = Centroids.from_mat(GLB_CENTROIDS_MAT)
-
- self.assertEqual(centroids.region_id[1062443], 35)
- self.assertEqual(centroids.region_id[170825], 28)
-
- def test_centroid_pass(self):
- """Read a centroid excel file correctly."""
- centroids = Centroids.from_excel(HAZ_TEMPLATE_XLS)
-
- n_centroids = 45
- self.assertEqual(centroids.coord.shape[0], n_centroids)
- self.assertEqual(centroids.coord.shape[1], 2)
- self.assertEqual(centroids.coord[0][0], -25.95)
- self.assertEqual(centroids.coord[0][1], 32.57)
- self.assertEqual(centroids.coord[n_centroids - 1][0], -24.7)
- self.assertEqual(centroids.coord[n_centroids - 1][1], 33.88)
-
- def test_base_grid(self):
- """Read new centroids using from_base_grid, then select by extent."""
- centroids = Centroids.from_base_grid(land=True, res_as=150)
- self.assertEqual(centroids.lat.size, 8858035)
- self.assertTrue(np.all(np.diff(centroids.lat) <= 0))
-
- count_sandwich = np.sum(centroids.region_id == 239)
- self.assertEqual(count_sandwich, 321)
-
- count_sgi = centroids.select(
- reg_id=239,
- extent=(-39, -34.7, -55.5, -53.6) # south georgia island
- ).size
- self.assertEqual(count_sgi, 296)
-
- # test negative latitudinal orientation by testing that northern hemisphere (Russia)
- # is listed before southern hemisphere (South Africa)
- russia_max_idx = (centroids.region_id == 643).nonzero()[0].max()
- safrica_min_idx = (centroids.region_id == 710).nonzero()[0].min()
- self.assertTrue(russia_max_idx < safrica_min_idx)
-
- def test_geodataframe(self):
- """Test that constructing a valid Centroids instance from gdf works."""
- gdf = gpd.GeoDataFrame(pd.read_excel(HAZ_TEMPLATE_XLS))
- gdf.geometry = gpd.points_from_xy(
- gdf['longitude'], gdf['latitude']
+from climada.util.constants import DEF_CRS, HAZ_DEMO_FL
+import climada.util.coordinates as u_coord
+from climada.entity import Exposures
+from rasterio import Affine
+
+
+DATA_DIR = CONFIG.hazard.test_data.dir()
+
+# Note: the coordinates are not directly on the cities, the region id and on land
+# otherwise do not work correctly. It is only a close point.
+LATLON = np.array([
+ [-21.1736, -175.1883], # Tonga, Nuku'alofa, TON, 776
+ [-18.133, 178.433], # Fidji, Suva, FJI, 242 IN WATER IN NATURAL EARTH
+ [-38.4689, 177.8642], # New-Zealand, Te Karaka, NZL, 554
+ [69.6833, 18.95], # Norway, Tromso, NOR, 578 IN WATER IN NATURAL EARTH
+ [78.84422, 20.82842], # Norway, Svalbard, NOR, 578
+ [1, 1], # Ocean, 0 (0,0 is onland in Natural earth for testing reasons)
+ [-77.85, 166.6778], # Antarctica, McMurdo station, ATA, 010
+ [-0.25, -78.5833] # Ecuador, Quito, ECU, 218
+])
+
+VEC_LAT = LATLON[:,0]
+VEC_LON = LATLON[:,1]
+
+ON_LAND = np.array([True, False, True, False, True, False, True, True])
+REGION_ID = np.array([776, 0, 554, 0, 578, 0, 10, 218])
+
+TEST_CRS = 'EPSG:4326'
+ALT_CRS = 'epsg:32632' # UTM zone 32N (Central Europe, 6-12°E)
+
+class TestCentroidsData(unittest.TestCase):
+ """ Test class for initialisation and value based creation of Centroids objects"""
+ def setUp(self):
+ self.lat = np.array([-10, 0, 10])
+ self.lon = np.array([-170, -150, -130])
+ self.region_id = np.array([1, 2, 3])
+ self.on_land = np.array([True, False, False])
+ self.crs = 'epsg:32632'
+ self.centr = Centroids(lat=VEC_LAT,lon=VEC_LON)
+
+ def test_centroids_check_pass(self):
+ """Test vector data in Centroids"""
+ centr = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS)
+
+ self.assertTrue(u_coord.equal_crs(centr.crs, CRS.from_user_input(ALT_CRS)))
+ self.assertEqual(
+ list(centr.total_bounds),
+ [VEC_LON.min(), VEC_LAT.min(), VEC_LON.max(), VEC_LAT.max()],
+ )
+
+ self.assertIsInstance(centr,Centroids)
+ self.assertIsInstance(centr.lat, np.ndarray)
+ self.assertIsInstance(centr.lon, np.ndarray)
+ self.assertIsInstance(centr.coord, np.ndarray)
+ self.assertTrue(np.array_equal(centr.lat, VEC_LAT))
+ self.assertTrue(np.array_equal(centr.lon, VEC_LON))
+ self.assertTrue(np.array_equal(centr.coord, np.array([VEC_LAT, VEC_LON]).transpose()))
+ self.assertEqual(centr.size, VEC_LON.size)
+
+ def test_init_pass(self):
+ # Creating Centroids with latitude and longitude arrays
+ # check instance - trivial...
+ # Checking attributes
+ np.testing.assert_array_equal(self.centr.lat, VEC_LAT)
+ np.testing.assert_array_equal(self.centr.lon, VEC_LON)
+ self.assertTrue(u_coord.equal_crs(self.centr.crs, DEF_CRS))
+
+ # Creating Centroids with additional attributes
+ centroids = Centroids(lat=VEC_LAT, lon=VEC_LON,
+ region_id=REGION_ID, on_land=ON_LAND)
+
+ # Checking additional attributes
+ np.testing.assert_array_equal(centroids.region_id, REGION_ID)
+ np.testing.assert_array_equal(centroids.on_land, ON_LAND)
+
+ def test_init_defaults(self):
+ ''' Checking default values for Centroids'''
+ centroids = Centroids(lat=VEC_LAT,lon=VEC_LON)
+ # Checking defaults: nothing set for region_id, on_land
+ self.assertFalse(centroids.region_id)
+ self.assertFalse(centroids.on_land)
+ # Guarantee a no-default TypeError for lon/lat
+ with self.assertRaises(TypeError): Centroids()
+
+ def test_init_properties(self):
+ """ Guarantee that Centroid objects have at least the properties: """
+ properties = ['gdf','lon','lat','geometry',
+ 'on_land','region_id','crs',
+ 'shape','size','total_bounds','coord']
+ centroids = Centroids(lat=[],lon=[])
+ [self.assertTrue(hasattr(centroids,prop)) for prop in properties]
+
+ def test_init_kwargs(self):
+ """ Test default crs and kwargs forwarding """
+ centr = Centroids(
+ lat=VEC_LAT,
+ lon=VEC_LON,
+ region_id=REGION_ID,
+ on_land=ON_LAND,
+ )
+ self.assertTrue(u_coord.equal_crs(centr.crs, DEF_CRS))
+ self.assertTrue(np.allclose(centr.region_id, REGION_ID))
+ self.assertTrue(np.allclose(centr.on_land, ON_LAND))
+
+ # make sure kwargs are properly forwarded to centroids.gdf
+ np.random.seed(1000)
+ randommask = np.random.choice([True,False],size=len(VEC_LON))
+ centroids = Centroids(lat=VEC_LAT,lon=VEC_LON,masked=randommask,ones=1)
+ self.assertTrue(hasattr(centroids.gdf,'masked'))
+ self.assertTrue(hasattr(centroids.gdf,'ones'))
+ np.testing.assert_array_equal(randommask,centroids.gdf.masked)
+ self.assertEqual(sum(centroids.gdf.ones),len(VEC_LON))
+
+ def test_from_meta_pass(self):
+ expected_lon = np.array([-30.0, -20.0, -10.0]*3)
+ expected_lat = np.repeat([30.0, 20.0, 10.0],3)
+ # Check metadata
+ meta = dict(
+ crs=DEF_CRS,
+ height=3,
+ width=3,
+ transform=Affine(
+ 10, 0, -35,
+ 0, -10, 35,
+ ),
+ )
+ centroids = Centroids.from_meta(meta)
+
+ # check created object
+ np.testing.assert_array_equal(centroids.lon,expected_lon)
+ np.testing.assert_array_equal(centroids.lat,expected_lat)
+ self.assertEqual(centroids.crs,DEF_CRS)
+ # generally we assume that from_meta does not set region_ids and on_land flags
+ self.assertFalse(centroids.region_id)
+ self.assertFalse(centroids.on_land)
+
+ def test_from_meta(self):
+ """Test from_meta"""
+ meta_ref = {
+ 'width': 10,
+ 'height': 8,
+ 'transform': rasterio.Affine(
+ 0.6, 0, -0.1,
+ 0, -0.6, 0.3,
+ ),
+ 'crs': DEF_CRS,
+ }
+
+ lon_ref = np.array([0.2, 0.8, 1.4, 2.0, 2.6, 3.2, 3.8, 4.4, 5.0, 5.6])
+ lat_ref = np.array([0.0, -0.6, -1.2, -1.8, -2.4, -3.0, -3.6, -4.2])
+ lon_ref, lat_ref = [ar.ravel() for ar in np.meshgrid(lon_ref, lat_ref)]
+
+ centr = Centroids.from_meta(meta_ref)
+ meta = centr.get_meta()
+ self.assertTrue(u_coord.equal_crs(meta_ref["crs"], meta["crs"]))
+ self.assertEqual(meta_ref["width"], meta["width"])
+ self.assertEqual(meta_ref["height"], meta["height"])
+ np.testing.assert_allclose(meta_ref["transform"], meta["transform"])
+
+ centr = Centroids.from_meta(
+ Centroids(lat=lat_ref, lon=lon_ref).get_meta()
)
- gdf['elevation'] = np.random.rand(gdf.geometry.size)
- gdf['region_id'] = np.zeros(gdf.geometry.size)
- gdf['region_id'][0] = np.NaN
- gdf['geom'] = gdf.geometry # this should have no effect on centroids
+ np.testing.assert_allclose(lat_ref, centr.lat)
+ np.testing.assert_allclose(lon_ref, centr.lon)
+
+ # `get_meta` enforces same resolution in x and y, and y-coordinates are decreasing.
+ # For other cases, `from_meta` needs to be checked manually.
+ meta_ref = {
+ 'width': 4,
+ 'height': 5,
+ 'transform': rasterio.Affine(
+ 0.5, 0, 0.2,
+ 0, 0.6, -0.7,
+ ),
+ 'crs': DEF_CRS,
+ }
+ lon_ref = np.array([0.45, 0.95, 1.45, 1.95])
+ lat_ref = np.array([-0.4, 0.2, 0.8, 1.4, 2.0])
+ lon_ref, lat_ref = [ar.ravel() for ar in np.meshgrid(lon_ref, lat_ref)]
+
+ centr = Centroids.from_meta(meta_ref)
+ np.testing.assert_allclose(lat_ref, centr.lat)
+ np.testing.assert_allclose(lon_ref, centr.lon)
+
+
+ def test_from_pnt_bounds(self):
+ """Test from_pnt_bounds"""
+ width, height = 26, 51
+ left, bottom, right, top = 5, 0, 10, 10
+
+ centr = Centroids.from_pnt_bounds((left, bottom, right, top), 0.2, crs=DEF_CRS)
+ self.assertTrue(u_coord.equal_crs(centr.crs, DEF_CRS))
+ self.assertEqual(centr.size, width * height)
+ np.testing.assert_allclose([5.0, 5.2, 5.0], centr.lon[[0, 1, width]], atol=0.1)
+ np.testing.assert_allclose([10.0, 10.0, 9.8], centr.lat[[0, 1, width]], atol=0.1)
+ # generally we assume that from_meta does not set region_ids and on_land flags
+ self.assertFalse(centr.region_id)
+ self.assertFalse(centr.on_land)
+
+class TestCentroidsTransformation(unittest.TestCase):
+ """ Test class for coordinate transformations of Centroid objects
+ and modifications using set_ methods"""
+ def setUp(self):
+ self.lat = np.array([-10, 0, 10])
+ self.lon = np.array([-170, -150, -130])
+ self.region_id = np.array([1, 2, 3])
+ self.on_land = np.array([True, False, False])
+ self.crs = 'epsg:32632'
+ self.centr = Centroids(lat=VEC_LAT,lon=VEC_LON,crs=TEST_CRS)
+
+ def test_to_default_crs(self):
+ # Creating Centroids with non-default CRS and
+ # inplace transformation afterwards
+ centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS)
+ self.assertTrue(u_coord.equal_crs(centroids.crs, ALT_CRS))
+ centroids.to_default_crs()
+ # make sure CRS is DEF_CRS after transformation
+ self.assertTrue(u_coord.equal_crs(centroids.crs, DEF_CRS))
+ # Checking that modification actually took place
+ [self.assertNotEqual(x-y,0) for x,y in zip(centroids.lon,VEC_LON)]
+ [self.assertNotEqual(x-y,0) for x,y in zip(centroids.lat,VEC_LAT) if not x == 0]
+
+ def test_to_default_crs_not_inplace(self):
+ centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS)
+ newcentr = centroids.to_default_crs(inplace=False)
+ # make sure that new object has been created
+ self.assertIsNot(centroids,newcentr)
+ self.assertIsInstance(newcentr,Centroids)
+ ## compare with inplace transformation
+ centroids.to_default_crs()
+ np.testing.assert_array_equal(centroids.lat,newcentr.lat)
+ np.testing.assert_array_equal(centroids.lon,newcentr.lon)
+
+ def test_to_crs(self):
+ # Creating Centroids with default CRS
+ centroids = Centroids(lat=self.lat, lon=self.lon, crs=DEF_CRS)
+
+ # Transforming to another CRS
+ new_crs = 'epsg:3857'
+ transformed_centroids = centroids.to_crs(new_crs)
+
+ self.assertIsNot(centroids,transformed_centroids)
+ self.assertFalse(centroids == transformed_centroids)
+
+ # Checking CRS string after transformation
+ self.assertTrue(u_coord.equal_crs(transformed_centroids.crs, new_crs))
+ self.assertTrue(u_coord.equal_crs(centroids.crs, DEF_CRS))
+
+ # Checking correctness of transformation
+ expected_lat = np.array([-1118889.974858, 0., 1118889.9748585])
+ expected_lon = np.array([-18924313.434857, -16697923.618991, -14471533.803126])
+ np.testing.assert_array_almost_equal(transformed_centroids.lat, expected_lat)
+ np.testing.assert_array_almost_equal(transformed_centroids.lon, expected_lon)
+
+ def test_to_crs_inplace(self):
+ centroids = Centroids(lat=self.lat,lon=self.lon,crs=DEF_CRS)
+ new_crs = 'epsg:3857'
+ transformed_centroids = centroids.to_crs(new_crs)
+
+ # inplace transforming to another CRS
+ centroids.to_crs(new_crs,inplace=True)
+
+ self.assertTrue(centroids == transformed_centroids)
+
+ expected_lat = np.array([-1118889.974858, 0., 1118889.9748585])
+ expected_lon = np.array([-18924313.434857, -16697923.618991, -14471533.803126])
+ np.testing.assert_array_almost_equal(centroids.lat, expected_lat)
+ np.testing.assert_array_almost_equal(centroids.lon, expected_lon)
+
+ def test_ne_crs_geom_pass(self):
+ """Test _ne_crs_geom"""
+ natural_earth_geom = self.centr._ne_crs_geom()
+ self.assertEqual(natural_earth_geom.crs, u_coord.NE_CRS)
+
+ centr = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS)
+ ne_geom = centr._ne_crs_geom()
+ self.assertTrue(u_coord.equal_crs(ne_geom.crs, u_coord.NE_CRS))
+ np.testing.assert_allclose(ne_geom.geometry[:].x.values, 4.5, atol=0.1)
+ np.testing.assert_allclose(ne_geom.geometry[:].y.values, 0.0, atol=0.001)
+
+ def test_set_on_land_pass(self):
+ """Test set_on_land"""
+ self.centr.set_on_land()
+ np.testing.assert_array_equal(self.centr.on_land, ON_LAND)
+
+ centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, on_land='natural_earth')
+ np.testing.assert_array_equal(centroids.on_land, ON_LAND)
+
+ def test_set_on_land_implementationerror(self):
+ centroids = Centroids(lat=self.lat,lon=self.lon)
+
+ with self.assertRaises(NotImplementedError):
+ centroids.set_on_land(source='satellite',overwrite=True)
+
+ def test_set_on_land_raster(self):
+ """Test set_on_land"""
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
+ centr_ras.set_on_land()
+ self.assertTrue(np.array_equal(centr_ras.on_land, np.ones(60 * 50, bool)))
+
+ def test_set_region_id_pass(self):
+ """Test set_region_id"""
+ self.centr.set_region_id()
+ np.testing.assert_array_equal(self.centr.region_id, REGION_ID)
+
+ centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, region_id='country')
+ np.testing.assert_array_equal(centroids.region_id, REGION_ID)
+
+ def test_set_region_id_raster(self):
+ """Test set_region_id raster file"""
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
+ centr_ras.set_region_id()
+ self.assertEqual(centr_ras.region_id.size, centr_ras.size)
+ self.assertTrue(np.array_equal(np.unique(centr_ras.region_id), np.array([862])))
+
+ def test_set_region_id_implementationerror(self):
+ centroids = Centroids(lat=self.lat,lon=self.lon)
+
+ with self.assertRaises(NotImplementedError):
+ centroids.set_region_id(level='continent',overwrite=True)
+
+ def test_set_geometry_points_pass(self):
+ """Test set_geometry_points"""
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
+ x_flat = np.arange(-69.3326495969998, -68.88264959699978, 0.009000000000000341)
+ y_flat = np.arange(10.423720966978939, 9.883720966978919, -0.009000000000000341)
+ x_grid, y_grid = np.meshgrid(x_flat, y_flat)
+ self.assertTrue(np.allclose(x_grid.flatten(), centr_ras.lon))
+ self.assertTrue(np.allclose(y_grid.flatten(), centr_ras.lat))
+
+
+class TestCentroidsReaderWriter(unittest.TestCase):
+ """Test class for file based creation of Centroid objects and output"""
+
+ def test_from_csv_def_crs(self):
+ """Read a centroid csv file correctly and use default CRS."""
+ # Create temporary csv file containing centroids data
+ tmpfile = Path('test_write_csv.csv')
+ lat = np.array([0, 90, -90, 0, 0])
+ lon = np.array([0, 0, 0, 180, -180])
+ df = pd.DataFrame({'lat': lat, 'lon': lon})
+ df.to_csv(tmpfile, index=False)
+
+ # Read centroids using from_csv method
+ centroids = Centroids.from_csv(tmpfile)
+
+ # test attributes
+ np.testing.assert_array_equal(centroids.lat, lat)
+ np.testing.assert_array_equal(centroids.lon, lon)
+ self.assertEqual(centroids.crs, DEF_CRS)
+
+ # delete file
+ tmpfile.unlink()
+
+ def test_from_csv(self):
+ """Read a centroid csv file which contains CRS information."""
+ tmpfile = Path('test_write_csv.csv')
+ lat = np.array([0, 20048966.1, -20048966, 0, 0])
+ lon = np.array([0, 0, 0, 20037508.34, -20037508.34])
+ region_id = np.array([1, 2, 3, 4, 5])
+ on_land = np.array([True, False, False, True, True])
+ df = pd.DataFrame({'lat': lat, 'lon': lon, 'region_id': region_id, 'on_land': on_land})
+ df['crs'] = CRS.from_user_input(3857).to_wkt()
+ df.to_csv(tmpfile, index=False)
+
+ # Read centroids using from_csv method
+ centroids = Centroids.from_csv(tmpfile)
+
+ # Test attributes
+ np.testing.assert_array_equal(centroids.lat, lat)
+ np.testing.assert_array_equal(centroids.lon, lon)
+ self.assertEqual(centroids.crs, 'epsg:3857')
+ np.testing.assert_array_equal(centroids.region_id, region_id)
+ np.testing.assert_array_equal(centroids.on_land, on_land)
+
+ # Delete file
+ tmpfile.unlink()
+
+ def test_write_read_csv(self):
+ """Write and read a Centroids CSV file correctly."""
+ # Create Centroids with latitude and longitude arrays
+ tmpfile = Path('test_write_csv.csv')
+ lat = np.array([10.0, 20.0, 30.0])
+ lon = np.array([-10.0, -20.0, -30.0])
+ region_id = np.array([1, 2, 3])
+ on_land = np.array([True, False, False])
+ centroids_out = Centroids(lat=lat, lon=lon, region_id=region_id, on_land=on_land)
+
+ # Write CSV file from Centroids using write_csv
+ centroids_out.write_csv(tmpfile)
+
+ # Read CSV file using read_csv
+ centroids_in = Centroids.from_csv(tmpfile)
+
+ # Test attributes
+ np.testing.assert_array_equal(centroids_in.lat, centroids_out.lat)
+ np.testing.assert_array_equal(centroids_in.lon, centroids_out.lon)
+ self.assertTrue(u_coord.equal_crs(centroids_in.crs, centroids_out.crs))
+ np.testing.assert_array_equal(centroids_in.region_id, centroids_out.region_id)
+ np.testing.assert_array_equal(centroids_in.on_land, centroids_out.on_land)
+
+ # delete file
+ tmpfile.unlink()
+
+ def test_from_excel_def_crs(self):
+ """Read a centroid excel file correctly and use default CRS."""
+ # Create temporary excel file containing centroids data
+ tmpfile = Path('test_write_excel.xlsx')
+ lat = np.array([0, 90, -90, 0, 0])
+ lon = np.array([0, 0, 0, 180, -180])
+ df = pd.DataFrame({'lat': lat, 'lon': lon})
+ df.to_excel(tmpfile, sheet_name='centroids', index=False)
+
+ # Read centroids using from_excel method
+ centroids = Centroids.from_excel(file_path=tmpfile)
+
+ # test attributes
+ np.testing.assert_array_equal(centroids.lat, lat)
+ np.testing.assert_array_equal(centroids.lon, lon)
+ self.assertEqual(centroids.crs, DEF_CRS)
+
+ # delete file
+ tmpfile.unlink()
+
+ def test_from_excel(self):
+ """Read a centroid excel file correctly which contains CRS information."""
+ # Create temporary excel file containing centroids data
+ tmpfile = Path('test_write_excel.xlsx')
+ lat = np.array([0, 20048966.1, -20048966, 0, 0])
+ lon = np.array([0, 0, 0, 20037508.34, -20037508.34])
+ region_id = np.array([1, 2, 3, 4, 5])
+ on_land = np.array([True, False, False, True, True])
+ df = pd.DataFrame({'lat': lat, 'lon': lon, 'region_id': region_id, 'on_land': on_land})
+ df['crs'] = CRS.from_user_input(3857).to_wkt()
+ df.to_excel(tmpfile, sheet_name='centroids', index=False)
+
+ # Read centroids using from_excel method
+ centroids = Centroids.from_excel(file_path=tmpfile)
+
+ # test attributes
+ np.testing.assert_array_equal(centroids.lat, lat)
+ np.testing.assert_array_equal(centroids.lon, lon)
+ self.assertEqual(centroids.crs, 'epsg:3857')
+ np.testing.assert_array_equal(centroids.region_id, region_id)
+ np.testing.assert_array_equal(centroids.on_land, on_land)
+
+ # delete file
+ tmpfile.unlink()
+
+ def test_write_read_excel(self):
+ """Write and read a Centroids Excel file correctly."""
+ # Create Centroids with latitude and longitude arrays
+ tmpfile = Path('test_write_excel.xlsx')
+ lat = np.array([10.0, 20.0, 30.0])
+ lon = np.array([-10.0, -20.0, -30.0])
+ region_id = np.array([1, 2, 3])
+ on_land = np.array([True, False, False])
+ centroids_out = Centroids(lat=lat, lon=lon, region_id=region_id, on_land=on_land)
+
+ # Write Excel file from Centroids using write_csv
+ centroids_out.write_excel(tmpfile)
+
+ # Read Excel file using read_csv
+ centroids_in = Centroids.from_excel(tmpfile)
+
+ # Test attributes
+ np.testing.assert_array_equal(centroids_in.lat, centroids_out.lat)
+ np.testing.assert_array_equal(centroids_in.lon, centroids_out.lon)
+ self.assertTrue(u_coord.equal_crs(centroids_in.crs, centroids_out.crs))
+ np.testing.assert_array_equal(centroids_in.region_id, centroids_out.region_id)
+ np.testing.assert_array_equal(centroids_in.on_land, centroids_out.on_land)
+
+ # delete file
+ tmpfile.unlink()
+
+ def test_from_raster_file(self):
+ """Test from_raster_file"""
+ width, height = 50, 60
+ o_lat, o_lon = (10.42822096697894, -69.33714959699981)
+ res_lat, res_lon = (-0.009000000000000341, 0.009000000000000341)
+
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, width, height))
+ self.assertTrue(u_coord.equal_crs(centr_ras.crs, DEF_CRS))
+ self.assertEqual(centr_ras.size, width * height)
+ np.testing.assert_allclose(
+ [-69.333, -69.324, -69.333], centr_ras.lon[[0, 1, width]], atol=0.001,
+ )
+ np.testing.assert_allclose(
+ [10.424, 10.424, 10.415], centr_ras.lat[[0, 1, width]], atol=0.001,
+ )
+
+ def test_from_vector_file(self):
+ """Test from_vector_file and values_from_vector_files"""
+ shp_file = shapereader.natural_earth(resolution='110m', category='cultural',
+ name='populated_places_simple')
+
+ centr = Centroids.from_vector_file(shp_file, dst_crs=DEF_CRS)
+ self.assertTrue(u_coord.equal_crs(centr.crs, DEF_CRS))
+ self.assertAlmostEqual(centr.lon[0], 12.453386544971766)
+ self.assertAlmostEqual(centr.lon[-1], 114.18306345846304)
+ self.assertAlmostEqual(centr.lat[0], 41.903282179960115)
+ self.assertAlmostEqual(centr.lat[-1], 22.30692675357551)
+
+ centr = Centroids.from_vector_file(shp_file, dst_crs=ALT_CRS)
+ self.assertTrue(u_coord.equal_crs(centr.crs, ALT_CRS))
+
+ def test_from_geodataframe(self):
+ """Test that constructing a valid Centroids instance from gdf works."""
+ crs = DEF_CRS
+ lat = np.arange(170, 180)
+ lon = np.arange(-50, -40)
+ region_id = np.arange(1, 11)
+ on_land = np.ones(10, dtype=bool)
+ extra = np.full(10, 'a')
+
+ gdf = gpd.GeoDataFrame({
+ 'geometry': gpd.points_from_xy(lon, lat),
+ 'region_id': region_id,
+ 'on_land': on_land,
+ 'extra': extra,
+ }, crs=crs)
centroids = Centroids.from_geodataframe(gdf)
- centroids.check()
-
- self.assertEqual(centroids.geometry.size, 45)
- self.assertEqual(centroids.lon[0], 32.57)
- self.assertEqual(centroids.lat[0], -25.95)
- self.assertEqual(centroids.elevation.size, 45)
- self.assertEqual(centroids.on_land.sum(), 44)
- self.assertIsInstance(centroids.geometry, gpd.GeoSeries)
- self.assertIsInstance(centroids.geometry.total_bounds, np.ndarray)
-
-
-class TestCentroidsWriter(unittest.TestCase):
-
- def test_write_hdf5(self):
- tmpfile = 'test_write_hdf5.out.hdf5'
- xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon = 5, 6.5, -0.08, 0.12, 4, 5
- centr = Centroids.from_pix_bounds(xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon)
- with self.assertLogs('climada.hazard.centroids.centr', level='INFO') as cm:
- centr.write_hdf5(tmpfile)
- self.assertEqual(1, len(cm.output))
- self.assertIn(f"Writing {tmpfile}", cm.output[0])
- centr.meta['nodata'] = None
- with self.assertLogs('climada.hazard.centroids.centr', level='INFO') as cm:
- centr.write_hdf5(tmpfile)
- self.assertEqual(2, len(cm.output))
- self.assertIn("Skip writing Centroids.meta['nodata'] for it is None.", cm.output[1])
- Path(tmpfile).unlink()
+ for name, array in zip(
+ ['lat', 'lon', 'region_id', 'on_land'],
+ [lat, lon, region_id, on_land],
+ ):
+ np.testing.assert_array_equal(array, getattr(centroids, name))
+ self.assertTrue('extra' in centroids.gdf.columns)
+ self.assertTrue(u_coord.equal_crs(centroids.crs, crs))
+
+ def test_from_geodataframe_invalid(self):
+
+ # Creating an invalid GeoDataFrame with geometries that are not points
+ invalid_geometry_gdf = gpd.GeoDataFrame({
+ 'geometry': [
+ shapely.Point((2,2)),
+ shapely.Polygon([(0, 0), (1, 1), (1, 0), (0, 0)]),
+ shapely.LineString([(0, 1), (1, 0)]),
+ ],
+ })
+
+ with self.assertRaises(ValueError):
+ # Trying to create Centroids from invalid GeoDataFrame
+ Centroids.from_geodataframe(invalid_geometry_gdf)
+
+ def test_from_exposures_with_region_id(self):
+ """
+ Test that the `from_exposures` method correctly extracts
+ centroids and region_id from an `Exposure` object with region_id.
+ """
+ # Create an Exposure object with region_id, on_land and custom crs
+ lat = np.array([10.0, 20.0, 30.0])
+ lon = np.array([-10.0, -20.0, -30.0])
+ value = np.array([1, 1, 1])
+ region_id = np.array([1, 2, 3])
+ on_land = [False, True, True]
+ crs = 'epsg:32632'
+ gdf = gpd.GeoDataFrame({
+ 'latitude': lat,
+ 'longitude': lon,
+ 'value': value,
+ 'region_id': region_id,
+ 'on_land': on_land,
+ })
+ exposures = Exposures(gdf, crs=crs)
+
+ # Extract centroids from exposures
+ centroids = Centroids.from_exposures(exposures)
+
+ # Check attributes
+ np.testing.assert_array_equal(centroids.lat, lat)
+ np.testing.assert_array_equal(centroids.lon, lon)
+ np.testing.assert_array_equal(centroids.region_id, region_id)
+ np.testing.assert_array_equal(centroids.on_land, on_land)
+ self.assertFalse(np.isin('value', centroids.gdf.columns))
+ self.assertEqual(centroids.crs, crs)
+
+ def test_from_exposures_without_region_id(self):
+ """
+ Test that the `from_exposures` method correctly extracts
+ centroids from an `Exposure` object without region_id.
+ """
+ # Create an Exposure object without region_id and variables to ignore
+ # and default crs
+ lat = np.array([10.0, 20.0, 30.0])
+ lon = np.array([-10.0, -20.0, -30.0])
+ value = np.array([1, 1, 1])
+ impf_TC = np.array([1, 2, 3])
+ centr_TC = np.array([1, 2, 3])
+ gdf = gpd.GeoDataFrame({
+ 'latitude': lat,
+ 'longitude': lon,
+ 'value': value,
+ 'impf_tc': impf_TC,
+ 'centr_TC': centr_TC,
+ })
+ exposures = Exposures(gdf)
+
+ # Extract centroids from exposures
+ centroids = Centroids.from_exposures(exposures)
+
+ # Check attributes
+ self.assertEqual(centroids.lat.tolist(), lat.tolist())
+ self.assertEqual(centroids.lon.tolist(), lon.tolist())
+ self.assertTrue(u_coord.equal_crs(centroids.crs, DEF_CRS))
+ self.assertEqual(centroids.region_id, None)
+ self.assertEqual(centroids.on_land, None)
+ np.testing.assert_equal(
+ np.isin(['value', 'impf_tc', 'centr_tc'], centroids.gdf.columns),
+ False,
+ )
+
+ def test_from_exposure_exceptions(self):
+ gdf = gpd.GeoDataFrame({
+ })
+ exposures = Exposures(gdf)
+ with self.assertRaises(ValueError):
+ Centroids.from_exposures(exposures)
+
+ def test_read_write_hdf5(self):
+ tmpfile = Path('test_write_hdf5.out.hdf5')
+ crs = DEF_CRS
+ centroids_w = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=crs)
+ centroids_w.write_hdf5(tmpfile)
+ centroids_r = Centroids.from_hdf5(tmpfile)
+ self.assertTrue(centroids_w == centroids_r)
+ tmpfile.unlink()
+
+ def test_from_hdf5_nonexistent_file(self):
+ """Test raising FileNotFoundError when creating Centroids object from a nonexistent HDF5 file"""
+ file_name = "/path/to/nonexistentfile.h5"
+ # prescribe that file does not exist
+ with patch("pathlib.Path.is_file", return_value=False):
+ with self.assertRaises(FileNotFoundError):
+ Centroids.from_hdf5(file_name)
class TestCentroidsMethods(unittest.TestCase):
+ """Test Centroids methods"""
+ def setUp(self):
+ self.centr = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=TEST_CRS)
+
+ def test_select_pass(self):
+ """Test Centroids.select method"""
+ region_id = np.zeros(VEC_LAT.size)
+ region_id[[2, 4]] = 10
+ centr = Centroids(lat=VEC_LAT, lon=VEC_LON, region_id=region_id)
+
+ fil_centr = centr.select(reg_id=10)
+ self.assertIsInstance(fil_centr,Centroids)
+ self.assertEqual(fil_centr.size, 2)
+ self.assertEqual(fil_centr.lat[0], VEC_LAT[2])
+ self.assertEqual(fil_centr.lat[1], VEC_LAT[4])
+ self.assertEqual(fil_centr.lon[0], VEC_LON[2])
+ self.assertEqual(fil_centr.lon[1], VEC_LON[4])
+ self.assertTrue(np.array_equal(fil_centr.region_id, np.ones(2) * 10))
+
+ def test_select_extent_pass(self):
+ """Test select extent"""
+ centr = Centroids(
+ lat=np.array([-5, -3, 0, 3, 5]),
+ lon=np.array([-180, -175, -170, 170, 175]),
+ region_id=np.zeros(5),
+ )
+ ext_centr = centr.select(extent=[-175, -170, -5, 5])
+ self.assertIsInstance(ext_centr,Centroids)
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-175, -170]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-3, 0]))
+
+ # Cross antimeridian, version 1
+ ext_centr = centr.select(extent=[170, -175, -5, 5])
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
+
+ # Cross antimeridian, version 2
+ ext_centr = centr.select(extent=[170, 185, -5, 5])
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
+
+ def test_append_pass(self):
+ """Append points"""
+ centr = self.centr
+ centr_bis = Centroids(lat=np.array([1, 2, 3]), lon=np.array([4, 5, 6]), crs=DEF_CRS)
+ with self.assertRaises(ValueError):
+ # Different crs
+ centr_bis.to_crs(ALT_CRS).append(centr)
+ centr_bis.append(centr)
+ self.assertAlmostEqual(centr_bis.lat[0], 1)
+ self.assertAlmostEqual(centr_bis.lat[1], 2)
+ self.assertAlmostEqual(centr_bis.lat[2], 3)
+ self.assertAlmostEqual(centr_bis.lon[0], 4)
+ self.assertAlmostEqual(centr_bis.lon[1], 5)
+ self.assertAlmostEqual(centr_bis.lon[2], 6)
+ self.assertTrue(np.array_equal(centr_bis.lat[3:], centr.lat))
+ self.assertTrue(np.array_equal(centr_bis.lon[3:], centr.lon))
+
+ def test_append(self):
+ lat2,lon2 = np.array([6,7,8,9,10]),np.array([6,7,8,9,10])
+ newcentr = Centroids(lat=lat2,lon=lon2)
+ newcentr.append(self.centr)
+ self.assertTrue(newcentr.size == len(self.centr.lon)+len(lon2))
+ np.testing.assert_array_equal(newcentr.lon,np.concatenate([lon2,self.centr.lon]))
+ np.testing.assert_array_equal(newcentr.lat,np.concatenate([lat2,self.centr.lat]))
+
+ def test_append_dif_crs(self):
+ lat2,lon2 = np.array([0,0,1,2,3,4,5]),np.array([0,0,1,2,3,4,5])
+ centr2 = Centroids(lat=lat2,lon=lon2,crs='epsg:3857')
+
+ # appending differing crs is not provided/possible
+ with self.assertRaises(ValueError):
+ self.centr.append(centr2)
+
+ def test_remove_duplicate_pass(self):
+ """Test remove_duplicate_points"""
+ centr = Centroids(
+ lat=np.hstack([VEC_LAT, VEC_LAT]),
+ lon=np.hstack([VEC_LON, VEC_LON]),
+ crs=TEST_CRS,
+ )
+ self.assertTrue(centr.gdf.shape[0] == 2 * self.centr.gdf.shape[0])
+ rem_centr = Centroids.remove_duplicate_points(centr)
+ self.assertIsInstance(rem_centr,Centroids)
+ self.assertTrue(self.centr == rem_centr)
+
+
+ def test_remove_duplicates_dif_on_land(self):
+ ### We currently expect that only the geometry of the gdf defines duplicates.
+ ### If one geometry is duplicated with differences in other attributes e.g. on_land
+ ### they get removed nevertheless. Only the first occurrence will be part of the new object
+ ### this test is only here to guarantee this behaviour
+ lat, lon = np.array([0,0,1,2,3,4,5]),np.array([0,0,1,2,3,4,5])
+ centr = Centroids(lat=lat,lon=lon,on_land=[True]+[False]*6)
+ centr_subset = centr.remove_duplicate_points()
+ # new object created
+ self.assertFalse(centr == centr_subset)
+ self.assertIsNot(centr,centr_subset)
+ # duplicates removed
+ self.assertTrue(centr_subset.size == len(lat)-1)
+ self.assertTrue(np.all(centr_subset.shape == (len(lat)-1,len(lon)-1)))
+ np.testing.assert_array_equal(centr_subset.lon,np.unique(lon))
+ np.testing.assert_array_equal(centr_subset.lat,np.unique(lat))
+ # only first on_land (True) is selected
+ self.assertTrue(centr_subset.on_land[0])
def test_union(self):
lat, lon = np.array([0, 1]), np.array([0, -1])
@@ -142,46 +794,208 @@ def test_union(self):
cent3 = Centroids(lat=lat3,lon=lon3)
cent = cent1.union(cent2)
- np.testing.assert_array_equal(cent.lat, [0, 1, 2, 3])
- np.testing.assert_array_equal(cent.lon, [0, -1, -2, 3])
- np.testing.assert_array_equal(cent.on_land, [True, True, False, False])
+ np.testing.assert_array_equal(cent.lat, np.concatenate([lat,lat2]))
+ np.testing.assert_array_equal(cent.lon, np.concatenate([lon,lon2]))
+ np.testing.assert_array_equal(cent.on_land, np.concatenate([on_land,on_land2]))
cent = cent1.union(cent1, cent2)
- np.testing.assert_array_equal(cent.lat, [0, 1, 2, 3])
- np.testing.assert_array_equal(cent.lon, [0, -1, -2, 3])
- np.testing.assert_array_equal(cent.on_land, [True, True, False, False])
+ np.testing.assert_array_equal(cent.lat, np.concatenate([lat,lat2]))
+ np.testing.assert_array_equal(cent.lon, np.concatenate([lon,lon2]))
+ np.testing.assert_array_equal(cent.on_land, np.concatenate([on_land,on_land2]))
cent = Centroids.union(cent1)
- np.testing.assert_array_equal(cent.lat, [0, 1])
- np.testing.assert_array_equal(cent.lon, [0, -1])
- np.testing.assert_array_equal(cent.on_land, [True, True])
+ np.testing.assert_array_equal(cent.lat, cent1.lat)
+ np.testing.assert_array_equal(cent.lon, cent1.lon)
+ np.testing.assert_array_equal(cent.on_land, cent1.on_land)
cent = cent1.union(cent1)
- np.testing.assert_array_equal(cent.lat, [0, 1])
- np.testing.assert_array_equal(cent.lon, [0, -1])
- np.testing.assert_array_equal(cent.on_land, [True, True])
+ np.testing.assert_array_equal(cent.lat, cent1.lat)
+ np.testing.assert_array_equal(cent.lon, cent1.lon)
+ np.testing.assert_array_equal(cent.on_land, cent1.on_land)
+ # if attributes are not part in one of the centroid objects it will be added as None in the union
cent = Centroids.union(cent1, cent2, cent3)
- np.testing.assert_array_equal(cent.lat, [0, 1, 2, 3, -1, -2])
- np.testing.assert_array_equal(cent.lon, [0, -1, -2, 3, 1, 2])
-
-
- def test_union_meta(self):
- cent1 = Centroids.from_pnt_bounds((-1, -1, 0, 0), res=1)
- cent2 = Centroids.from_pnt_bounds((0, 0, 1, 1), res=1)
- cent3 = Centroids.from_lat_lon(np.array([1]), np.array([1]))
-
- cent = cent1.union(cent2)
- np.testing.assert_array_equal(cent.lat, [0, 0, -1, -1, 1, 1, 0])
- np.testing.assert_array_equal(cent.lon, [-1, 0, -1, 0, 0, 1, 1])
+ np.testing.assert_array_equal(cent.lat, np.concatenate([lat,lat2,lat3]))
+ np.testing.assert_array_equal(cent.lon, np.concatenate([lon,lon2,lon3]))
+ np.testing.assert_array_equal(cent.on_land, np.concatenate([on_land,on_land2,[None,None]]))
+
+ def test_select_pass(self):
+ """Test Centroids.select method"""
+ region_id = np.zeros(VEC_LAT.size)
+ region_id[[2, 4]] = 10
+ centr = Centroids(lat=VEC_LAT, lon=VEC_LON, region_id=region_id)
+
+ fil_centr = centr.select(reg_id=10)
+ self.assertEqual(fil_centr.size, 2)
+ self.assertEqual(fil_centr.lat[0], VEC_LAT[2])
+ self.assertEqual(fil_centr.lat[1], VEC_LAT[4])
+ self.assertEqual(fil_centr.lon[0], VEC_LON[2])
+ self.assertEqual(fil_centr.lon[1], VEC_LON[4])
+ self.assertTrue(np.array_equal(fil_centr.region_id, np.ones(2) * 10))
+
+ def test_select_extent_pass(self):
+ """Test select extent"""
+ centr = Centroids(
+ lat=np.array([-5, -3, 0, 3, 5]),
+ lon=np.array([-180, -175, -170, 170, 175]),
+ region_id=np.zeros(5),
+ )
+ ext_centr = centr.select(extent=[-175, -170, -5, 5])
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-175, -170]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-3, 0]))
+
+ # Cross antimeridian, version 1
+ ext_centr = centr.select(extent=[170, -175, -5, 5])
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
+
+ # Cross antimeridian, version 2
+ ext_centr = centr.select(extent=[170, 185, -5, 5])
+ np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
+ np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
+
+ def test_get_meta(self):
+ """
+ Test that the `get_meta` method correctly generates metadata
+ for a raster with a specified resolution.
+ """
+ # Create centroids with specified resolution
+ lon = np.array([-10.0, -20.0, -30.0])
+ lat = np.array([10.0, 20.0, 30.0])
+ centroids = Centroids(lat=lat, lon=lon, crs=DEF_CRS)
+
+ # Get metadata
+ meta = centroids.get_meta()
+
+ # Check metadata
+ expected_meta = dict(
+ crs=DEF_CRS,
+ height=3,
+ width=3,
+ transform=Affine(
+ 10, 0, -35,
+ 0, -10, 35,
+ ),
+ )
+ self.assertEqual(meta['height'], expected_meta['height'])
+ self.assertEqual(meta['width'], expected_meta['width'])
+ self.assertTrue(u_coord.equal_crs(meta['crs'], expected_meta['crs']))
+ self.assertTrue(meta['transform'].almost_equals(expected_meta['transform']))
+
+ def test_get_closest_point(self):
+ """Test get_closest_point"""
+ for n, (lat, lon) in enumerate(LATLON):
+ x, y, idx = self.centr.get_closest_point(lon * 0.99, lat * 1.01)
+ self.assertAlmostEqual(x, lon)
+ self.assertAlmostEqual(y, lat)
+ self.assertEqual(idx, n)
+ self.assertEqual(self.centr.lon[n], x)
+ self.assertEqual(self.centr.lat[n], y)
+
+ def test_get_closest_point(self):
+ """Test get_closest_point"""
+ for y_sign in [1, -1]:
+ meta = {
+ 'width': 10,
+ 'height': 20,
+ 'transform': rasterio.Affine(0.5, 0, 0.1, 0, y_sign * 0.6, y_sign * (-0.3)),
+ 'crs': DEF_CRS,
+ }
+ centr_ras = Centroids.from_meta(meta=meta)
+
+ test_data = np.array([
+ [0.4, 0.1, 0.35, 0.0, 0],
+ [-0.1, 0.2, 0.35, 0.0, 0],
+ [2.2, 0.1, 2.35, 0.0, 4],
+ [1.4, 2.5, 1.35, 2.4, 42],
+ [5.5, -0.1, 4.85, 0.0, 9],
+ ])
+ test_data[:,[1,3]] *= y_sign
+ for x_in, y_in, x_out, y_out, idx_out in test_data:
+ x, y, idx = centr_ras.get_closest_point(x_in, y_in)
+ self.assertEqual(x, x_out)
+ self.assertEqual(y, y_out)
+ self.assertEqual(idx, idx_out)
+ self.assertEqual(centr_ras.lon[idx], x)
+ self.assertEqual(centr_ras.lat[idx], y)
+
+ centr_ras = Centroids(lat=np.array([0, 0.2, 0.7]), lon=np.array([-0.4, 0.2, 1.1]))
+ x, y, idx = centr_ras.get_closest_point(0.1, 0.0)
+ self.assertEqual(x, 0.2)
+ self.assertEqual(y, 0.2)
+ self.assertEqual(idx, 1)
+
+ def test_dist_coast_pass(self):
+ """Test set_dist_coast"""
+ dist_coast = self.centr.get_dist_coast()
+ # Just checking that the output doesnt change over time.
+ REF_VALUES = np.array([
+ 1605.243, 603.261, 26112.239, 2228.629, 7207.817,
+ 156271.372, 661.114, 158184.4,
+ ])
+ #
+ self.assertIsInstance(dist_coast,np.ndarray)
+ np.testing.assert_array_almost_equal(dist_coast, REF_VALUES, decimal=3)
+
+ def test_dist_coast_pass_raster(self):
+ """Test set_region_id"""
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
+ dist_coast = centr_ras.get_dist_coast()
+ self.assertLess(abs(dist_coast[0] - 117000), 1000)
+ self.assertLess(abs(dist_coast[-1] - 104000), 1000)
+
+ def test_area_pass(self):
+ """Test set_area"""
+ ulx, xres, lrx = 60, 1, 90
+ uly, yres, lry = 0, 1, 20
+ xx, yy = np.meshgrid(np.arange(ulx + xres / 2, lrx, xres),
+ np.arange(uly + yres / 2, lry, yres))
+ vec_data = gpd.GeoDataFrame({
+ 'geometry': [Point(xflat, yflat) for xflat, yflat in zip(xx.flatten(), yy.flatten())],
+ 'lon': xx.flatten(),
+ 'lat': yy.flatten(),
+ }, crs={'proj': 'cea'})
+ centr = Centroids.from_geodataframe(vec_data)
+ area_pixel = centr.get_area_pixel()
+ self.assertTrue(np.allclose(area_pixel, np.ones(centr.size)))
+
+ def test_area_pass_raster(self):
+ """Test set_area"""
+ window_size = (0, 0, 2, 3)
+ centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(*window_size))
+ area_pixel = centr_ras.get_area_pixel()
+
+ # Result in the crs of the test file (ESPG:4326)
+ # This is a wrong result as it should be projected to CEA (for correct area)
+ res = 0.009000000000000341
+ self.assertFalse(
+ np.allclose(area_pixel, np.ones(window_size[2] * window_size[3]) * res**2)
+ )
- cent = cent3.union(cent1)
- np.testing.assert_array_equal(cent.lat, [1, 0, 0, -1, -1])
- np.testing.assert_array_equal(cent.lon, [1, -1, 0, -1, 0])
+ # Correct result in CEA results in unequal pixel area
+ test_area = np.array([
+ 981010.32497514, 981010.3249724 , 981037.92674855,
+ 981037.92674582, 981065.50487659, 981065.50487385,
+ ])
+ np.testing.assert_allclose(area_pixel, test_area)
+
+ def test_equal_pass(self):
+ """Test equal"""
+ centr_list = [
+ Centroids(lat=VEC_LAT, lon=VEC_LON, crs=DEF_CRS),
+ Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS),
+ Centroids(lat=VEC_LAT + 1, lon=VEC_LON + 1)
+ ]
+ for centr1, centr2 in itertools.combinations(centr_list, 2):
+ self.assertFalse(centr2 == centr1)
+ self.assertFalse(centr1 == centr2)
+ self.assertTrue(centr1 == centr1)
+ self.assertTrue(centr2 == centr2)
# Execute Tests
if __name__ == "__main__":
- TESTS = unittest.TestLoader().loadTestsFromTestCase(TestCentroidsReader)
+ TESTS = unittest.TestLoader().loadTestsFromTestCase(TestCentroidsData)
+ TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCentroidsReaderWriter))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCentroidsMethods))
unittest.TextTestRunner(verbosity=2).run(TESTS)
diff --git a/climada/hazard/centroids/test/test_vec_ras.py b/climada/hazard/centroids/test/test_vec_ras.py
deleted file mode 100644
index a4ec17f7e7..0000000000
--- a/climada/hazard/centroids/test/test_vec_ras.py
+++ /dev/null
@@ -1,726 +0,0 @@
-"""
-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 .
-
----
-
-Test CentroidsVector and CentroidsRaster classes.
-"""
-import unittest
-
-from cartopy.io import shapereader
-import geopandas as gpd
-import numpy as np
-from pyproj.crs import CRS
-import rasterio
-from rasterio.windows import Window
-from shapely.geometry.point import Point
-from shapely.geometry.polygon import Polygon
-
-from climada import CONFIG
-from climada.hazard.centroids.centr import Centroids
-from climada.util.constants import HAZ_DEMO_FL, DEF_CRS
-import climada.util.coordinates as u_coord
-
-DATA_DIR = CONFIG.hazard.test_data.dir()
-
-VEC_LON = np.array([
- -59.6250000000000, -59.6250000000000, -59.6250000000000, -59.5416666666667, -59.5416666666667,
- -59.4583333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333,
- -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333,
- -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333,
- -60.2083333333333, -60.1250000000000, -60.1250000000000, -60.1250000000000, -60.1250000000000,
- -60.1250000000000, -60.1250000000000, -60.1250000000000, -60.1250000000000, -60.1250000000000,
- -60.1250000000000, -60.1250000000000, -60.1250000000000, -60.1250000000000, -60.1250000000000,
- -60.1250000000000, -60.1250000000000, -60.0416666666667, -60.0416666666667, -60.0416666666667,
- -60.0416666666667, -60.0416666666667, -60.0416666666667, -60.0416666666667, -60.0416666666667,
- -60.0416666666667, -60.0416666666667, -60.0416666666667, -60.0416666666667, -60.0416666666667,
- -60.0416666666667, -60.0416666666667, -60.0416666666667, -59.9583333333333, -59.9583333333333,
- -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.9583333333333,
- -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.9583333333333,
- -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.9583333333333, -59.8750000000000,
- -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000,
- -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000,
- -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000, -59.8750000000000,
- -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667,
- -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667,
- -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667,
- -59.7916666666667, -59.7083333333333, -59.7083333333333, -59.7083333333333, -59.7083333333333,
- -59.7083333333333, -59.7083333333333, -59.7083333333333, -59.7083333333333, -59.7083333333333,
- -59.7083333333333, -59.7083333333333, -59.7083333333333, -59.7083333333333, -59.7083333333333,
- -59.7083333333333, -59.7083333333333, -59.6250000000000, -59.6250000000000, -59.6250000000000,
- -59.6250000000000, -59.6250000000000, -59.6250000000000, -59.6250000000000, -59.6250000000000,
- -59.6250000000000, -59.6250000000000, -59.6250000000000, -59.6250000000000, -59.6250000000000,
- -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.5416666666667,
- -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.5416666666667,
- -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.5416666666667, -59.4583333333333,
- -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.4583333333333,
- -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.4583333333333,
- -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.4583333333333, -59.3750000000000,
- -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000,
- -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000,
- -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000, -59.3750000000000,
- -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667,
- -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667,
- -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667, -59.2916666666667,
- -59.2916666666667, -59.2083333333333, -59.2083333333333, -59.2083333333333, -59.2083333333333,
- -59.2083333333333, -59.2083333333333, -59.2083333333333, -59.2083333333333, -59.2083333333333,
- -59.2083333333333, -59.2083333333333, -59.2083333333333, -59.2083333333333, -59.2083333333333,
- -59.2083333333333, -59.2083333333333, -59.1250000000000, -59.1250000000000, -59.1250000000000,
- -59.1250000000000, -59.1250000000000, -59.1250000000000, -59.1250000000000, -59.1250000000000,
- -59.1250000000000, -59.1250000000000, -59.1250000000000, -59.1250000000000, -59.1250000000000,
- -59.1250000000000, -59.1250000000000, -59.1250000000000, -59.0416666666667, -59.0416666666667,
- -59.0416666666667, -59.0416666666667, -59.0416666666667, -59.0416666666667, -59.0416666666667,
- -59.0416666666667, -59.0416666666667, -59.0416666666667, -59.0416666666667, -59.0416666666667,
- -59.0416666666667, -59.0416666666667, -59.0416666666667, -58.9583333333333, -58.9583333333333,
- -58.9583333333333, -58.9583333333333, -58.9583333333333, -58.9583333333333, -58.9583333333333,
- -58.9583333333333, -58.9583333333333, -58.9583333333333, -58.9583333333333, -58.9583333333333,
- -58.9583333333333, -58.9583333333333, -61.0416666666667, -61.0416666666667, -61.0416666666667,
- -61.0416666666667, -61.0416666666667, -61.0416666666667, -61.0416666666667, -60.6250000000000,
- -60.6250000000000, -60.6250000000000, -60.6250000000000, -60.6250000000000, -60.6250000000000,
- -60.6250000000000, -60.2083333333333, -60.2083333333333, -60.2083333333333, -60.2083333333333,
- -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.7916666666667, -59.3750000000000,
- -59.3750000000000, -59.3750000000000, -59.3750000000000, -58.9583333333333, -58.9583333333333,
- -58.9583333333333, -58.9583333333333, -58.5416666666667, -58.5416666666667, -58.5416666666667,
- -58.5416666666667, -58.5416666666667, -58.5416666666667, -58.5416666666667, -58.1250000000000,
- -58.1250000000000, -58.1250000000000, -58.1250000000000, -58.1250000000000, -58.1250000000000,
- -58.1250000000000,
-])
-
-VEC_LAT = np.array([
- 13.125, 13.20833333, 13.29166667, 13.125, 13.20833333, 13.125, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667, 13.375,
- 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333,
- 12.79166667, 12.875, 12.95833333, 13.04166667, 13.375, 13.45833333, 13.54166667, 13.625,
- 13.70833333, 13.79166667, 12.54166667, 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333,
- 13.04166667, 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667,
- 12.54166667, 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667,
- 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667,
- 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667,
- 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667,
- 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 13.79166667, 12.54166667,
- 12.625, 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333,
- 13.29166667, 13.375, 13.45833333, 13.54166667, 13.625, 13.70833333, 12.54166667, 12.625,
- 12.70833333, 12.79166667, 12.875, 12.95833333, 13.04166667, 13.125, 13.20833333, 13.29166667,
- 13.375, 13.45833333, 13.54166667, 13.625, 11.875, 12.29166667, 12.70833333, 13.125,
- 13.54166667, 13.95833333, 14.375, 11.875, 12.29166667, 12.70833333, 13.125, 13.54166667,
- 13.95833333, 14.375, 11.875, 12.29166667, 13.95833333, 14.375, 11.875, 12.29166667,
- 13.95833333, 14.375, 11.875, 12.29166667, 13.95833333, 14.375, 11.875, 12.29166667,
- 13.95833333, 14.375, 11.875, 12.29166667, 12.70833333, 13.125, 13.54166667, 13.95833333,
- 14.375, 11.875, 12.29166667, 12.70833333, 13.125, 13.54166667, 13.95833333, 14.375
-])
-
-
-def data_vector():
- vec_data = gpd.GeoDataFrame({
- 'geometry': [Point(lon, lat) for lon, lat in zip(VEC_LON, VEC_LAT)],
- 'lon': VEC_LON,
- 'lat': VEC_LAT
- }, crs='epsg:32632')
- return vec_data.lat.values, vec_data.lon.values, vec_data.geometry
-
-class TestVector(unittest.TestCase):
- """Test CentroidsVector class"""
-
- def test_from_lat_lon_pass(self):
- """Test from_lat_lon"""
- centr = Centroids.from_lat_lon(VEC_LAT, VEC_LON)
- self.assertTrue(np.allclose(centr.lat, VEC_LAT))
- self.assertTrue(np.allclose(centr.lon, VEC_LON))
- self.assertTrue(u_coord.equal_crs(centr.crs, DEF_CRS))
- self.assertTrue(u_coord.equal_crs(centr.geometry.crs, DEF_CRS))
- self.assertEqual(centr.geometry.size, 0)
-
- centr.set_area_pixel()
- self.assertEqual(centr.geometry.size, centr.lat.size)
-
- def test_ne_crs_geom_pass(self):
- """Test _ne_crs_geom"""
- lat, lon, geometry = data_vector()
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- xy_vec = centr._ne_crs_geom()
- lon, lat = xy_vec.geometry[:].x.values, xy_vec.geometry[:].y.values
- self.assertAlmostEqual(4.51072194, lon[0])
- self.assertAlmostEqual(0.00011838, lat[0])
- self.assertAlmostEqual(4.5107354, lon[-1])
- self.assertAlmostEqual(0.0001297, lat[-1])
-
- def test_dist_coast_pass(self):
- """Test set_dist_coast"""
- lat, lon, geometry = data_vector()
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr.geometry.crs = 'epsg:4326'
- centr.set_dist_coast()
- self.assertAlmostEqual(2594.2070842031694, centr.dist_coast[1])
- self.assertAlmostEqual(166295.87602398323, centr.dist_coast[-2])
-
- def test_region_id_pass(self):
- """Test set_region_id"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr.set_region_id()
- self.assertEqual(np.count_nonzero(centr.region_id), 6)
- self.assertEqual(centr.region_id[0], 52) # 052 for barbados
-
- def test_on_land(self):
- """Test set_on_land"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr.set_on_land()
- centr.set_region_id()
- centr.region_id[centr.region_id > 0] = 1
- self.assertTrue(np.array_equal(centr.on_land.astype(int),
- centr.region_id))
-
- def test_remove_duplicate_pass(self):
- """Test remove_duplicate_points"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- # create duplicates manually:
- centr.geometry.values[100] = centr.geometry.values[101]
- centr.geometry.values[120] = centr.geometry.values[101]
- centr.geometry.values[5] = Point([-59.7, 12.5])
- centr.geometry.values[133] = Point([-59.7, 12.5])
- centr.geometry.values[121] = Point([-59.7, 12.5])
- centr.lon = centr.geometry.apply(lambda pt: pt.x).values
- centr.lat = centr.geometry.apply(lambda pt: pt.y).values
- self.assertEqual(centr.size, 296)
- rem_centr = centr.remove_duplicate_points()
- self.assertEqual(centr.size, 296)
- self.assertEqual(rem_centr.size, 292)
- rem2_centr = rem_centr.remove_duplicate_points()
- self.assertEqual(rem_centr.size, 292)
- self.assertEqual(rem2_centr.size, 292)
-
- def test_area_pass(self):
- """Test set_area"""
- ulx, xres, lrx = 60, 1, 90
- uly, yres, lry = 0, 1, 20
- xx, yy = np.meshgrid(np.arange(ulx + xres / 2, lrx, xres),
- np.arange(uly + yres / 2, lry, yres))
- vec_data = gpd.GeoDataFrame({
- 'geometry': [Point(xflat, yflat) for xflat, yflat in zip(xx.flatten(), yy.flatten())],
- 'lon': xx.flatten(),
- 'lat': yy.flatten(),
- }, crs={'proj': 'cea'})
-
- centr = Centroids.from_lat_lon(vec_data.lat.values, vec_data.lon.values)
- centr.geometry = vec_data.geometry
- centr.set_area_pixel()
- self.assertTrue(np.allclose(centr.area_pixel, np.ones(centr.size)))
-
- def test_size_pass(self):
- """Test size property"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- self.assertEqual(centr.size, 296)
-
- def test_get_closest_point(self):
- """Test get_closest_point"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- x, y, idx = centr.get_closest_point(-58.13, 14.38)
- self.assertAlmostEqual(x, -58.125)
- self.assertAlmostEqual(y, 14.375)
- self.assertEqual(idx, 295)
- self.assertEqual(centr.lon[idx], x)
- self.assertEqual(centr.lat[idx], y)
-
- def test_set_lat_lon_to_meta_pass(self):
- """Test set_lat_lon_to_meta"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
-
- centr.set_lat_lon_to_meta()
- self.assertTrue(u_coord.equal_crs(centr.meta['crs'], 'epsg:4326'))
- self.assertEqual(centr.meta['width'], 36)
- self.assertEqual(centr.meta['height'], 31)
- self.assertEqual(centr.meta['transform'][1], 0.0)
- self.assertEqual(centr.meta['transform'][3], 0.0)
- self.assertAlmostEqual(centr.meta['transform'][0], 0.08333333)
- self.assertAlmostEqual(centr.meta['transform'][2], -61.08333333)
- self.assertAlmostEqual(centr.meta['transform'][4], 0.08333333)
- self.assertAlmostEqual(centr.meta['transform'][5], 11.83333333)
-
- def test_get_pixel_polygons_pass(self):
- """Test calc_pixels_polygons"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- poly = centr.calc_pixels_polygons()
- self.assertIsInstance(poly[0], Polygon)
- self.assertTrue(np.allclose(poly.centroid[:].y.values, centr.lat))
- self.assertTrue(np.allclose(poly.centroid[:].x.values, centr.lon))
-
- def test_area_approx(self):
- """Test set_area_approx"""
- lat, lon, geometry = data_vector()
- geometry.crs = 'epsg:4326'
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- with self.assertRaises(ValueError):
- centr.set_area_approx()
-
- def test_append_pass(self):
- """Append points"""
- lat, lon, geometry = data_vector()
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr_bis = Centroids.from_lat_lon(np.array([1, 2, 3]), np.array([4, 5, 6]))
- with self.assertRaises(ValueError):
- centr_bis.append(centr)
- centr.geometry.crs = 'epsg:4326'
- centr_bis.append(centr)
- self.assertAlmostEqual(centr_bis.lat[0], 1)
- self.assertAlmostEqual(centr_bis.lat[1], 2)
- self.assertAlmostEqual(centr_bis.lat[2], 3)
- self.assertAlmostEqual(centr_bis.lon[0], 4)
- self.assertAlmostEqual(centr_bis.lon[1], 5)
- self.assertAlmostEqual(centr_bis.lon[2], 6)
- self.assertTrue(np.array_equal(centr_bis.lat[3:], centr.lat))
- self.assertTrue(np.array_equal(centr_bis.lon[3:], centr.lon))
-
- def test_equal_pass(self):
- """Test equal"""
- lat, lon, geometry = data_vector()
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr_bis = Centroids.from_lat_lon(np.array([1, 2, 3]), np.array([4, 5, 6]))
- self.assertFalse(centr.equal(centr_bis))
- self.assertFalse(centr_bis.equal(centr))
- self.assertTrue(centr_bis.equal(centr_bis))
- self.assertTrue(centr.equal(centr))
-
-
-class TestRaster(unittest.TestCase):
- """Test CentroidsRaster class"""
-
- def test_from_pix_bounds_pass(self):
- """Test from_pix_bounds"""
- xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon = 10, 5, -0.5, 0.2, 20, 25
- centr = Centroids.from_pix_bounds(xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon)
- self.assertTrue(u_coord.equal_crs(centr.meta['crs'], DEF_CRS))
- self.assertEqual(centr.meta['width'], n_lon)
- self.assertEqual(centr.meta['height'], n_lat)
- self.assertAlmostEqual(centr.meta['transform'][0], d_lon)
- self.assertAlmostEqual(centr.meta['transform'][1], 0.0)
- self.assertAlmostEqual(centr.meta['transform'][2], xo_lon)
- self.assertAlmostEqual(centr.meta['transform'][3], 0.0)
- self.assertAlmostEqual(centr.meta['transform'][4], d_lat)
- self.assertAlmostEqual(centr.meta['transform'][5], xf_lat)
- self.assertTrue('lat' in centr.__dict__.keys())
- self.assertTrue('lon' in centr.__dict__.keys())
-
- def test_from_pnt_bounds_pass(self):
- """Test from_pnt_bounds"""
- left, bottom, right, top = 5, 0, 10, 10
- centr = Centroids.from_pnt_bounds((left, bottom, right, top), 0.2)
- self.assertTrue(u_coord.equal_crs(centr.meta['crs'], DEF_CRS))
- self.assertEqual(centr.meta['width'], 26)
- self.assertEqual(centr.meta['height'], 51)
- self.assertAlmostEqual(centr.meta['transform'][0], 0.2)
- self.assertAlmostEqual(centr.meta['transform'][1], 0.0)
- self.assertAlmostEqual(centr.meta['transform'][2], 5 - 0.2 / 2)
- self.assertAlmostEqual(centr.meta['transform'][3], 0.0)
- self.assertAlmostEqual(centr.meta['transform'][4], -0.2)
- self.assertAlmostEqual(centr.meta['transform'][5], 10 + 0.2 / 2)
- self.assertTrue('lat' in centr.__dict__.keys())
- self.assertTrue('lon' in centr.__dict__.keys())
-
- def test_read_all_pass(self):
- """Test centr_ras data"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- self.assertAlmostEqual(centr_ras.meta['crs'], DEF_CRS)
- self.assertAlmostEqual(centr_ras.meta['transform'].c, -69.33714959699981)
- self.assertAlmostEqual(centr_ras.meta['transform'].a, 0.009000000000000341)
- self.assertAlmostEqual(centr_ras.meta['transform'].b, 0.0)
- self.assertAlmostEqual(centr_ras.meta['transform'].f, 10.42822096697894)
- self.assertAlmostEqual(centr_ras.meta['transform'].d, 0.0)
- self.assertAlmostEqual(centr_ras.meta['transform'].e, -0.009000000000000341)
- self.assertEqual(centr_ras.meta['height'], 60)
- self.assertEqual(centr_ras.meta['width'], 50)
-
- inten_ras = centr_ras.values_from_raster_files([HAZ_DEMO_FL], window=Window(0, 0, 50, 60))
- self.assertEqual(inten_ras.shape, (1, 60 * 50))
-
- def test_ne_crs_geom_pass(self):
- """Test _ne_crs_geom"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.meta['crs'] = 'epsg:32632'
-
- xy_vec = centr_ras._ne_crs_geom()
- x_vec, y_vec = xy_vec.geometry[:].x.values, xy_vec.geometry[:].y.values
- self.assertAlmostEqual(4.51063496489, x_vec[0])
- self.assertAlmostEqual(9.40153761711e-05, y_vec[0])
- self.assertAlmostEqual(4.51063891581, x_vec[-1])
- self.assertAlmostEqual(8.92260922066e-05, y_vec[-1])
-
- def test_region_id_pass(self):
- """Test set_dist_coast"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.set_region_id()
- self.assertEqual(centr_ras.region_id.size, centr_ras.size)
- self.assertTrue(np.array_equal(np.unique(centr_ras.region_id), np.array([862])))
-
- def test_set_geometry_points_pass(self):
- """Test set_geometry_points"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.set_geometry_points()
- x_flat = np.arange(-69.3326495969998, -68.88264959699978, 0.009000000000000341)
- y_flat = np.arange(10.423720966978939, 9.883720966978919, -0.009000000000000341)
- x_grid, y_grid = np.meshgrid(x_flat, y_flat)
- self.assertTrue(np.allclose(x_grid.flatten(), centr_ras.lon))
- self.assertTrue(np.allclose(y_grid.flatten(), centr_ras.lat))
-
- def test_dist_coast_pass(self):
- """Test set_region_id"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.set_dist_coast()
- centr_ras.check()
- self.assertTrue(abs(centr_ras.dist_coast[0] - 117000) < 1000)
- self.assertTrue(abs(centr_ras.dist_coast[-1] - 104000) < 1000)
-
- def test_on_land(self):
- """Test set_on_land"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.set_on_land()
- centr_ras.check()
- self.assertTrue(np.array_equal(centr_ras.on_land, np.ones(60 * 50, bool)))
-
- def test_area_pass(self):
- """Test set_area"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.meta['crs'] = {'proj': 'cea'}
- centr_ras.set_area_pixel()
- centr_ras.check()
- self.assertTrue(
- np.allclose(centr_ras.area_pixel,
- np.ones(60 * 50) * 0.009000000000000341 * 0.009000000000000341))
-
- def test_area_approx(self):
- """Test set_area_approx"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_ras.set_area_approx()
- approx_dim = (centr_ras.meta['transform'][0] * 111 * 1000
- * centr_ras.meta['transform'][0] * 111 * 1000)
- self.assertEqual(centr_ras.area_pixel.size, centr_ras.size)
- self.assertEqual(np.unique(centr_ras.area_pixel).size, 60)
- self.assertTrue(np.array_equal((approx_dim / np.unique(centr_ras.area_pixel)).astype(int),
- np.ones(60)))
-
- def test_size_pass(self):
- """Test size property"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- self.assertEqual(centr_ras.size, 50 * 60)
-
- def test_get_closest_point(self):
- """Test get_closest_point"""
- for y_sign in [1, -1]:
- meta = {
- 'width': 10,
- 'height': 20,
- 'transform': rasterio.Affine(0.5, 0, 0.1, 0, y_sign * 0.6, y_sign * (-0.3)),
- 'crs': DEF_CRS,
- }
- centr_ras = Centroids(meta=meta)
-
- test_data = np.array([
- [0.4, 0.1, 0.35, 0.0, 0],
- [-0.1, 0.2, 0.35, 0.0, 0],
- [2.2, 0.1, 2.35, 0.0, 4],
- [1.4, 2.5, 1.35, 2.4, 42],
- [5.5, -0.1, 4.85, 0.0, 9],
- ])
- test_data[:,[1,3]] *= y_sign
- for x_in, y_in, x_out, y_out, idx_out in test_data:
- x, y, idx = centr_ras.get_closest_point(x_in, y_in)
- self.assertEqual(x, x_out)
- self.assertEqual(y, y_out)
- self.assertEqual(idx, idx_out)
- self.assertEqual(centr_ras.lon[idx], x)
- self.assertEqual(centr_ras.lat[idx], y)
-
- centr_ras = Centroids.from_lat_lon(np.array([0, 0.2, 0.7]), np.array([-0.4, 0.2, 1.1]))
- x, y, idx = centr_ras.get_closest_point(0.1, 0.0)
- self.assertEqual(x, 0.2)
- self.assertEqual(y, 0.2)
- self.assertEqual(idx, 1)
-
- def test_set_meta_to_lat_lon_pass(self):
- """Test set_meta_to_lat_lon by using its inverse set_lat_lon_to_meta"""
- lat, lon, geometry = data_vector()
-
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
-
- centr.set_lat_lon_to_meta()
- meta = centr.meta
- centr.set_meta_to_lat_lon()
- self.assertEqual(centr.meta, meta)
- self.assertAlmostEqual(lat.max(), centr.lat.max(), 6)
- self.assertAlmostEqual(lat.min(), centr.lat.min(), 6)
- self.assertAlmostEqual(lon.max(), centr.lon.max(), 6)
- self.assertAlmostEqual(lon.min(), centr.lon.min(), 6)
- self.assertAlmostEqual(np.diff(centr.lon).max(), meta['transform'][0])
- self.assertAlmostEqual(np.diff(centr.lat).max(), meta['transform'][4])
- self.assertTrue(u_coord.equal_crs(geometry.crs, centr.geometry.crs))
-
- def test_append_equal_pass(self):
- """Append raster"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_bis = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_bis.append(centr_ras)
- self.assertAlmostEqual(centr_bis.meta['crs'], DEF_CRS)
- self.assertAlmostEqual(centr_bis.meta['transform'].c, -69.33714959699981)
- self.assertAlmostEqual(centr_bis.meta['transform'].a, 0.009000000000000341)
- self.assertAlmostEqual(centr_bis.meta['transform'].b, 0.0)
- self.assertAlmostEqual(centr_bis.meta['transform'].f, 10.42822096697894)
- self.assertAlmostEqual(centr_bis.meta['transform'].d, 0.0)
- self.assertAlmostEqual(centr_bis.meta['transform'].e, -0.009000000000000341)
- self.assertEqual(centr_bis.meta['height'], 60)
- self.assertEqual(centr_bis.meta['width'], 50)
-
- def test_equal_pass(self):
- """Test equal"""
- centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
- centr_bis = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(51, 61, 10, 10))
- self.assertFalse(centr_ras.equal(centr_bis))
- self.assertFalse(centr_bis.equal(centr_ras))
- self.assertTrue(centr_ras.equal(centr_ras))
- self.assertTrue(centr_bis.equal(centr_bis))
-
-
-class TestCentroids(unittest.TestCase):
- """Test Centroids class"""
-
- def test_centroids_check_pass(self):
- """Test vector data in Centroids"""
- data_vec = data_vector()
- lat, lon, geometry = data_vec
- centr = Centroids(lat=lat, lon=lon, geometry=geometry)
- centr.check()
-
- self.assertTrue(u_coord.equal_crs(centr.crs, data_vec[2].crs))
- self.assertIsInstance(centr.total_bounds, tuple)
- for i in range(4):
- self.assertEqual(centr.total_bounds[i], data_vec[2].total_bounds[i])
-
- self.assertIsInstance(centr.lat, np.ndarray)
- self.assertIsInstance(centr.lon, np.ndarray)
- self.assertIsInstance(centr.coord, np.ndarray)
- self.assertTrue(np.array_equal(centr.lat, VEC_LAT))
- self.assertTrue(np.array_equal(centr.lon, VEC_LON))
- self.assertTrue(np.array_equal(centr.coord, np.array([VEC_LAT, VEC_LON]).transpose()))
- self.assertEqual(centr.size, VEC_LON.size)
-
- centr.area_pixel = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.area_pixel = np.array([])
- centr.dist_coast = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.dist_coast = np.array([])
- centr.on_land = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.on_land = np.array([])
- centr.region_id = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.region_id = np.array([])
- centr.lat = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.lat = np.array([])
- centr.lon = np.array([1])
- with self.assertRaises(ValueError):
- centr.check()
-
- centr.lon = np.array([])
- centr.geometry = gpd.GeoSeries(Point(0,0))
- with self.assertRaises(ValueError) as raised:
- centr.check()
- self.assertEqual(str(raised.exception), 'Wrong geometry size: 0 != 1.')
-
- meta = {
- 'width': 10,
- 'height': 20,
- 'transform': rasterio.Affine(0.1, 0, 0, 0, 0.1, 0),
- 'crs': DEF_CRS,
- }
- cen = Centroids(meta=meta)
- with self.assertRaises(ValueError):
- cen.check()
-
-class TestReader(unittest.TestCase):
- """Test Centroids setter vector and raster methods"""
- def test_from_vector_file(self):
- """Test from_vector_file and values_from_vector_files"""
- shp_file = shapereader.natural_earth(resolution='110m', category='cultural',
- name='populated_places_simple')
- centr = Centroids.from_vector_file(shp_file)
- inten = centr.values_from_vector_files([shp_file], val_names=['pop_min', 'pop_max'])
-
- self.assertTrue(u_coord.equal_crs(centr.geometry.crs, u_coord.NE_EPSG))
- self.assertEqual(centr.geometry.size, centr.lat.size)
- self.assertTrue(u_coord.equal_crs(centr.geometry.crs, u_coord.NE_EPSG))
- self.assertAlmostEqual(centr.lon[0], 12.453386544971766)
- self.assertAlmostEqual(centr.lon[-1], 114.18306345846304)
- self.assertAlmostEqual(centr.lat[0], 41.903282179960115)
- self.assertAlmostEqual(centr.lat[-1], 22.30692675357551)
-
- self.assertEqual(inten.shape, (2, 243))
- # population min
- self.assertEqual(inten[0, 0], 832)
- self.assertEqual(inten[0, -1], 4551579)
- # population max
- self.assertEqual(inten[1, 0], 832)
- self.assertEqual(inten[1, -1], 7206000)
-
- # Test reading values from file with incompatible geometry
- shp_file = shapereader.natural_earth(resolution='10m', category='cultural',
- name='populated_places_simple')
- with self.assertRaises(ValueError) as cm:
- centr.values_from_vector_files([shp_file], val_names=['pop_min', 'pop_max'])
- self.assertIn(
- "Vector data inconsistent with contained vector", str(cm.exception)
- )
-
- def test_from_raster_file_wrong_fail(self):
- """Test from_raster_file with wrong centroids"""
- centr = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(10, 20, 50, 60))
- self.assertAlmostEqual(centr.meta['crs'], DEF_CRS)
- self.assertAlmostEqual(centr.meta['transform'].c, -69.2471495969998)
- self.assertAlmostEqual(centr.meta['transform'].a, 0.009000000000000341)
- self.assertAlmostEqual(centr.meta['transform'].b, 0.0)
- self.assertAlmostEqual(centr.meta['transform'].f, 10.248220966978932)
- self.assertAlmostEqual(centr.meta['transform'].d, 0.0)
- self.assertAlmostEqual(centr.meta['transform'].e, -0.009000000000000341)
- self.assertEqual(centr.meta['height'], 60)
- self.assertEqual(centr.meta['width'], 50)
-
- inten_ras = centr.values_from_raster_files([HAZ_DEMO_FL], window=Window(10, 20, 50, 60))
- self.assertEqual(inten_ras.shape, (1, 60 * 50))
- self.assertAlmostEqual(inten_ras.reshape((60, 50)).tocsr()[25, 12], 0.056825936)
-
- with self.assertRaises(ValueError):
- centr.values_from_raster_files([HAZ_DEMO_FL], window=Window(10, 20, 52, 60))
-
- def test_write_read_raster_h5(self):
- """Write and read hdf5 format"""
- file_name = str(DATA_DIR.joinpath('test_centr.h5'))
-
- xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon = 10, 5, -0.5, 0.2, 20, 25
- centr = Centroids.from_pix_bounds(xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon)
- centr.write_hdf5(file_name)
-
- centr_read = Centroids.from_hdf5(file_name)
- self.assertTrue(centr_read.meta)
- self.assertFalse(centr_read.lat.size)
- self.assertFalse(centr_read.lon.size)
- self.assertEqual(centr_read.meta['width'], centr.meta['width'])
- self.assertEqual(centr_read.meta['height'], centr.meta['height'])
- self.assertAlmostEqual(centr_read.meta['transform'].a, centr.meta['transform'].a)
- self.assertAlmostEqual(centr_read.meta['transform'].b, centr.meta['transform'].b)
- self.assertAlmostEqual(centr_read.meta['transform'].c, centr.meta['transform'].c)
- self.assertAlmostEqual(centr_read.meta['transform'].d, centr.meta['transform'].d)
- self.assertAlmostEqual(centr_read.meta['transform'].e, centr.meta['transform'].e)
- self.assertAlmostEqual(centr_read.meta['transform'].f, centr.meta['transform'].f)
- self.assertTrue(u_coord.equal_crs(centr_read.meta['crs'], centr.meta['crs']))
-
- def test_write_read_points_h5(self):
- file_name = str(DATA_DIR.joinpath('test_centr.h5'))
-
- centr = Centroids.from_lat_lon(VEC_LAT, VEC_LON)
- centr.write_hdf5(file_name)
-
- centr_read = Centroids.from_hdf5(file_name)
- self.assertFalse(centr_read.meta)
- self.assertTrue(centr_read.lat.size)
- self.assertTrue(centr_read.lon.size)
- self.assertTrue(np.allclose(centr_read.lat, centr.lat))
- self.assertTrue(np.allclose(centr_read.lon, centr.lon))
- self.assertTrue(u_coord.equal_crs(centr_read.crs, centr.crs))
-
-class TestCentroidsFuncs(unittest.TestCase):
- """Test Centroids methods"""
- def test_select_pass(self):
- """Test Centroids.select method"""
- centr = Centroids.from_lat_lon(VEC_LAT, VEC_LON)
-
- centr.region_id = np.zeros(VEC_LAT.size)
- centr.region_id[[100, 200]] = 10
-
- fil_centr = centr.select(10)
- self.assertEqual(fil_centr.size, 2)
- self.assertEqual(fil_centr.lat[0], VEC_LAT[100])
- self.assertEqual(fil_centr.lat[1], VEC_LAT[200])
- self.assertEqual(fil_centr.lon[0], VEC_LON[100])
- self.assertEqual(fil_centr.lon[1], VEC_LON[200])
- self.assertTrue(np.array_equal(fil_centr.region_id, np.ones(2) * 10))
-
- def test_select_extent_pass(self):
- """Test select extent"""
- centr = Centroids.from_lat_lon(
- np.array([-5, -3, 0, 3, 5]), np.array([-180, -175, -170, 170, 175]))
- centr.check()
- centr.region_id = np.zeros(len(centr.lat))
- ext_centr = centr.select(extent=[-175, -170, -5, 5])
- np.testing.assert_array_equal(ext_centr.lon, np.array([-175, -170]))
- np.testing.assert_array_equal(ext_centr.lat, np.array([-3, 0]))
-
- # Cross antimeridian, version 1
- ext_centr = centr.select(extent=[170, -175, -5, 5])
- np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
- np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
-
- # Cross antimeridian, version 2
- ext_centr = centr.select(extent=[170, 185, -5, 5])
- np.testing.assert_array_equal(ext_centr.lon, np.array([-180, -175, 170, 175]))
- np.testing.assert_array_equal(ext_centr.lat, np.array([-5, -3, 3, 5]))
-
-# Execute Tests
-if __name__ == "__main__":
- TESTS = unittest.TestLoader().loadTestsFromTestCase(TestRaster)
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestVector))
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCentroids))
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestReader))
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCentroidsFuncs))
- unittest.TextTestRunner(verbosity=2).run(TESTS)
diff --git a/climada/hazard/storm_europe.py b/climada/hazard/storm_europe.py
index ffa7a9676b..a0c5b77718 100644
--- a/climada/hazard/storm_europe.py
+++ b/climada/hazard/storm_europe.py
@@ -557,9 +557,7 @@ def _centroids_from_nc(file_name):
if create_meshgrid:
lats, lons = np.array([np.repeat(lats, len(lons)),
np.tile(lons, len(lats))])
- cent = Centroids.from_lat_lon(lats, lons)
- cent.set_area_pixel()
- cent.set_on_land()
+ cent = Centroids(lat=lats, lon=lons, on_land='natural_earth')
return cent
@@ -658,26 +656,27 @@ def calc_ssi(self, method='dawkins', intensity=None, on_land=True,
intensity = intensity.multiply(intensity > self.intensity_thres)
cent = self.centroids
+ area_pixel = cent.get_area_pixel()
if sel_cen is not None:
pass
elif on_land is True:
- sel_cen = cent.on_land
+ sel_cen = cent.on_land.astype(bool)
else: # select all centroids
- sel_cen = np.ones_like(cent.area_pixel, dtype=bool)
+ sel_cen = np.ones_like(area_pixel, dtype=bool)
ssi = np.zeros(intensity.shape[0])
if method == 'dawkins':
- area_c = cent.area_pixel / 1000 / 1000 * sel_cen
+ area_c = area_pixel / 1000 / 1000 * sel_cen
for i, inten_i in enumerate(intensity):
- ssi_i = area_c * inten_i.power(3).todense().T
+ ssi_i = inten_i.power(3).dot(area_c)
# matrix crossproduct (row x column vector)
ssi[i] = ssi_i.item(0)
elif method == 'wisc_gust':
for i, inten_i in enumerate(intensity[:, sel_cen]):
- area = np.sum(cent.area_pixel[inten_i.indices]) / 1000 / 1000
+ area = np.sum(area_pixel[inten_i.indices]) / 1000 / 1000
inten_mean = np.mean(inten_i)
ssi[i] = area * np.power(inten_mean, 3)
@@ -774,8 +773,7 @@ def generate_prob_storms(self, reg_id=528, spatial_shift=4, ssi_args=None,
"""
# bool vector selecting the targeted centroids
if reg_id is not None:
- if self.centroids.region_id.size == 0:
- self.centroids.set_region_id()
+ self.centroids.set_region_id(overwrite=False)
if not isinstance(reg_id, list):
reg_id = [reg_id]
sel_cen = np.isin(self.centroids.region_id, reg_id)
diff --git a/climada/hazard/tc_tracks.py b/climada/hazard/tc_tracks.py
index d4e28ca632..89a44b3744 100644
--- a/climada/hazard/tc_tracks.py
+++ b/climada/hazard/tc_tracks.py
@@ -1214,7 +1214,7 @@ def generate_centroids(self, res_deg, buffer_deg):
lat = np.arange(bounds[1] + 0.5 * res_deg, bounds[3], res_deg)
lon = np.arange(bounds[0] + 0.5 * res_deg, bounds[2], res_deg)
lon, lat = [ar.ravel() for ar in np.meshgrid(lon, lat)]
- return Centroids.from_lat_lon(lat, lon)
+ return Centroids(lat=lat, lon=lon)
def plot(self, axis=None, figsize=(9, 13), legend=True, adapt_fontsize=True, **kwargs):
"""Track over earth. Historical events are blue, probabilistic black.
diff --git a/climada/hazard/test/data/fp_centroids-test.xls b/climada/hazard/test/data/fp_centroids-test.xls
index 0df4fae2ba..270daf2f55 100644
Binary files a/climada/hazard/test/data/fp_centroids-test.xls and b/climada/hazard/test/data/fp_centroids-test.xls differ
diff --git a/climada/hazard/test/test_base.py b/climada/hazard/test/test_base.py
index 8b586d6a55..f322412856 100644
--- a/climada/hazard/test/test_base.py
+++ b/climada/hazard/test/test_base.py
@@ -41,10 +41,7 @@
"""
Directory for writing (and subsequent reading) of temporary files created during tests.
"""
-HAZ_TEST_MAT :Path = get_test_file('atl_prob_no_name', file_format='matlab')
-"""
-Hazard test file from Git repository. Fraction is 1. Format: matlab.
-"""
+
HAZ_TEST_TC :Path = get_test_file('test_tc_florida')
"""
Hazard test file from Data API: Hurricanes from 1851 to 2011 over Florida with 100 centroids.
@@ -66,8 +63,8 @@ def dummy_hazard():
"TC",
intensity=intensity,
fraction=fraction,
- centroids=Centroids.from_lat_lon(
- np.array([1, 3, 5]), np.array([2, 4, 6])),
+ centroids=Centroids(
+ lat=np.array([1, 3, 5]), lon=np.array([2, 4, 6])),
event_id=np.array([1, 2, 3, 4]),
event_name=['ev1', 'ev2', 'ev3', 'ev4'],
date=np.array([1, 2, 3, 4]),
@@ -82,8 +79,11 @@ class TestLoader(unittest.TestCase):
def setUp(self):
"""Test fixure: Build a valid hazard"""
- centroids = Centroids.from_lat_lon(np.array([1, 3]), np.array([2, 3]))
- centroids.region_id = np.array([1, 2])
+ centroids = Centroids(
+ lat=np.array([1, 3]),
+ lon=np.array([2, 3]),
+ region_id=np.array([1, 2]),
+ )
self.hazard = Hazard(
"TC",
centroids=centroids,
@@ -114,12 +114,6 @@ def test_init_empty_fraction(self):
np.testing.assert_array_equal(hazard.fraction.shape, hazard.intensity.shape)
self.assertEqual(hazard.fraction.nnz, 0) # No nonzero entries
- def test_check_wrongCentroids_fail(self):
- """Wrong hazard definition"""
- self.hazard.centroids.region_id = np.array([1, 2, 3, 4])
- with self.assertRaises(ValueError):
- self.hazard.check()
-
def test_check_wrongFreq_fail(self):
"""Wrong hazard definition"""
self.hazard.frequency = np.array([1, 2])
@@ -247,7 +241,7 @@ def test_same_events_same(self):
duplicate events, initial events are obtained with 0 intensity and
fraction in new appended centroids."""
haz1 = dummy_hazard()
- centroids = Centroids.from_lat_lon(np.array([7, 9, 11]), np.array([8, 10, 12]))
+ centroids = Centroids(lat=np.array([7, 9, 11]), lon=np.array([8, 10, 12]))
fraction = sparse.csr_matrix([[0.22, 0.32, 0.44],
[0.11, 0.11, 0.11],
[0.32, 0.11, 0.99],
@@ -507,7 +501,7 @@ def test_select_date_invalid_pass(self):
def test_select_reg_id_pass(self):
"""Test select region of centroids."""
haz = dummy_hazard()
- haz.centroids.region_id = np.array([5, 7, 9])
+ haz.centroids.gdf['region_id'] = np.array([5, 7, 9])
sel_haz = haz.select(date=(2, 4), orig=False, reg_id=9)
self.assertTrue(np.array_equal(sel_haz.centroids.coord.squeeze(),
@@ -724,8 +718,8 @@ def test_all_different_extend(self):
haz2 = Hazard('TC',
date=np.ones((4,)),
orig=np.ones((4,)),
- centroids=Centroids.from_lat_lon(
- np.array([7, 9, 11]), np.array([8, 10, 12])),
+ centroids=Centroids(
+ lat=np.array([7, 9, 11]), lon=np.array([8, 10, 12])),
event_id=np.array([5, 6, 7, 8]),
event_name=['ev5', 'ev6', 'ev7', 'ev8'],
frequency=np.array([0.9, 0.75, 0.75, 0.22]),
@@ -775,8 +769,8 @@ def test_same_events_append(self):
[8.33, 4.11, 4.4],
[9.33, 9.22, 1.77]])
haz2 = Hazard('TC',
- centroids=Centroids.from_lat_lon(
- np.array([7, 9, 11]), np.array([8, 10, 12])),
+ centroids=Centroids(
+ lat=np.array([7, 9, 11]), lon=np.array([8, 10, 12])),
event_id=haz1.event_id,
event_name=haz1.event_name.copy(),
frequency=haz1.frequency,
@@ -822,8 +816,8 @@ def test_concat_pass(self):
"""Test concatenate function."""
haz_1 = Hazard("TC",
- centroids=Centroids.from_lat_lon(
- np.array([1, 3, 5]), np.array([2, 4, 6])),
+ centroids=Centroids(
+ lat=np.array([1, 3, 5]), lon=np.array([2, 4, 6])),
event_id=np.array([1]),
event_name=['ev1'],
date=np.array([1]),
@@ -832,11 +826,10 @@ def test_concat_pass(self):
frequency_unit='1/week',
fraction=sparse.csr_matrix([[0.02, 0.03, 0.04]]),
intensity=sparse.csr_matrix([[0.2, 0.3, 0.4]]),
- units='m/s',)
+ units='m/s')
haz_2 = Hazard("TC",
- centroids=Centroids.from_lat_lon(
- np.array([1, 3, 5]), np.array([2, 4, 6])),
+ centroids=Centroids(lat=np.array([1, 3, 5]), lon=np.array([2, 4, 6])),
event_id=np.array([1]),
event_name=['ev2'],
date=np.array([2]),
@@ -845,7 +838,7 @@ def test_concat_pass(self):
frequency_unit='1/week',
fraction=sparse.csr_matrix([[1.02, 1.03, 1.04]]),
intensity=sparse.csr_matrix([[1.2, 1.3, 1.4]]),
- units='m/s',)
+ units='m/s')
haz = Hazard.concat([haz_1, haz_2])
@@ -931,7 +924,7 @@ def test_change_centroids(self):
self.assertTrue(np.array_equal(haz_2.orig, [True]))
"""Test error for projection"""
- lat3, lon3 = np.array([0.5, 3]), np.array([-0.5, 3])
+ lat3, lon3 = np.array([0.5, 3, 1]), np.array([-0.5, 3, 1])
on_land3 = np.array([True, True, False])
cent3 = Centroids(lat=lat3, lon=lon3, on_land=on_land3)
@@ -961,8 +954,6 @@ def test_change_centroids_raster(self):
"""Test with raster centroids"""
cent4 = Centroids.from_pnt_bounds(points_bounds=(-1, 0, 0, 1), res=1)
- cent4.lat, cent1.lon = np.array([0, 1]), np.array([0, -1])
- cent4.on_land = np.array([True, True])
haz_4 = haz_1.change_centroids(cent4)
@@ -1089,63 +1080,6 @@ def test_hazard_pass(self):
self.assertEqual(hazard.haz_type, 'TC')
-class TestReaderMat(unittest.TestCase):
- """Test reader functionality of the ExposuresExcel class"""
-
- def test_hazard_pass(self):
- """Read a hazard matlab file correctly."""
- # Read demo matlab file
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
-
- # Check results
- n_events = 14450
- n_centroids = 100
-
- self.assertEqual(hazard.units, 'm/s')
-
- self.assertEqual(hazard.centroids.coord.shape, (n_centroids, 2))
-
- self.assertEqual(hazard.event_id.dtype, int)
- self.assertEqual(hazard.event_id.shape, (n_events,))
-
- self.assertEqual(hazard.frequency.dtype, float)
- self.assertEqual(hazard.frequency.shape, (n_events,))
-
- self.assertEqual(hazard.frequency_unit, DEF_FREQ_UNIT)
-
- self.assertEqual(hazard.intensity.dtype, float)
- self.assertEqual(hazard.intensity.shape, (n_events, n_centroids))
- self.assertEqual(hazard.intensity[12, 46], 12.071393519949979)
- self.assertEqual(hazard.intensity[13676, 49], 17.228323602220616)
-
- self.assertEqual(hazard.fraction.dtype, float)
- self.assertEqual(hazard.fraction.shape, (n_events, n_centroids))
- self.assertEqual(hazard.fraction[8454, 98], 1)
- self.assertEqual(hazard.fraction[85, 54], 0)
-
- self.assertEqual(len(hazard.event_name), n_events)
- self.assertEqual(hazard.event_name[124], 125)
-
- self.assertEqual(len(hazard.date), n_events)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[0]).year, 1851)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[0]).month, 6)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[0]).day, 25)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[78]).year, 1852)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[78]).month, 9)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[78]).day, 22)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[-1]).year, 2011)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[-1]).month, 11)
- self.assertEqual(dt.datetime.fromordinal(hazard.date[-1]).day, 6)
-
- self.assertTrue(hazard.orig[0])
- self.assertTrue(hazard.orig[11580])
- self.assertTrue(hazard.orig[4940])
- self.assertFalse(hazard.orig[3551])
- self.assertFalse(hazard.orig[10651])
- self.assertFalse(hazard.orig[4818])
-
- self.assertEqual(hazard.haz_type, 'TC')
-
class TestHDF5(unittest.TestCase):
"""Test reader functionality of the ExposuresExcel class"""
@@ -1185,8 +1119,8 @@ def test_reproject_vector_pass(self):
event_name=['1'],
intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])),
fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2),
- centroids=Centroids.from_lat_lon(
- np.array([1, 2, 3]), np.array([1, 2, 3])),)
+ centroids=Centroids(
+ lat=np.array([1, 2, 3]), lon=np.array([1, 2, 3])),)
haz_fl.check()
haz_fl.reproject_vector(dst_crs='epsg:2202')
@@ -1198,69 +1132,6 @@ def test_reproject_vector_pass(self):
self.assertTrue(np.allclose(haz_fl.intensity.toarray(), np.array([0.5, 0.2, 0.1])))
self.assertTrue(np.allclose(haz_fl.fraction.toarray(), np.array([0.5, 0.2, 0.1]) / 2))
- def test_vector_to_raster_pass(self):
- """Test vector_to_raster"""
- haz_fl = Hazard('FL',
- event_id=np.array([1]),
- date=np.array([1]),
- frequency=np.array([1]),
- orig=np.array([1]),
- event_name=['1'],
- intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])),
- fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2),
- centroids=Centroids.from_lat_lon(
- np.array([1, 2, 3]), np.array([1, 2, 3])),)
- haz_fl.check()
-
- haz_fl.vector_to_raster()
- self.assertTrue(u_coord.equal_crs(haz_fl.centroids.meta['crs'], 'epsg:4326'))
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][0], 1.0)
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][1], 0)
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][2], 0.5)
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][3], 0)
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][4], -1.0)
- self.assertAlmostEqual(haz_fl.centroids.meta['transform'][5], 3.5)
- self.assertEqual(haz_fl.centroids.meta['height'], 3)
- self.assertEqual(haz_fl.centroids.meta['width'], 3)
- self.assertEqual(haz_fl.centroids.lat.size, 0)
- self.assertEqual(haz_fl.centroids.lon.size, 0)
- self.assertTrue(haz_fl.intensity.min() >= 0)
- self.assertTrue(haz_fl.intensity.max() <= 0.5)
- self.assertTrue(haz_fl.fraction.min() >= 0)
- self.assertTrue(haz_fl.fraction.max() <= 0.5 / 2)
-
-
-class TestClear(unittest.TestCase):
- """Test clear method"""
-
- def test_clear(self):
- """Clear method clears everything"""
- haz1 = Hazard.from_excel(HAZ_TEMPLATE_XLS, haz_type='TC')
- haz1.units = "m"
- haz1.frequency_unit = "1/m"
- haz1.foo = np.arange(10)
- haz1.clear()
- self.assertEqual(haz1.haz_type, '')
- self.assertEqual(haz1.units, '')
- self.assertEqual(haz1.frequency_unit, DEF_FREQ_UNIT)
- self.assertEqual(haz1.centroids.size, 0)
- self.assertEqual(len(haz1.event_name), 0)
- for attr in vars(haz1).keys():
- if attr not in ['haz_type', 'units', 'event_name', 'pool', 'frequency_unit']:
- self.assertEqual(getattr(haz1, attr).size, 0)
- self.assertIsNone(haz1.pool)
-
- def test_clear_pool(self):
- """Clear method should not clear a process pool"""
- haz1 = Hazard.from_excel(HAZ_TEMPLATE_XLS, haz_type='TC')
- pool = Pool(nodes=2)
- haz1.pool = pool
- haz1.check()
- haz1.clear()
- self.assertEqual(haz1.pool, pool)
- pool.close()
- pool.join()
- pool.clear()
def dummy_step_impf(haz):
from climada.entity import ImpactFunc
@@ -1375,13 +1246,11 @@ def test_get_fraction(self):
TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLoader)
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestHDF5))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestReaderExcel))
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestReaderMat))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRemoveDupl))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestSelect))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestStats))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestYearset))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAppend))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCentroids))
- TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestClear))
TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestImpactFuncs))
unittest.TextTestRunner(verbosity=2).run(TESTS)
diff --git a/climada/hazard/test/test_storm_europe.py b/climada/hazard/test/test_storm_europe.py
index e6927e1852..f919cbaa45 100644
--- a/climada/hazard/test/test_storm_europe.py
+++ b/climada/hazard/test/test_storm_europe.py
@@ -27,9 +27,8 @@
from climada import CONFIG
from climada.hazard.storm_europe import StormEurope, generate_WS_forecast_hazard
-from climada.hazard.centroids.centr import DEF_VAR_EXCEL, Centroids
+from climada.hazard.centroids.centr import Centroids
from climada.util.constants import WS_DEMO_NC
-from climada.util.api_client import Client
DATA_DIR = CONFIG.hazard.test_data.dir()
@@ -64,12 +63,9 @@ def test_read_with_ref(self):
def test_read_with_cent(self):
"""Test from_footprints while passing in a Centroids object"""
- var_names = copy.deepcopy(DEF_VAR_EXCEL)
- var_names['sheet_name'] = 'fp_centroids-test'
- var_names['col_name']['region_id'] = 'iso_n3'
test_centroids = Centroids.from_excel(
- DATA_DIR.joinpath('fp_centroids-test.xls'),
- var_names=var_names
+ file_path=DATA_DIR.joinpath('fp_centroids-test.xls'),
+ sheet_name='fp_centroids-test'
)
storms = StormEurope.from_footprints(WS_DEMO_NC, centroids=test_centroids)
diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py
index 41fc799722..0452c6b78c 100644
--- a/climada/hazard/test/test_trop_cyclone.py
+++ b/climada/hazard/test/test_trop_cyclone.py
@@ -39,6 +39,7 @@
from climada.hazard.centroids.centr import Centroids
import climada.hazard.test as hazard_test
+
DATA_DIR = Path(hazard_test.__file__).parent.joinpath('data')
TEST_TRACK = DATA_DIR.joinpath("trac_brb_test.csv")
@@ -517,7 +518,7 @@ def test_apply_criterion_track(self):
np.allclose(tc.intensity[2, :].toarray(), tc_cc.intensity[2, :].toarray()))
self.assertTrue(np.allclose(tc.frequency, tc_cc.frequency))
- def test_apply_criterion_track(self):
+ def test_apply_criterion_track2(self):
"""Test _apply_criterion function."""
criterion = [{'basin': 'NA', 'category': [1, 2, 3, 4, 5],
'year': 2100, 'change': 1.045, 'variable': 'intensity'}
diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py
index 38eb6d00d9..686a88b214 100644
--- a/climada/hazard/trop_cyclone.py
+++ b/climada/hazard/trop_cyclone.py
@@ -178,7 +178,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 = 'H08',
ignore_distance_to_coast: bool = False,
@@ -223,12 +223,16 @@ def from_tracks(
"ER11" (Emanuel and Rotunno 2011).
Default: "H08".
ignore_distance_to_coast : boolean, optional
- If True, centroids far from coast are not ignored. Default: False.
+ If True, centroids far from coast are not ignored.
+ If False, the centroids' distances to the coast are calculated with the
+ `Centroids.get_dist_coast()` method (unless there is "dist_coast" column in the
+ centroids' GeoDataFrame) and centroids far from coast are ignored.
+ Default: False.
store_windfields : boolean, optional
- If True, the Hazard object gets a list ``windfields`` of sparse matrices. For each track,
- the full velocity vectors at each centroid and track position are stored in a sparse
- matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full ndarray
- of shape (npositions, ncentroids, 2). Default: False.
+ If True, the Hazard object gets a list ``windfields`` of sparse matrices. For each
+ track, the full velocity vectors at each centroid and track position are stored in a
+ sparse matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full
+ ndarray of shape (npositions, ncentroids, 2). Default: False.
metric : str, optional
Specify an approximation method to use for earth distances:
@@ -262,21 +266,18 @@ def from_tracks(
TropCyclone
"""
num_tracks = tracks.size
- if centroids is None:
- centroids = Centroids.from_base_grid(res_as=360, land=False)
-
- 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()
+ if 'dist_coast' not in centroids.gdf.columns:
+ dist_coast = centroids.get_dist_coast()
+ else:
+ dist_coast = centroids.gdf.dist_coast.values
[idx_centr_filter] = (
- (centroids.dist_coast <= max_dist_inland_km * 1000)
+ (dist_coast <= max_dist_inland_km * 1000)
& (np.abs(centroids.lat) <= max_latitude)
).nonzero()
@@ -539,7 +540,8 @@ def from_single_track(
If True, store windfields. Default: False.
metric : str, optional
Specify an approximation method to use for earth distances: "equirect" (faster) or
- "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``.
+ "geosphere" (more accurate). See ``dist_approx`` function in
+ ``climada.util.coordinates``.
Default: "equirect".
intensity_thres : float, optional
Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5
@@ -1430,9 +1432,9 @@ def _x_holland_2010(
mask_centr_close : np.ndarray of shape (nnodes, ncentroids)
Mask indicating for each track node which centroids are within reach of the windfield.
v_n : np.ndarray of shape (nnodes,) or float, optional
- Peripheral wind speeds (in m/s) at radius ``r_n`` outside of radius of maximum winds ``r_max``.
- In absence of a second wind speed measurement, this value defaults to 17 m/s following
- Holland et al. 2010 (at a radius of 300 km).
+ Peripheral wind speeds (in m/s) at radius ``r_n`` outside of radius of maximum winds
+ ``r_max``. In absence of a second wind speed measurement, this value defaults to 17 m/s
+ following Holland et al. 2010 (at a radius of 300 km).
r_n_km : np.ndarray of shape (nnodes,) or float, optional
Radius (in km) where the peripheral wind speed ``v_n`` is measured (or assumed).
In absence of a second wind speed measurement, this value defaults to 300 km following
@@ -1556,8 +1558,8 @@ def _stat_holland_1980(
Parameters
----------
si_track : xr.Dataset
- Output of ``tctrack_to_si`` with "hol_b" (see, e.g., _B_holland_1980) data variable. The data
- variables used by this function are "lat", "cp", "rad", "cen", "env", and "hol_b".
+ Output of ``tctrack_to_si`` with "hol_b" (see, e.g., _B_holland_1980) data variable. The
+ data variables used by this function are "lat", "cp", "rad", "cen", "env", and "hol_b".
d_centr : np.ndarray of shape (nnodes, ncentroids)
Distance (in m) between centroids and track nodes.
mask_centr_close : np.ndarray of shape (nnodes, ncentroids)
@@ -1621,8 +1623,8 @@ def _stat_er_2011(
Parameters
----------
si_track : xr.Dataset
- Output of ``tctrack_to_si``. The data variables used by this function are "lat", "cp", "rad",
- and "vmax".
+ Output of ``tctrack_to_si``. The data variables used by this function are "lat", "cp",
+ "rad", and "vmax".
d_centr : np.ndarray of shape (nnodes, ncentroids)
Distance (in m) between centroids and track nodes.
mask_centr_close : np.ndarray of shape (nnodes, ncentroids)
diff --git a/climada/test/test_api_client.py b/climada/test/test_api_client.py
index 92421f3992..6c6c0e72e1 100644
--- a/climada/test/test_api_client.py
+++ b/climada/test/test_api_client.py
@@ -212,6 +212,10 @@ def test_get_litpop_fail(self):
client.get_litpop(['AUT', 'CHE'])
self.assertIn(" can only query single countries. Download the data for multiple countries individually and concatenate ",
str(cm.exception))
+
+ def test_get_centroids_plot(self):
+ client = Client()
+ client.get_centroids(country='COM').plot()
def test_get_dataset_file(self):
client = Client()
diff --git a/climada/test/test_calibration.py b/climada/test/test_calibration.py
index 96c3353c6b..0986580d48 100644
--- a/climada/test/test_calibration.py
+++ b/climada/test/test_calibration.py
@@ -32,7 +32,7 @@
from climada.test import get_test_file
-HAZ_TEST_MAT = get_test_file('atl_prob_no_name', file_format='matlab')
+HAZ_TEST_TC = get_test_file('test_tc_florida', file_format='hdf5')
DATA_FOLDER = CONFIG.test_data.dir()
@@ -47,7 +47,7 @@ def test_calib_instance(self):
ent.check()
# Read default hazard file
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
# get impact function from set
imp_func = ent.impact_funcs.get_func(hazard.haz_type,
@@ -59,10 +59,10 @@ def test_calib_instance(self):
# create input frame
df_in = pd.DataFrame.from_dict({'v_threshold': [25.7],
'other_param': [2],
- 'hazard': [HAZ_TEST_MAT]})
+ 'hazard': [HAZ_TEST_TC]})
df_in_yearly = pd.DataFrame.from_dict({'v_threshold': [25.7],
'other_param': [2],
- 'hazard': [HAZ_TEST_MAT]})
+ 'hazard': [HAZ_TEST_TC]})
# Compute the impact over the whole exposures
df_out = calib_instance(hazard, ent.exposures, imp_func, df_in)
diff --git a/climada/test/test_hazard.py b/climada/test/test_hazard.py
index 853ce261b3..6166e9ecd3 100644
--- a/climada/test/test_hazard.py
+++ b/climada/test/test_hazard.py
@@ -30,15 +30,18 @@
from climada.hazard.base import Hazard
from climada.hazard.centroids import Centroids
from climada.hazard.storm_europe import StormEurope
-from climada.util.constants import (HAZ_DEMO_FL, WS_DEMO_NC)
+from climada.util.constants import (HAZ_DEMO_FL, WS_DEMO_NC, DEF_CRS)
from climada.util.api_client import Client
from climada.util import coordinates as u_coord
from climada.test import get_test_file
DATA_DIR = CONFIG.test_data.dir()
-# Hazard test file from Git repository. Fraction is 1. Format: matlab.
-HAZ_TEST_MAT :Path = get_test_file('atl_prob_no_name', file_format='matlab')
+HAZ_TEST_TC :Path = get_test_file('test_tc_florida', file_format='hdf5')
+"""
+Hazard test file from Data API: Hurricanes from 1851 to 2011 over Florida with 100 centroids.
+Fraction is empty. Format: HDF5.
+"""
class TestCentroids(unittest.TestCase):
"""Test centroids functionalities"""
@@ -47,18 +50,18 @@ def test_read_write_raster_pass(self):
"""Test write_raster: Hazard from raster data"""
haz_fl = Hazard.from_raster([HAZ_DEMO_FL])
haz_fl.haz_type = 'FL'
- haz_fl.check()
self.assertEqual(haz_fl.intensity.shape, (1, 1032226))
self.assertEqual(haz_fl.intensity.min(), -9999.0)
self.assertAlmostEqual(haz_fl.intensity.max(), 4.662774085998535)
- haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'))
+ haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), variable='intensity')
haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')])
haz_fl.haz_type = 'FL'
self.assertTrue(np.allclose(haz_fl.intensity.toarray(), haz_read.intensity.toarray()))
self.assertEqual(np.unique(np.array(haz_fl.fraction.toarray())).size, 2)
+ DATA_DIR.joinpath('test_write_hazard.tif').unlink()
def test_read_raster_pool_pass(self):
"""Test from_raster constructor with pool"""
@@ -74,51 +77,118 @@ def test_read_raster_pool_pass(self):
pool.join()
def test_read_write_vector_pass(self):
- """Test write_raster: Hazard from vector data"""
- haz_fl = Hazard('FL',
- event_id=np.array([1]),
- date=np.array([1]),
- frequency=np.array([1]),
- orig=np.array([1]),
- event_name=['1'],
- intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])),
- fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2),
- centroids=Centroids.from_lat_lon(
- np.array([1, 2, 3]), np.array([1, 2, 3])),)
- haz_fl.check()
+ """Test write_raster: Rasterize intensity from vector data"""
+ haz_fl = Hazard(
+ 'FL',
+ event_id=np.array([1]),
+ date=np.array([1]),
+ frequency=np.array([1]),
+ orig=np.array([1]),
+ event_name=['1'],
+ intensity=sparse.csr_matrix(np.array([0.11, 0.22, 0.33, 0.31])),
+ fraction=sparse.csr_matrix(np.array([0, 1, 2, 3]) ),
+ centroids=Centroids(
+ lon=np.array([1, 2, 3, 3]), lat=np.array([1, 2, 3, 1]), crs=DEF_CRS
+ )
+ )
- haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'))
+ haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), variable='intensity')
haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')], haz_type='FL')
self.assertEqual(haz_read.intensity.shape, (1, 9))
- self.assertTrue(np.allclose(np.unique(np.array(haz_read.intensity.toarray())),
- np.array([0.0, 0.1, 0.2, 0.5])))
-
- def test_write_fraction_pass(self):
- """Test write_raster with fraction"""
- haz_fl = Hazard('FL',
- event_id=np.array([1]),
- date=np.array([1]),
- frequency=np.array([1]),
- orig=np.array([1]),
- event_name=['1'],
- intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])),
- fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2),
- centroids=Centroids.from_lat_lon(
- np.array([1, 2, 3]), np.array([1, 2, 3])),)
- haz_fl.check()
- haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), intensity=False)
+ output_raster = np.array([
+ [1, 3], [2, 3], [3, 3],
+ [1, 2], [2, 2], [3, 2],
+ [1, 1], [2, 1], [3, 1]
+ ])
+ output_instensity = np.array([
+ 0, 0, 0.33,
+ 0, 0.22, 0,
+ 0.11, 0, 0.31
+ ])
+
+ np.testing.assert_array_equal(
+ haz_read.centroids.lon,
+ output_raster[:, 0]
+ )
+ np.testing.assert_array_equal(
+ haz_read.centroids.lat,
+ output_raster[:, 1]
+ )
+ np.testing.assert_array_almost_equal(
+ haz_read.intensity.toarray().flatten(),
+ output_instensity
+ )
- haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')],
- files_fraction=[DATA_DIR.joinpath('test_write_hazard.tif')],
- haz_type='FL')
- self.assertEqual(haz_read.intensity.shape, (1, 9))
+ DATA_DIR.joinpath('test_write_hazard.tif').unlink()
+
+ def test_read_write_vector_fraction_pass(self):
+ """Test write_raster: Rasterize fraction from vector data"""
+ haz_fl = Hazard(
+ 'FL',
+ event_id=np.array([1]),
+ date=np.array([1]),
+ frequency=np.array([1]),
+ orig=np.array([1]),
+ event_name=['1'],
+ intensity=sparse.csr_matrix(np.array([-0.11, -0.22, -0.33, -0.31])),
+ fraction=sparse.csr_matrix(np.array([0.11, 0.22, 0.33, 0.31])),
+ centroids=Centroids(
+ lon=np.array([1, 2, 3, 3]), lat=np.array([1, 2, 3, 1]), crs=DEF_CRS
+ )
+ )
+
+ intensity_file = DATA_DIR.joinpath('test_write_hazard_intensity.tif')
+ fraction_file = DATA_DIR.joinpath('test_write_hazard_fraction.tif')
+
+ haz_fl.write_raster(fraction_file, variable='fraction')
+ haz_fl.write_raster(intensity_file, variable='intensity')
+
+ haz_read = Hazard.from_raster(
+ [intensity_file], [fraction_file], haz_type='FL'
+ )
self.assertEqual(haz_read.fraction.shape, (1, 9))
- self.assertTrue(np.allclose(np.unique(np.array(haz_read.fraction.toarray())),
- np.array([0.0, 0.05, 0.1, 0.25])))
- self.assertTrue(np.allclose(np.unique(np.array(haz_read.intensity.toarray())),
- np.array([0.0, 0.05, 0.1, 0.25])))
+ self.assertEqual(haz_read.intensity.shape, (1, 9))
+
+
+ output_raster = np.array([
+ [1, 3], [2, 3], [3, 3],
+ [1, 2], [2, 2], [3, 2],
+ [1, 1], [2, 1], [3, 1]
+ ])
+ output_fraction = np.array([
+ 0, 0, 0.33,
+ 0, 0.22, 0,
+ 0.11, 0, 0.31
+ ])
+
+ output_intensity = np.array([
+ 0, 0, -0.33,
+ 0, -0.22, 0,
+ -0.11, 0, -0.31
+ ])
+
+ np.testing.assert_array_equal(
+ haz_read.centroids.lon,
+ output_raster[:, 0]
+ )
+ np.testing.assert_array_equal(
+ haz_read.centroids.lat,
+ output_raster[:, 1]
+ )
+ np.testing.assert_array_almost_equal(
+ haz_read.fraction.toarray().flatten(),
+ output_fraction
+ )
+ np.testing.assert_array_almost_equal(
+ haz_read.intensity.toarray().flatten(),
+ output_intensity
+ )
+
+ DATA_DIR.joinpath(intensity_file).unlink()
+ DATA_DIR.joinpath(fraction_file).unlink()
+
class TestStormEurope(unittest.TestCase):
@@ -270,7 +340,7 @@ def test_write_read_pass(self):
file_name = str(DATA_DIR.joinpath("test_haz.h5"))
# Read demo matlab file
- hazard = Hazard.from_mat(HAZ_TEST_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
hazard.event_name = list(map(str, hazard.event_name))
for todense_flag in [False, True]:
if todense_flag:
@@ -308,62 +378,16 @@ def test_write_read_pass(self):
)
self.assertIsInstance(haz_read.fraction, sparse.csr_matrix)
- def test_raster_to_vector_pass(self):
- """Test raster_to_vector method"""
+ def test_read_raster_pass(self):
+ """Test from_raster method"""
haz_fl = Hazard.from_raster([HAZ_DEMO_FL], haz_type="FL")
- haz_fl.check()
- meta_orig = haz_fl.centroids.meta
inten_orig = haz_fl.intensity
fract_orig = haz_fl.fraction
- haz_fl.raster_to_vector()
-
- self.assertEqual(haz_fl.centroids.meta, dict())
- self.assertAlmostEqual(
- haz_fl.centroids.lat.min(),
- meta_orig["transform"][5]
- + meta_orig["height"] * meta_orig["transform"][4]
- - meta_orig["transform"][4] / 2,
- )
- self.assertAlmostEqual(
- haz_fl.centroids.lat.max(),
- meta_orig["transform"][5] + meta_orig["transform"][4] / 2,
- )
- self.assertAlmostEqual(
- haz_fl.centroids.lon.max(),
- meta_orig["transform"][2]
- + meta_orig["width"] * meta_orig["transform"][0]
- - meta_orig["transform"][0] / 2,
- )
- self.assertAlmostEqual(
- haz_fl.centroids.lon.min(),
- meta_orig["transform"][2] + meta_orig["transform"][0] / 2,
- )
- self.assertTrue(u_coord.equal_crs(haz_fl.centroids.crs, meta_orig["crs"]))
self.assertTrue(np.allclose(haz_fl.intensity.data, inten_orig.data))
self.assertTrue(np.allclose(haz_fl.fraction.data, fract_orig.data))
- def test_reproject_raster_pass(self):
- """Test reproject_raster reference."""
-
- haz_fl = Hazard.from_raster([HAZ_DEMO_FL])
- haz_fl.check()
-
- haz_fl.reproject_raster(dst_crs="epsg:2202")
-
- self.assertEqual(haz_fl.intensity.shape, (1, 1046408))
- self.assertIsInstance(haz_fl.intensity, sparse.csr_matrix)
- self.assertIsInstance(haz_fl.fraction, sparse.csr_matrix)
- self.assertEqual(haz_fl.fraction.shape, (1, 1046408))
- self.assertTrue(u_coord.equal_crs(haz_fl.centroids.meta["crs"], "epsg:2202"))
- self.assertEqual(haz_fl.centroids.meta["width"], 968)
- self.assertEqual(haz_fl.centroids.meta["height"], 1081)
- self.assertEqual(haz_fl.fraction.min(), 0)
- self.assertEqual(haz_fl.fraction.max(), 1)
- self.assertEqual(haz_fl.intensity.min(), -9999)
- self.assertTrue(haz_fl.intensity.max() < 4.7)
-
# Execute Tests
if __name__ == "__main__":
diff --git a/climada/test/test_plot.py b/climada/test/test_plot.py
index 607fefdeab..dcfb608f98 100644
--- a/climada/test/test_plot.py
+++ b/climada/test/test_plot.py
@@ -26,19 +26,28 @@
import matplotlib.pyplot as plt
import pandas as pd
import contextily as ctx
+from pathlib import Path
from climada.engine.unsequa import UncOutput
from climada.engine import ImpactCalc, ImpactFreqCurve, CostBenefit
from climada.entity import (Entity, ImpactFuncSet, Exposures, DiscRates, ImpfTropCyclone, Measure,
MeasureSet)
from climada.hazard import Hazard, Centroids
-from climada.util.constants import HAZ_DEMO_MAT, ENT_DEMO_TODAY, TEST_UNC_OUTPUT_COSTBEN
+from climada.util.constants import ENT_DEMO_TODAY, TEST_UNC_OUTPUT_COSTBEN, HAZ_DEMO_FL
from climada.util.api_client import Client
+from climada.test import get_test_file
test_unc_output_costben = Client().get_dataset_file(name=TEST_UNC_OUTPUT_COSTBEN, status='test_dataset')
+
+HAZ_TEST_TC :Path = get_test_file('test_tc_florida')
+"""
+Hazard test file from Data API: Hurricanes from 1851 to 2011 over Florida with 100 centroids.
+Fraction is empty. Format: HDF5.
+"""
+
class TestPlotter(unittest.TestCase):
"""Test plot functions."""
@@ -47,7 +56,7 @@ def setUp(self):
def test_hazard_intensity_pass(self):
"""Generate all possible plots of the hazard intensity."""
- hazard = Hazard.from_mat(HAZ_DEMO_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
hazard.event_name = [""] * hazard.event_id.size
hazard.event_name[35] = "NNN_1185106_gen5"
hazard.event_name[3898] = "NNN_1190604_gen8"
@@ -81,34 +90,22 @@ def test_hazard_intensity_pass(self):
def test_hazard_fraction_pass(self):
"""Generate all possible plots of the hazard fraction."""
- hazard = Hazard.from_mat(HAZ_DEMO_MAT)
+ hazard = Hazard.from_raster([HAZ_DEMO_FL, HAZ_DEMO_FL])
hazard.event_name = [""] * hazard.event_id.size
- hazard.event_name[35] = "NNN_1185106_gen5"
- hazard.event_name[11897] = "GORDON_gen7"
- myax = hazard.plot_fraction(event=36)
- self.assertIn('Event ID 36: NNN_1185106_gen5', myax.get_title())
-
- myax = hazard.plot_fraction(event=-1)
- self.assertIn('1-largest Event. ID 11898: GORDON_gen7', myax.get_title())
+ hazard.event_name[0] = "NNN_1185106_gen5"
+ myax = hazard.plot_fraction(event=1)
+ self.assertIn('Event ID 1: NNN_1185106_gen5', myax.get_title())
- myax = hazard.plot_fraction(centr=59)
- self.assertIn('Centroid 59: (30.0, -79.0)', myax.get_title())
-
- myax = hazard.plot_fraction(centr=-1)
- self.assertIn('1-largest Centroid. 79: (30.0, -77.0)', myax.get_title())
+ myax = hazard.plot_fraction(centr=1)
+ self.assertIn('Centroid 1: (10.424, -69.324)', myax.get_title())
def test_hazard_rp_intensity(self):
""""Plot exceedance intensity maps for different return periods"""
- hazard = Hazard.from_mat(HAZ_DEMO_MAT)
+ hazard = Hazard.from_hdf5(HAZ_TEST_TC)
(axis1, axis2), _ = hazard.plot_rp_intensity([25, 50])
self.assertEqual('Return period: 25 years', axis1.get_title())
self.assertEqual('Return period: 50 years', axis2.get_title())
- def test_hazard_centroids(self):
- """Plot centroids scatter points over earth."""
- centroids = Centroids.from_base_grid()
- centroids.plot()
-
def test_exposures_value_pass(self):
"""Plot exposures values."""
myexp = pd.read_excel(ENT_DEMO_TODAY)
@@ -140,7 +137,7 @@ def test_impact_pass(self):
"""Plot impact exceedence frequency curves."""
myent = Entity.from_excel(ENT_DEMO_TODAY)
myent.exposures.check()
- myhaz = Hazard.from_mat(HAZ_DEMO_MAT)
+ myhaz = Hazard.from_hdf5(HAZ_TEST_TC)
myhaz.event_name = [""] * myhaz.event_id.size
myimp = ImpactCalc(myent.exposures, myent.impact_funcs, myhaz).impact()
ifc = myimp.calc_freq_curve()
diff --git a/climada/util/coordinates.py b/climada/util/coordinates.py
index 741b405864..8b9113c5bf 100644
--- a/climada/util/coordinates.py
+++ b/climada/util/coordinates.py
@@ -68,9 +68,6 @@
NE_CRS = f"epsg:{NE_EPSG}"
"""Natural Earth CRS"""
-TMP_ELEVATION_FILE = SYSTEM_DIR.joinpath('tmp_elevation.tif')
-"""Path of elevation file written in set_elevation"""
-
DEM_NODATA = -9999
"""Value to use for no data values in DEM, i.e see points"""
@@ -756,7 +753,7 @@ def nat_earth_resolution(resolution):
raise ValueError('Natural Earth does not accept resolution %s m.' % resolution)
return str(resolution) + 'm'
-def get_country_geometries(country_names=None, extent=None, resolution=10):
+def get_country_geometries(country_names=None, extent=None, resolution=10, center_crs=True):
"""Natural Earth country boundaries within given extent
If no arguments are given, simply returns the whole natural earth dataset.
@@ -765,8 +762,9 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
starts including the projection information. (They are saving a whopping 147 bytes by omitting
it.) Same goes for UTF.
- If extent is provided, longitude values in 'geom' will all lie within 'extent' longitude
- range. Therefore setting extent to e.g. [160, 200, -20, 20] will provide longitude values
+ If extent is provided and center_crs is True, longitude values in 'geom' will all lie
+ within 'extent' longitude range.
+ Therefore setting extent to e.g. [160, 200, -20, 20] will provide longitude values
between 160 and 200 degrees.
Parameters
@@ -778,6 +776,10 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
Extent, assumed to be in the same CRS as the natural earth data.
resolution : float, optional
10, 50 or 110. Resolution in m. Default: 10m
+ center_crs : bool
+ if True, the crs of the countries is centered such that
+ longitude values in 'geom' will all lie within 'extent' longitude range.
+ Default is True.
Returns
-------
@@ -839,7 +841,7 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
bbox = gpd.GeoSeries(bbox, crs=DEF_CRS)
bbox = gpd.GeoDataFrame({'geometry': bbox}, crs=DEF_CRS)
out = gpd.overlay(out, bbox, how="intersection")
- if ~lon_normalized:
+ if ~lon_normalized and center_crs:
lon_mid = 0.5 * (extent[0] + extent[1])
# reset the CRS attribute after rewrapping (we don't really change the CRS)
out = (
@@ -984,7 +986,7 @@ def assign_coordinates(*args, **kwargs):
return match_coordinates(*args, **kwargs)
def match_coordinates(coords, coords_to_assign, distance="euclidean",
- threshold=NEAREST_NEIGHBOR_THRESHOLD, **kwargs):
+ threshold=NEAREST_NEIGHBOR_THRESHOLD, **kwargs):
"""To each coordinate in `coords`, assign a matching coordinate in `coords_to_assign`
If there is no exact match for some entry, an attempt is made to assign the geographically
@@ -1082,7 +1084,7 @@ def match_coordinates(coords, coords_to_assign, distance="euclidean",
def match_centroids(coord_gdf, centroids, distance='euclidean',
threshold=NEAREST_NEIGHBOR_THRESHOLD):
"""Assign to each gdf coordinate point its closest centroids's coordinate.
- If disatances > threshold in points' distances, -1 is returned.
+ If distances > threshold in points' distances, -1 is returned.
If centroids are in a raster and coordinate point is outside of it ``-1`` is assigned
Parameters
@@ -1133,15 +1135,10 @@ def match_centroids(coord_gdf, centroids, distance='euclidean',
# no error is raised and it is assumed that the user set the crs correctly
pass
- if centroids.meta:
- assigned = match_grid_points(
- coord_gdf.longitude.values, coord_gdf.latitude.values,
- centroids.meta['width'], centroids.meta['height'],
- centroids.meta['transform'])
- else:
- assigned = match_coordinates(
- np.stack([coord_gdf.latitude.values, coord_gdf.longitude.values], axis=1),
- centroids.coord, distance=distance, threshold=threshold)
+ assigned = match_coordinates(
+ np.stack([coord_gdf['latitude'].values, coord_gdf['longitude'].values], axis=1),
+ centroids.coord, distance=distance, threshold=threshold,
+ )
return assigned
@numba.njit
@@ -1549,7 +1546,8 @@ def get_country_code(lat, lon, gridded=False):
region_id = region_id.astype(int)
else:
(lon_min, lat_min, lon_max, lat_max) = latlon_bounds(lat, lon, 0.001)
- countries = get_country_geometries(extent=(lon_min, lon_max, lat_min, lat_max))
+ countries = get_country_geometries(
+ extent=(lon_min, lon_max, lat_min, lat_max), center_crs=False)
with warnings.catch_warnings():
# in order to suppress the following
# UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely
diff --git a/climada/util/test/test_coordinates.py b/climada/util/test/test_coordinates.py
index 15b8b9119a..7026e629d7 100644
--- a/climada/util/test/test_coordinates.py
+++ b/climada/util/test/test_coordinates.py
@@ -550,30 +550,33 @@ def test_match_centroids(self):
'width': 20, 'height': 10,
'transform': rasterio.Affine(1.5, 0.0, -20, 0.0, -1.4, 8)
}
- centroids = Centroids(meta=meta)
+ centroids = Centroids.from_meta(meta=meta)
df = pd.DataFrame({
'longitude': np.array([
-20.1, -20.0, -19.8, -19.0, -18.6, -18.4,
-19.0, -19.0, -19.0, -19.0,
-20.1, 0.0, 10.1, 10.1, 10.1, 0.0, -20.2, -20.3,
-6.4, 9.8, 0.0,
- ]),
+ ]),
'latitude': np.array([
7.3, 7.3, 7.3, 7.3, 7.3, 7.3,
8.1, 7.9, 6.7, 6.5,
8.1, 8.2, 8.3, 0.0, -6.1, -6.2, -6.3, 0.0,
-1.9, -1.7, 0.0,
- ])
+ ]),
})
- gdf = gpd.GeoDataFrame(df,geometry=gpd.points_from_xy(df.longitude, df.latitude),
- crs=DEF_CRS)
- assigned = u_coord.match_centroids(gdf,centroids)
+ gdf = gpd.GeoDataFrame(
+ df,
+ geometry=gpd.points_from_xy(df.longitude, df.latitude),
+ crs=DEF_CRS,
+ )
+ assigned = u_coord.match_centroids(gdf, centroids)
expected_result = [
# constant y-value, varying x-value
- -1, 0, 0, 0, 0, 1,
+ 0, 0, 0, 0, 0, 1,
# constant x-value, varying y-value
- -1, 0, 0, 20,
+ 0, 0, 0, 20,
# out of bounds: topleft, top, topright, right, bottomright, bottom, bottomleft, left
-1, -1, -1, -1, -1, -1, -1, -1,
# some explicit points within the raster
@@ -595,9 +598,9 @@ def test_match_centroids(self):
centroids = Centroids(
lat=coords_to_assign[:, 0],
lon=coords_to_assign[:, 1],
- geometry = gpd.GeoSeries(crs=DEF_CRS)
+ crs=DEF_CRS
)
- centroids_empty = Centroids()
+ centroids_empty = Centroids(lat=np.array([]), lon=np.array([]))
expected_results = [
# test with different thresholds (in km)
@@ -630,7 +633,7 @@ def test_match_centroids(self):
centroids = Centroids(
lat=[1100000,1200000],
lon=[2500000,2600000],
- geometry = gpd.GeoSeries(crs='EPSG:2056')
+ crs='EPSG:2056'
)
with self.assertRaises(ValueError) as cm:
diff --git a/doc/tutorial/climada_hazard_Hazard.ipynb b/doc/tutorial/climada_hazard_Hazard.ipynb
index 2375b5d97c..6785a80b7e 100644
--- a/doc/tutorial/climada_hazard_Hazard.ipynb
+++ b/doc/tutorial/climada_hazard_Hazard.ipynb
@@ -65,8 +65,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- " \n",
- "## Part 1: Read hazards from raster data\n",
+ "\n",
+ "## Part 1: Read hazards from raster data\n",
"\n",
"Raster data can be read in any format accepted by [rasterio](https://rasterio.readthedocs.io/en/stable/) using `Hazard`'s `from_raster()` method. The raster information might refer to the `intensity` or `fraction`of the hazard. Different configuration options such as transforming the coordinates, changing the CRS and reading only a selected area or band are available through the `from_raster()` arguments as follows:"
]
@@ -195,7 +195,7 @@
"haz.check()\n",
"print('\\n Solution 1:')\n",
"print('centroids CRS:', haz.centroids.crs)\n",
- "print('raster info:', haz.centroids.meta)\n",
+ "print('raster info:', haz.centroids.get_meta())\n",
"\n",
"# 2. Transformations of the coordinates can be set using the transform option and Affine\n",
"from rasterio import Affine\n",
@@ -205,7 +205,7 @@
" height=500, width=501)\n",
"haz.check()\n",
"print('\\n Solution 2:')\n",
- "print('raster info:', haz.centroids.meta)\n",
+ "print('raster info:', haz.centroids.get_meta())\n",
"print('intensity size:', haz.intensity.shape)\n",
"\n",
"# 3. A partial part of the raster can be loaded using the window or geometry\n",
@@ -213,7 +213,7 @@
"haz = Hazard.from_raster([HAZ_DEMO_FL], haz_type='FL', window=Window(10, 10, 20, 30))\n",
"haz.check()\n",
"print('\\n Solution 3:')\n",
- "print('raster info:', haz.centroids.meta)\n",
+ "print('raster info:', haz.centroids.get_meta())\n",
"print('intensity size:', haz.intensity.shape)"
]
},
@@ -221,13 +221,13 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- " \n",
- "## Part 2: Read hazards from other data\n",
+ "\n",
+ "## Part 2: Read hazards from other data\n",
"\n",
"- excel: Hazards can be read from Excel files following the template in `climada_python/climada/data/system/hazard_template.xlsx` using the `from_excel()` method. \n",
"- MATLAB: Hazards generated with CLIMADA's MATLAB version (.mat format) can be read using `from_mat()`.\n",
"- vector data: Use `Hazard`'s `from_vector`-constructor to read shape data (all formats supported by [fiona](https://fiona.readthedocs.io/en/latest/manual.html)).\n",
- "- hdf5: Hazards generated with the CLIMADA in Python (.h5 format) can be read using `from_hdf5()`."
+ "- hdf5: Hazards generated with the CLIMADA in Python (.h5 format) can be read using `from_hdf5()`.\n"
]
},
{
@@ -255,8 +255,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- " \n",
- "## Part 3: Define hazards manually\n",
+ "\n",
+ "## Part 3: Define hazards manually\n",
+ "\n",
"A `Hazard` can be defined by filling its values one by one, as follows:"
]
},
@@ -314,7 +315,7 @@
"haz = Hazard(haz_type='TC',\n",
" intensity=intensity,\n",
" fraction=fraction,\n",
- " centroids=Centroids.from_lat_lon(lat, lon), # default crs used\n",
+ " centroids=Centroids(lat=lat, lon=lon), # default crs used\n",
" units='m',\n",
" event_id=np.arange(n_ev, dtype=int),\n",
" event_name=['ev_12', 'ev_21', 'Maria', 'ev_35',\n",
@@ -328,17 +329,86 @@
"haz.centroids.plot();"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Or the `Hazard` can be defined with a grid:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# using from_pnt_bounds\n",
+ "\n",
+ "# bounds\n",
+ "left, bottom, right, top = -72, -3.0, -52.0, 22 # the bounds refer to the bounds of the center of the pixel\n",
+ "# resolution\n",
+ "res = 0.5\n",
+ "centroids = Centroids.from_pnt_bounds((left, bottom, right, top), res) # default crs used"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "True"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# the same can be done with the method `from_meta`, by definition of a raster meta object\n",
+ "\n",
+ "import rasterio\n",
+ "from climada.util.constants import DEF_CRS\n",
+ "\n",
+ "# raster info:\n",
+ "# border upper left corner (of the pixel, not of the center of the pixel)\n",
+ "max_lat = top + res/2\n",
+ "min_lon = left - res/2\n",
+ "# resolution in lat and lon\n",
+ "d_lat = -res # negative because starting in upper corner\n",
+ "d_lon = res # same step as d_lat\n",
+ "# number of points\n",
+ "n_lat, n_lon = centroids.shape\n",
+ "\n",
+ "# meta: raster specification\n",
+ "meta = {\n",
+ " 'dtype': 'float32',\n",
+ " 'width': n_lon,\n",
+ " 'height': n_lat,\n",
+ " 'crs': DEF_CRS,\n",
+ " 'transform': rasterio.Affine(\n",
+ " a=d_lon, b=0.0, c=min_lon,\n",
+ " d=0.0, e=d_lat, f=max_lat),\n",
+ "}\n",
+ "\n",
+ "centroids_from_meta = Centroids.from_meta(meta) # default crs used\n",
+ "\n",
+ "centroids_from_meta == centroids"
+ ]
+ },
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
- "Check centroids borders: (-72.0, -3.0, -52.0, 22.0)\n",
- "Check centroids borders: (-72.25, -3.25, -51.75, 22.25)\n"
+ "Check centroids borders: [-72. -3. -52. 22.]\n"
]
},
{
@@ -355,23 +425,13 @@
}
],
"source": [
- "# setting raster\n",
+ "# create a Hazard object with random events\n",
+ "\n",
"import numpy as np\n",
"from scipy import sparse\n",
"\n",
- "# raster info:\n",
- "# border upper left corner (of the pixel, not of the center of the pixel)\n",
- "xf_lat = 22\n",
- "xo_lon = -72\n",
- "# resolution in lat and lon\n",
- "d_lat = -0.5 # negative because starting in upper corner\n",
- "d_lon = 0.5 # same step as d_lat\n",
- "# number of points\n",
- "n_lat = 50\n",
- "n_lon = 40\n",
- "\n",
"n_ev = 10 # number of events\n",
- "centroids = Centroids.from_pix_bounds(xf_lat, xo_lon, d_lat, d_lon, n_lat, n_lon) # default crs used\n",
+ "\n",
"intensity = sparse.csr_matrix(np.random.random((n_ev, centroids.size)))\n",
"fraction = intensity.copy()\n",
"fraction.data.fill(1)\n",
@@ -391,12 +451,7 @@
"\n",
"haz.check()\n",
"print('Check centroids borders:', haz.centroids.total_bounds)\n",
- "haz.centroids.plot()\n",
- "\n",
- "# using from_pnt_bounds, the bounds refer to the bounds of the center of the pixel\n",
- "left, bottom, right, top = xo_lon, -3.0, -52.0, xf_lat\n",
- "haz.centroids = Centroids.from_pnt_bounds((left, bottom, right, top), 0.5) # default crs used\n",
- "print('Check centroids borders:', haz.centroids.total_bounds)"
+ "haz.centroids.plot();"
]
},
{
@@ -418,7 +473,7 @@
"- `select()` returns a new hazard with the selected region, date and/or synthetic or historical filter.\n",
"- `remove_duplicates()` removes events with same name and date.\n",
"- `local_exceedance_inten()` returns a matrix with the exceedence frequency at every frequency and provided return periods. This is the one used in `plot_rp_intensity()`.\n",
- "- `reproject_raster()`, `reproject_vector()`, `raster_to_vector()`, `vector_to_raster()` are methods to change centroids' CRS and between raster and vector data. \n",
+ "- `reproject_vector()` is a method to change the centroids' CRS. \n",
"\n",
"`Centroids` methods:\n",
"- centroids properties such as area per pixel, distance to coast, country ISO code, on land mask or elevation are available through different `set_XX()`methods.\n",
@@ -496,8 +551,7 @@
"print('Year with most hurricanes between 1995 and 2001:', max_year)\n",
"\n",
"# 4. What is the number of centroids with distance to coast smaller than 1km?\n",
- "hist_tc.centroids.set_dist_coast()\n",
- "num_cen_coast = np.argwhere(hist_tc.centroids.dist_coast < 1000).size\n",
+ "num_cen_coast = np.argwhere(hist_tc.centroids.get_dist_coast() < 1000).size\n",
"print('Number of centroids close to coast: ', num_cen_coast)"
]
},
@@ -719,7 +773,7 @@
"# 7. one figure with two plots: maximum intensities and selected centroid with all intensities:\n",
"from climada.util.plot import make_map\n",
"import matplotlib.pyplot as plt\n",
- "plt.ioff()\n",
+ "\n",
"fig, ax1, fontsize = make_map(1) # map\n",
"ax2 = fig.add_subplot(2, 1, 2) # add regular axes\n",
"haz_tc_fl.plot_intensity(axis=ax1, event=0) # plot original resolution\n",
diff --git a/script/jenkins/petals_regression_test/Jenkinsfile b/script/jenkins/petals_regression_test/Jenkinsfile
index 56a8b0241a..433771e3ab 100644
--- a/script/jenkins/petals_regression_test/Jenkinsfile
+++ b/script/jenkins/petals_regression_test/Jenkinsfile
@@ -4,20 +4,7 @@ pipeline {
stages {
stage('integ_test') {
steps {
- sh '''#!/bin/bash
- export PATH=$PATH:$CONDAPATH
- mamba env update -n climada_env -f ~/jobs/petals_install_env/workspace/requirements/env_climada.yml
-
- source activate climada_env
- pip install -e ~/jobs/petals_install_env/workspace/
-
- cp ~/jobs/petals_install_env/workspace/climada.conf climada.conf
- python script/jenkins/set_config.py test_directory ~/jobs/petals_install_env/workspace/climada_petals
-
- PYTHONPATH=.:$PYTHONPATH pytest --junitxml=tests_xml/tests.xml ~/jobs/petals_install_env/workspace/climada_petals
-
- git checkout climada.conf
- '''
+ sh 'bash script/jenkins/petals_regression_test/run_integ_test.sh'
}
}
}
diff --git a/script/jenkins/petals_regression_test/run_integ_test.sh b/script/jenkins/petals_regression_test/run_integ_test.sh
new file mode 100644
index 0000000000..9f9844bf60
--- /dev/null
+++ b/script/jenkins/petals_regression_test/run_integ_test.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+export PATH=$PATH:$CONDAPATH
+mamba env update -n climada_env -f ~/jobs/petals_install_env/workspace/requirements/env_climada.yml
+
+source activate climada_env
+
+REGTESTENV=~/jobs/petals_compatibility/petals_env
+BRANCH=`git name-rev --name-only HEAD | cut -f 3- -d /`
+PETALS_DIR=`test -e $REGTESTENV/$BRANCH && cat $REGTESTENV/$BRANCH || echo ~/jobs/petals_branches/branches/develop/workspace`
+
+python -m venv --system-site-packages tvenv
+source tvenv/bin/activate
+
+pip install -e $PETALS_DIR
+
+cp $PETALS_DIR/climada.conf climada.conf
+python script/jenkins/set_config.py test_directory $PETALS_DIR/climada_petals
+
+PYTHONPATH=.:$PYTHONPATH pytest --junitxml=tests_xml/tests.xml $PETALS_DIR/climada_petals
+
+git checkout climada.conf
+
+deactivate
+#rm -r tvenv