diff --git a/climada/hazard/tc_tracks.py b/climada/hazard/tc_tracks.py index f6d970cfc..71b9c6955 100644 --- a/climada/hazard/tc_tracks.py +++ b/climada/hazard/tc_tracks.py @@ -46,14 +46,13 @@ import pandas as pd import pathos import scipy.io.matlab as matlab -import shapely.ops +import shapely import statsmodels.api as sm import xarray as xr from matplotlib.collections import LineCollection from matplotlib.colors import BoundaryNorm, ListedColormap from matplotlib.lines import Line2D from shapely.geometry import LineString, MultiLineString, Point, Polygon -from shapely.ops import unary_union from sklearn.metrics import DistanceMetric from tqdm import tqdm @@ -264,7 +263,7 @@ class BasinBoundsStorm(Enum): [(10.0, -60.0), (135.0, -60.0), (135.0, -5.0), (10.0, -5.0), (10.0, -60.0)] ) - SP = unary_union( + SP = shapely.union_all( [ Polygon( # west side of antimeridian [ @@ -642,7 +641,7 @@ def tracks_in_exp(self, exposure, buffer=1.0): raise Exception("this is not an Exposures object") exp_buffer = exposure.gdf.buffer(distance=buffer, resolution=0) - exp_buffer = exp_buffer.unary_union + exp_buffer = exp_buffer.union_all() tc_tracks_lines = self.to_geodataframe().buffer(distance=buffer) select_tracks = tc_tracks_lines.intersects(exp_buffer) diff --git a/climada/hazard/test/test_tc_tracks.py b/climada/hazard/test/test_tc_tracks.py index 52535c34c..7ebf74a26 100644 --- a/climada/hazard/test/test_tc_tracks.py +++ b/climada/hazard/test/test_tc_tracks.py @@ -714,14 +714,14 @@ def test_to_geodataframe_points(self): tc_track = tc.TCTracks.from_processed_ibtracs_csv(TEST_TRACK) gdf_points = tc_track.to_geodataframe(as_points=True) - self.assertIsInstance(gdf_points.unary_union.bounds, tuple) + self.assertIsInstance(gdf_points.union_all().bounds, tuple) self.assertEqual(gdf_points.shape[0], len(tc_track.data[0]["time"])) self.assertEqual( gdf_points.shape[1], len(tc_track.data[0].variables) + len(tc_track.data[0].attrs) - 1, ) self.assertAlmostEqual( - gdf_points.buffer(3).unary_union.area, 348.79972062947854 + gdf_points.buffer(3).union_all().area, 348.79972062947854 ) self.assertIsInstance( gdf_points.iloc[0].time, pd._libs.tslibs.timestamps.Timestamp diff --git a/climada/util/coordinates.py b/climada/util/coordinates.py index f76fc11b7..2c6320383 100644 --- a/climada/util/coordinates.py +++ b/climada/util/coordinates.py @@ -42,8 +42,7 @@ import rasterio.warp import scipy.interpolate import scipy.spatial -import shapely.ops -import shapely.vectorized +import shapely import shapely.wkt from cartopy.io import shapereader from pyproj.crs import CRS as PCRS @@ -753,7 +752,7 @@ def get_land_geometry(country_names=None, extent=None, resolution=10): """ geom = get_country_geometries(country_names, extent, resolution) # combine all into a single multipolygon - geom = geom.geometry.unary_union + geom = geom.geometry.union_all() if not isinstance(geom, MultiPolygon): geom = MultiPolygon([geom]) return geom @@ -804,7 +803,7 @@ def coord_on_land(lat, lon, land_geom=None): lon_mid = 0.5 * (land_bounds[0] + land_bounds[2]) lon_normalize(lons, center=lon_mid) - return shapely.vectorized.contains(land_geom, lons, lat) + return shapely.contains_xy(land_geom, lons, lat) def nat_earth_resolution(resolution): @@ -917,7 +916,7 @@ def get_country_geometries( lon_left, lon_right = lon_normalize(np.array(extent[:2])) extent_left = (lon_left, 180, extent[2], extent[3]) extent_right = (-180, lon_right, extent[2], extent[3]) - bbox = shapely.ops.unary_union( + bbox = shapely.union_all( [box(*toggle_extent_bounds(e)) for e in [extent_left, extent_right]] ) bbox = gpd.GeoSeries(bbox, crs=DEF_CRS) @@ -1529,7 +1528,7 @@ def _nearest_neighbor_approx( if num_warn: LOGGER.warning( - "Distance to closest centroid is greater than %s" "km for %s coordinates.", + "Distance to closest centroid is greater than %skm for %s coordinates.", threshold, num_warn, ) @@ -1723,8 +1722,7 @@ def _nearest_neighbor_antimeridian(centroids, coordinates, threshold, assigned, lon_max = max(centroids[:, 1].max(), coordinates[:, 1].max()) if lon_max - lon_min > 360: raise ValueError( - "Longitudinal coordinates need to be normalized" - "to a common 360 degree range" + "Longitudinal coordinates need to be normalizedto a common 360 degree range" ) mid_lon = 0.5 * (lon_max + lon_min) antimeridian = mid_lon + 180 @@ -1970,18 +1968,16 @@ def get_country_code(lat, lon, gridded=False): countries["area"] = countries.geometry.area countries = countries.sort_values(by=["area"], ascending=False) region_id = np.full((lon.size,), -1, dtype=int) - total_land = countries.geometry.unary_union + total_land = countries.geometry.union_all() ocean_mask = ( region_id.all() if total_land is None - else ~shapely.vectorized.contains(total_land, lon, lat) + else ~shapely.contains_xy(total_land, lon, lat) ) region_id[ocean_mask] = 0 for country in countries.itertuples(): unset = (region_id == -1).nonzero()[0] - select = shapely.vectorized.contains( - country.geometry, lon[unset], lat[unset] - ) + select = shapely.contains_xy(country.geometry, lon[unset], lat[unset]) region_id[unset[select]] = natearth_country_to_int(country) region_id[region_id == -1] = 0 return region_id @@ -3209,8 +3205,9 @@ def subraster_from_bounds(transform, bounds): # align the window bounds to the raster by rounding col_min, col_max = np.round(window.col_off), np.round(window.col_off + window.width) - row_min, row_max = np.round(window.row_off), np.round( - window.row_off + window.height + row_min, row_max = ( + np.round(window.row_off), + np.round(window.row_off + window.height), ) window = rasterio.windows.Window( col_min, row_min, col_max - col_min, row_max - row_min diff --git a/climada/util/lines_polys_handler.py b/climada/util/lines_polys_handler.py index ee2058a68..8bcfd7899 100755 --- a/climada/util/lines_polys_handler.py +++ b/climada/util/lines_polys_handler.py @@ -764,7 +764,7 @@ def _interp_one_poly_grid(poly, x_grid, y_grid): if poly.is_empty: return shgeom.MultiPoint([]) - in_geom = sh.vectorized.contains(poly, x_grid, y_grid) + in_geom = sh.contains_xy(poly, x_grid, y_grid) if sum(in_geom.flatten()) > 1: return shgeom.MultiPoint(list(zip(x_grid[in_geom], y_grid[in_geom]))) @@ -796,7 +796,7 @@ def _interp_one_poly(poly, res): height, width, trafo = u_coord.pts_to_raster_meta(poly.bounds, (res, res)) x_grid, y_grid = u_coord.raster_to_meshgrid(trafo, width, height) - in_geom = sh.vectorized.contains(poly, x_grid, y_grid) + in_geom = sh.contains_xy(poly, x_grid, y_grid) if sum(in_geom.flatten()) > 1: return shgeom.MultiPoint(list(zip(x_grid[in_geom], y_grid[in_geom]))) @@ -835,7 +835,7 @@ def _interp_one_poly_m(poly, res, orig_crs): height, width, trafo = u_coord.pts_to_raster_meta(poly_m.bounds, (res, res)) x_grid, y_grid = u_coord.raster_to_meshgrid(trafo, width, height) - in_geom = sh.vectorized.contains(poly_m, x_grid, y_grid) + in_geom = sh.contains_xy(poly_m, x_grid, y_grid) if sum(in_geom.flatten()) > 1: x_poly, y_poly = reproject_grid( x_grid[in_geom], y_grid[in_geom], m_crs, orig_crs