From 6a6965823c017ee634a3d98d09efa31eb94e49a1 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Thu, 22 Sep 2022 23:08:28 +0200 Subject: [PATCH] Make GeoPandas work with Shapely 2.0 (#2275) --- ci/envs/310-dev.yaml | 3 +- ci/envs/39-latest-conda-forge.yaml | 6 +- doc/source/getting_started/install.rst | 26 +++- geopandas/_compat.py | 22 ++- geopandas/_vectorized.py | 206 +++++++++++++++++++------ geopandas/array.py | 16 +- geopandas/geoseries.py | 16 +- geopandas/io/sql.py | 15 +- geopandas/io/tests/test_pickle.py | 7 +- geopandas/sindex.py | 58 ++++--- geopandas/tests/test_array.py | 10 +- geopandas/tests/test_geodataframe.py | 4 +- geopandas/tests/test_geom_methods.py | 2 +- geopandas/tests/test_sindex.py | 30 ++-- geopandas/tools/sjoin.py | 6 +- geopandas/tools/tests/test_sjoin.py | 4 +- 16 files changed, 313 insertions(+), 118 deletions(-) diff --git a/ci/envs/310-dev.yaml b/ci/envs/310-dev.yaml index 0c250d6..8efed12 100644 --- a/ci/envs/310-dev.yaml +++ b/ci/envs/310-dev.yaml @@ -5,7 +5,6 @@ dependencies: - python=3.10 - cython # required - - shapely - pyproj - geos - packaging @@ -30,7 +29,7 @@ dependencies: - fiona - git+https://github.com/pandas-dev/pandas.git@main - git+https://github.com/matplotlib/matplotlib.git@main - # - git+https://github.com/Toblerity/Shapely.git@main + - git+https://github.com/shapely/shapely.git@main - git+https://github.com/pygeos/pygeos.git@master - git+https://github.com/python-visualization/folium.git@main - git+https://github.com/geopandas/xyzservices.git@main diff --git a/ci/envs/39-latest-conda-forge.yaml b/ci/envs/39-latest-conda-forge.yaml index 49c7f64..6a61ddb 100644 --- a/ci/envs/39-latest-conda-forge.yaml +++ b/ci/envs/39-latest-conda-forge.yaml @@ -1,14 +1,16 @@ name: test channels: - conda-forge + - conda-forge/label/shapely_dev dependencies: - python=3.9 # required - pandas=1.3 - - shapely + - shapely=2 - fiona - pyproj - - pygeos + # use this build to have one with only shapely 2.0 installed + # - pygeos - packaging # testing - pytest diff --git a/doc/source/getting_started/install.rst b/doc/source/getting_started/install.rst index db6b125..07ff8d6 100644 --- a/doc/source/getting_started/install.rst +++ b/doc/source/getting_started/install.rst @@ -161,9 +161,19 @@ For plotting, these additional packages may be used: Using the optional PyGEOS dependency ------------------------------------ +.. attention:: + + The upcoming Shapely 2.0 release will absorb all improvements from PyGEOS. + If you are considering trying out those improvements, you can also test + the prerelease of Shapely instead. + See https://shapely.readthedocs.io/en/latest/release/2.x.html#version-2-0-0 + for the release notes of Shapely 2.0, and https://github.com/shapely/shapely/discussions/1464 + on how to install this and give feedback. + Work is ongoing to improve the performance of GeoPandas. Currently, the fast implementations of basic spatial operations live in the `PyGEOS`_ -package (but work is under way to contribute those improvements to Shapely). +package (but work is under way to contribute those improvements to Shapely, +coming to Shapely 2.0). Starting with GeoPandas 0.8, it is possible to optionally use those experimental speedups by installing PyGEOS. This can be done with conda (using the conda-forge channel) or pip:: @@ -182,11 +192,23 @@ More specifically, whether the speedups are used or not is determined by: - You can still toggle the use of PyGEOS when it is available, by: - Setting an environment variable (``USE_PYGEOS=0/1``). Note this variable - is only checked at first import of GeoPandas. + is only checked at first import of GeoPandas. You can set this environment + variable before starting the python process, or in your code right before + importing geopandas: + + .. code-block:: python + + import os + os.environ["USE_PYGEOS"] = "0" + import geopandas + - Setting an option: ``geopandas.options.use_pygeos = True/False``. Note, although this variable can be set during an interactive session, it will only work if the GeoDataFrames you use are created (e.g. reading a file with ``read_file``) after changing this value. + Attention: changing this option will no longer work in all cases when + having Shapely >=2.0 installed. In that case, use the environment variable + (see option above). .. warning:: diff --git a/geopandas/_compat.py b/geopandas/_compat.py index dd67212..f7673fa 100644 --- a/geopandas/_compat.py +++ b/geopandas/_compat.py @@ -29,13 +29,15 @@ PANDAS_GE_14 = Version(pd.__version__) >= Version("1.4.0rc0") SHAPELY_GE_18 = Version(shapely.__version__) >= Version("1.8") SHAPELY_GE_182 = Version(shapely.__version__) >= Version("1.8.2") -SHAPELY_GE_20 = Version(shapely.__version__) >= Version("2.0") +SHAPELY_GE_20 = Version(shapely.__version__) >= Version("2.0.0.dev0") +SHAPELY_G_20a1 = Version(shapely.__version__) > Version("2.0a1") GEOS_GE_390 = shapely.geos.geos_version >= (3, 9, 0) HAS_PYGEOS = None USE_PYGEOS = None +USE_SHAPELY_20 = None PYGEOS_SHAPELY_COMPAT = None PYGEOS_GE_09 = None @@ -74,6 +76,7 @@ def set_use_pygeos(val=None): Alternatively, pass a value here to force a True/False value. """ global USE_PYGEOS + global USE_SHAPELY_20 global PYGEOS_SHAPELY_COMPAT if val is not None: @@ -94,11 +97,18 @@ def set_use_pygeos(val=None): # validate the pygeos version if not Version(pygeos.__version__) >= Version("0.8"): - raise ImportError( - "PyGEOS >= 0.8 is required, version {0} is installed".format( - pygeos.__version__ + if SHAPELY_GE_20: + USE_PYGEOS = False + warnings.warn( + "The PyGEOS version is too old, and Shapely >= 2 is installed, " + "thus using Shapely by default and not PyGEOS." + ) + else: + raise ImportError( + "PyGEOS >= 0.8 is required, version {0} is installed".format( + pygeos.__version__ + ) ) - ) # Check whether Shapely and PyGEOS use the same GEOS version. # Based on PyGEOS from_shapely implementation. @@ -123,6 +133,8 @@ def set_use_pygeos(val=None): except ImportError: raise ImportError(INSTALL_PYGEOS_ERROR) + USE_SHAPELY_20 = (not USE_PYGEOS) and SHAPELY_GE_20 + set_use_pygeos() diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py index 159e578..90b6dd1 100644 --- a/geopandas/_vectorized.py +++ b/geopandas/_vectorized.py @@ -9,6 +9,7 @@ import warnings import numpy as np import pandas as pd +import shapely import shapely.geometry import shapely.geos import shapely.ops @@ -38,8 +39,11 @@ _names = { "GEOMETRYCOLLECTION": "GeometryCollection", } -if compat.USE_PYGEOS: - type_mapping = {p.value: _names[p.name] for p in pygeos.GeometryType} +if compat.USE_SHAPELY_20 or compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + type_mapping = {p.value: _names[p.name] for p in shapely.GeometryType} + else: + type_mapping = {p.value: _names[p.name] for p in pygeos.GeometryType} geometry_type_ids = list(type_mapping.keys()) geometry_type_values = np.array(list(type_mapping.values()), dtype=object) else: @@ -68,8 +72,11 @@ def _pygeos_to_shapely(geom): return None if compat.PYGEOS_SHAPELY_COMPAT: - geom = shapely.geos.lgeos.GEOSGeom_clone(geom._ptr) - return shapely.geometry.base.geom_factory(geom) + # we can only use this compatible fast path for shapely < 2, because + # shapely 2+ doesn't expose clone + if not compat.SHAPELY_GE_20: + geom = shapely.geos.lgeos.GEOSGeom_clone(geom._ptr) + return shapely.geometry.base.geom_factory(geom) # fallback going through WKB if pygeos.is_empty(geom) and pygeos.get_type_id(geom) == 0: @@ -161,11 +168,11 @@ def from_wkb(data): """ Convert a list or array of WKB objects to a np.ndarray[geoms]. """ + if compat.USE_SHAPELY_20: + return shapely.from_wkb(data) if compat.USE_PYGEOS: return pygeos.from_wkb(data) - import shapely.wkb - out = [] for geom in data: @@ -182,7 +189,9 @@ def from_wkb(data): def to_wkb(data, hex=False, **kwargs): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.to_wkb(data, hex=hex, **kwargs) + elif compat.USE_PYGEOS: return pygeos.to_wkb(data, hex=hex, **kwargs) else: if hex: @@ -196,11 +205,11 @@ def from_wkt(data): """ Convert a list or array of WKT objects to a np.ndarray[geoms]. """ + if compat.USE_SHAPELY_20: + return shapely.from_wkt(data) if compat.USE_PYGEOS: return pygeos.from_wkt(data) - import shapely.wkt - out = [] for geom in data: @@ -219,7 +228,9 @@ def from_wkt(data): def to_wkt(data, **kwargs): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.to_wkt(data, **kwargs) + elif compat.USE_PYGEOS: return pygeos.to_wkt(data, **kwargs) else: out = [geom.wkt if geom is not None else None for geom in data] @@ -246,7 +257,9 @@ def points_from_xy(x, y, z=None): if z is not None: z = np.asarray(z, dtype="float64") - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.points(x, y, z) + elif compat.USE_PYGEOS: return pygeos.points(x, y, z) else: out = _points_from_xy(x, y, z) @@ -464,21 +477,27 @@ def _unary_op(op, left, null_value=False): def is_valid(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.is_valid(data) + elif compat.USE_PYGEOS: return pygeos.is_valid(data) else: return _unary_op("is_valid", data, null_value=False) def is_empty(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.is_empty(data) + elif compat.USE_PYGEOS: return pygeos.is_empty(data) else: return _unary_op("is_empty", data, null_value=False) def is_simple(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.is_simple(data) + elif compat.USE_PYGEOS: return pygeos.is_simple(data) else: return _unary_op("is_simple", data, null_value=False) @@ -510,21 +529,28 @@ def is_ring(data): def is_closed(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.is_closed(data) + elif compat.USE_PYGEOS: return pygeos.is_closed(data) else: return _unary_op("is_closed", data, null_value=False) def has_z(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.has_z(data) + elif compat.USE_PYGEOS: return pygeos.has_z(data) else: return _unary_op("has_z", data, null_value=False) def geom_type(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + res = shapely.get_type_id(data) + return geometry_type_values[np.searchsorted(geometry_type_ids, res)] + elif compat.USE_PYGEOS: res = pygeos.get_type_id(data) return geometry_type_values[np.searchsorted(geometry_type_ids, res)] else: @@ -532,14 +558,18 @@ def geom_type(data): def area(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.area(data) + elif compat.USE_PYGEOS: return pygeos.area(data) else: return _unary_op("area", data, null_value=np.nan) def length(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.length(data) + elif compat.USE_PYGEOS: return pygeos.length(data) else: return _unary_op("length", data, null_value=np.nan) @@ -561,35 +591,45 @@ def _unary_geo(op, left, *args, **kwargs): def boundary(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.boundary(data) + elif compat.USE_PYGEOS: return pygeos.boundary(data) else: return _unary_geo("boundary", data) def centroid(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.centroid(data) + elif compat.USE_PYGEOS: return pygeos.centroid(data) else: return _unary_geo("centroid", data) def convex_hull(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.convex_hull(data) + elif compat.USE_PYGEOS: return pygeos.convex_hull(data) else: return _unary_geo("convex_hull", data) def envelope(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.envelope(data) + elif compat.USE_PYGEOS: return pygeos.envelope(data) else: return _unary_geo("envelope", data) def exterior(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.get_exterior_ring(data) + elif compat.USE_PYGEOS: return pygeos.get_exterior_ring(data) else: return _unary_geo("exterior", data) @@ -639,14 +679,18 @@ def representative_point(data): def covers(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.covers(data, other) + elif compat.USE_PYGEOS: return _binary_method("covers", data, other) else: return _binary_predicate("covers", data, other) def covered_by(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.covered_by(data, other) + elif compat.USE_PYGEOS: return _binary_method("covered_by", data, other) else: raise NotImplementedError( @@ -655,70 +699,88 @@ def covered_by(data, other): def contains(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.contains(data, other) + elif compat.USE_PYGEOS: return _binary_method("contains", data, other) else: return _binary_predicate("contains", data, other) def crosses(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.crosses(data, other) + elif compat.USE_PYGEOS: return _binary_method("crosses", data, other) else: return _binary_predicate("crosses", data, other) def disjoint(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.disjoint(data, other) + elif compat.USE_PYGEOS: return _binary_method("disjoint", data, other) else: return _binary_predicate("disjoint", data, other) def equals(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.equals(data, other) + elif compat.USE_PYGEOS: return _binary_method("equals", data, other) else: return _binary_predicate("equals", data, other) def intersects(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.intersects(data, other) + elif compat.USE_PYGEOS: return _binary_method("intersects", data, other) else: return _binary_predicate("intersects", data, other) def overlaps(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.overlaps(data, other) + elif compat.USE_PYGEOS: return _binary_method("overlaps", data, other) else: return _binary_predicate("overlaps", data, other) def touches(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.touches(data, other) + elif compat.USE_PYGEOS: return _binary_method("touches", data, other) else: return _binary_predicate("touches", data, other) def within(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.within(data, other) + elif compat.USE_PYGEOS: return _binary_method("within", data, other) else: return _binary_predicate("within", data, other) def equals_exact(data, other, tolerance): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.equals_exact(data, other, tolerance=tolerance) + elif compat.USE_PYGEOS: return _binary_method("equals_exact", data, other, tolerance=tolerance) else: return _binary_predicate("equals_exact", data, other, tolerance=tolerance) def almost_equals(self, other, decimal): - if compat.USE_PYGEOS: + if compat.USE_PYGEOS or compat.USE_SHAPELY_20: return self.equals_exact(other, 0.5 * 10 ** (-decimal)) else: return _binary_predicate("almost_equals", self, other, decimal=decimal) @@ -744,28 +806,36 @@ def clip_by_rect(data, xmin, ymin, xmax, ymax): def difference(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.difference(data, other) + elif compat.USE_PYGEOS: return _binary_method("difference", data, other) else: return _binary_geo("difference", data, other) def intersection(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.intersection(data, other) + elif compat.USE_PYGEOS: return _binary_method("intersection", data, other) else: return _binary_geo("intersection", data, other) def symmetric_difference(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.symmetric_difference(data, other) + elif compat.USE_PYGEOS: return _binary_method("symmetric_difference", data, other) else: return _binary_geo("symmetric_difference", data, other) def union(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.union(data, other) + elif compat.USE_PYGEOS: return _binary_method("union", data, other) else: return _binary_geo("union", data, other) @@ -777,14 +847,23 @@ def union(data, other): def distance(data, other): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.distance(data, other) + elif compat.USE_PYGEOS: return _binary_method("distance", data, other) else: return _binary_op_float("distance", data, other) def buffer(data, distance, resolution=16, **kwargs): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + if compat.SHAPELY_G_20a1: + return shapely.buffer(data, distance, quad_segs=resolution, **kwargs) + else: + # TODO: temporary keep this (so geopandas works with latest released + # shapely, currently alpha1) until shapely beta1 is out + return shapely.buffer(data, distance, quadsegs=resolution, **kwargs) + elif compat.USE_PYGEOS: return pygeos.buffer(data, distance, quadsegs=resolution, **kwargs) else: out = np.empty(len(data), dtype=object) @@ -815,7 +894,9 @@ def buffer(data, distance, resolution=16, **kwargs): def interpolate(data, distance, normalized=False): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.line_interpolate_point(data, distance, normalized=normalized) + elif compat.USE_PYGEOS: try: return pygeos.line_interpolate_point(data, distance, normalized=normalized) except TypeError: # support for pygeos<0.9 @@ -843,7 +924,9 @@ def interpolate(data, distance, normalized=False): def simplify(data, tolerance, preserve_topology=True): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.simplify(data, tolerance, preserve_topology=preserve_topology) + elif compat.USE_PYGEOS: # preserve_topology has different default as pygeos! return pygeos.simplify(data, tolerance, preserve_topology=preserve_topology) else: @@ -874,7 +957,9 @@ def _shapely_normalize(geom): def normalize(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.normalize(data) + elif compat.USE_PYGEOS: return pygeos.normalize(data) else: out = np.empty(len(data), dtype=object) @@ -886,7 +971,9 @@ def normalize(data): def project(data, other, normalized=False): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.line_locate_point(data, other, normalized=normalized) + elif compat.USE_PYGEOS: try: return pygeos.line_locate_point(data, other, normalized=normalized) except TypeError: # support for pygeos<0.9 @@ -896,6 +983,8 @@ def project(data, other, normalized=False): def relate(data, other): + if compat.USE_SHAPELY_20: + return shapely.relate(data, other) data = to_shapely(data) if isinstance(other, np.ndarray): other = to_shapely(other) @@ -903,7 +992,9 @@ def relate(data, other): def unary_union(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.union_all(data) + elif compat.USE_PYGEOS: return _pygeos_to_shapely(pygeos.union_all(data)) else: data = [g for g in data if g is not None] @@ -919,21 +1010,27 @@ def unary_union(data): def get_x(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.get_x(data) + elif compat.USE_PYGEOS: return pygeos.get_x(data) else: return _unary_op("x", data, null_value=np.nan) def get_y(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.get_y(data) + elif compat.USE_PYGEOS: return pygeos.get_y(data) else: return _unary_op("y", data, null_value=np.nan) def get_z(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.get_z(data) + elif compat.USE_PYGEOS: return pygeos.get_z(data) else: data = [geom.z if geom.has_z else np.nan for geom in data] @@ -941,7 +1038,9 @@ def get_z(data): def bounds(data): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.bounds(data) + elif compat.USE_PYGEOS: return pygeos.bounds(data) # ensure that for empty arrays, the result has the correct shape if len(data) == 0: @@ -965,6 +1064,11 @@ def bounds(data): def transform(data, func): + if compat.USE_SHAPELY_20: + coords = shapely.get_coordinates(data) + new_coords = func(coords[:, 0], coords[:, 1]) + result = shapely.set_coordinates(data.copy(), np.array(new_coords).T) + return result if compat.USE_PYGEOS: coords = pygeos.get_coordinates(data) new_coords = func(coords[:, 0], coords[:, 1]) diff --git a/geopandas/array.py b/geopandas/array.py index 9b47d6a..e30e4a8 100644 --- a/geopandas/array.py +++ b/geopandas/array.py @@ -108,7 +108,9 @@ def _geom_to_shapely(geom): """ Convert internal representation (PyGEOS or Shapely) to external Shapely object. """ - if not compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return geom + elif not compat.USE_PYGEOS: return geom else: return vectorized._pygeos_to_shapely(geom) @@ -118,7 +120,9 @@ def _shapely_to_geom(geom): """ Convert external Shapely object to internal representation (PyGEOS or Shapely). """ - if not compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return geom + elif not compat.USE_PYGEOS: return geom else: return vectorized._shapely_to_pygeos(geom) @@ -410,7 +414,9 @@ class GeometryArray(ExtensionArray): # ) def __getstate__(self): - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return (shapely.to_wkb(self.data), self._crs) + elif compat.USE_PYGEOS: return (pygeos.to_wkb(self.data), self._crs) else: return self.__dict__ @@ -1065,7 +1071,9 @@ class GeometryArray(ExtensionArray): """ Boolean NumPy array indicating if each value is missing """ - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + return shapely.is_missing(self.data) + elif compat.USE_PYGEOS: return pygeos.is_missing(self.data) else: return np.array([g is None for g in self.data], dtype="bool") diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py index 1747cb8..8dc054e 100644 --- a/geopandas/geoseries.py +++ b/geopandas/geoseries.py @@ -7,6 +7,7 @@ from pandas import Series, MultiIndex, DataFrame from pandas.core.internals import SingleBlockManager from pyproj import CRS +import shapely from shapely.geometry.base import BaseGeometry from geopandas.base import GeoPandasBase, _delegate_property @@ -871,12 +872,17 @@ class GeoSeries(GeoPandasBase, Series): ) index_parts = True - if compat.USE_PYGEOS and compat.PYGEOS_GE_09: - import pygeos # noqa + if compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_09): + if compat.USE_SHAPELY_20: + geometries, outer_idx = shapely.get_parts( + self.values.data, return_index=True + ) + else: + import pygeos # noqa - geometries, outer_idx = pygeos.get_parts( - self.values.data, return_index=True - ) + geometries, outer_idx = pygeos.get_parts( + self.values.data, return_index=True + ) if len(outer_idx): # Generate inner index as a range per value of outer_idx diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py index 93b8c18..1e8db4d 100644 --- a/geopandas/io/sql.py +++ b/geopandas/io/sql.py @@ -3,6 +3,7 @@ from contextlib import contextmanager import pandas as pd +import shapely import shapely.wkb from geopandas import GeoDataFrame @@ -85,7 +86,10 @@ def _df_to_geodf(df, geom_col="geom", crs=None): df[geom_col] = geoms = geoms.apply(load_geom) if crs is None: - srid = shapely.geos.lgeos.GEOSGetSRID(geoms.iat[0]._geom) + if compat.SHAPELY_GE_20: + srid = shapely.get_srid(geoms.iat[0]) + else: + srid = shapely.geos.lgeos.GEOSGetSRID(geoms.iat[0]._geom) # if no defined SRID in geodatabase, returns SRID of 0 if srid != 0: crs = "epsg:{}".format(srid) @@ -284,7 +288,14 @@ def _convert_linearring_to_linestring(gdf, geom_name): def _convert_to_ewkb(gdf, geom_name, srid): """Convert geometries to ewkb.""" - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20: + geoms = shapely.to_wkb( + shapely.set_srid(gdf[geom_name].values.data, srid=srid), + hex=True, + include_srid=True, + ) + + elif compat.USE_PYGEOS: from pygeos import set_srid, to_wkb geoms = to_wkb( diff --git a/geopandas/io/tests/test_pickle.py b/geopandas/io/tests/test_pickle.py index 85bbcba..5e3dbc0 100644 --- a/geopandas/io/tests/test_pickle.py +++ b/geopandas/io/tests/test_pickle.py @@ -48,9 +48,12 @@ def with_use_pygeos(option): @pytest.mark.skipif( - compat.USE_PYGEOS or (Version(pyproj.__version__) < Version("2.4")), + compat.USE_SHAPELY_20 + or compat.USE_PYGEOS + or (Version(pyproj.__version__) < Version("2.4")), reason=( - "pygeos-based unpickling currently only works for pygeos-written files; " + "shapely 2.0/pygeos-based unpickling currently only works for " + "shapely-2.0/pygeos-written files; " "old pyproj versions can't read pickles from newer pyproj versions" ), ) diff --git a/geopandas/sindex.py b/geopandas/sindex.py index 85044ca..b37755b 100644 --- a/geopandas/sindex.py +++ b/geopandas/sindex.py @@ -14,7 +14,7 @@ def _get_sindex_class(): Required to comply with _compat.USE_PYGEOS. The selection order goes PyGEOS > RTree > Error. """ - if compat.USE_PYGEOS: + if compat.USE_SHAPELY_20 or compat.USE_PYGEOS: return PyGEOSSTRTreeIndex if compat.HAS_RTREE: return RTreeIndex @@ -627,13 +627,19 @@ if compat.HAS_RTREE: return self.size -if compat.HAS_PYGEOS: +if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS: from . import geoseries # noqa from . import array # noqa - import pygeos # noqa - _PYGEOS_PREDICATES = {p.name for p in pygeos.strtree.BinaryPredicate} | set([None]) + if compat.USE_SHAPELY_20: + import shapely as mod # noqa + + _PYGEOS_PREDICATES = {p.name for p in mod.strtree.BinaryPredicate} | set([None]) + else: + import pygeos as mod # noqa + + _PYGEOS_PREDICATES = {p.name for p in mod.strtree.BinaryPredicate} | set([None]) class PyGEOSSTRTreeIndex(BaseSpatialIndex): """A simple wrapper around pygeos's STRTree. @@ -651,9 +657,9 @@ if compat.HAS_PYGEOS: # https://github.com/pygeos/pygeos/issues/146 # https://github.com/pygeos/pygeos/issues/147 non_empty = geometry.copy() - non_empty[pygeos.is_empty(non_empty)] = None + non_empty[mod.is_empty(non_empty)] = None # set empty geometries to None to maintain indexing - self._tree = pygeos.STRtree(non_empty) + self._tree = mod.STRtree(non_empty) # store geometries, including empty geometries for user access self.geometries = geometry.copy() @@ -719,6 +725,8 @@ if compat.HAS_PYGEOS: return geometry.data elif isinstance(geometry, BaseGeometry): return array._shapely_to_geom(geometry) + elif geometry is None: + return None elif isinstance(geometry, list): return np.asarray( [ @@ -742,7 +750,10 @@ if compat.HAS_PYGEOS: geometry = self._as_geometry_array(geometry) - res = self._tree.query_bulk(geometry, predicate) + if compat.USE_SHAPELY_20: + res = self._tree.query(geometry, predicate) + else: + res = self._tree.query_bulk(geometry, predicate) if sort: # sort by first array (geometry) and then second (tree) @@ -756,23 +767,34 @@ if compat.HAS_PYGEOS: def nearest( self, geometry, return_all=True, max_distance=None, return_distance=False ): - if not compat.PYGEOS_GE_010: - raise NotImplementedError("sindex.nearest requires pygeos >= 0.10") + if not (compat.USE_SHAPELY_20 or compat.PYGEOS_GE_010): + raise NotImplementedError( + "sindex.nearest requires shapely >= 2.0 or pygeos >= 0.10" + ) geometry = self._as_geometry_array(geometry) + if isinstance(geometry, BaseGeometry) or geometry is None: + geometry = [geometry] - if not return_all and max_distance is None and not return_distance: - return self._tree.nearest(geometry) - - result = self._tree.nearest_all( - geometry, max_distance=max_distance, return_distance=return_distance - ) + if compat.USE_SHAPELY_20: + result = self._tree.query_nearest( + geometry, + max_distance=max_distance, + return_distance=return_distance, + all_matches=return_all, + ) + else: + if not return_all and max_distance is None and not return_distance: + return self._tree.nearest(geometry) + result = self._tree.nearest_all( + geometry, max_distance=max_distance, return_distance=return_distance + ) if return_distance: indices, distances = result else: indices = result - if not return_all: + if not return_all and not compat.USE_SHAPELY_20: # first subarray of geometry indices is sorted, so we can use this # trick to get the first of each index value mask = np.diff(indices[0, :]).astype("bool") @@ -806,9 +828,9 @@ if compat.HAS_PYGEOS: # need to convert tuple of bounds to a geometry object if len(coordinates) == 4: - indexes = self._tree.query(pygeos.box(*coordinates)) + indexes = self._tree.query(mod.box(*coordinates)) elif len(coordinates) == 2: - indexes = self._tree.query(pygeos.points(*coordinates)) + indexes = self._tree.query(mod.points(*coordinates)) else: raise TypeError( "Invalid coordinates, must be iterable in format " diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py index 00750ec..cbbde3e 100644 --- a/geopandas/tests/test_array.py +++ b/geopandas/tests/test_array.py @@ -11,7 +11,11 @@ import shapely.geometry from shapely.geometry.base import CAP_STYLE, JOIN_STYLE import shapely.wkb import shapely.wkt -from shapely._buildcfg import geos_version + +try: + from shapely import geos_version +except ImportError: + from shapely._buildcfg import geos_version import geopandas from geopandas.array import ( @@ -142,7 +146,7 @@ def test_from_wkb(): # missing values # TODO(pygeos) does not support empty strings, np.nan, or pd.NA missing_values = [None] - if not compat.USE_PYGEOS: + if not (compat.USE_SHAPELY_20 or compat.USE_PYGEOS): missing_values.extend([b"", np.nan]) missing_values.append(pd.NA) @@ -216,7 +220,7 @@ def test_from_wkt(string_type): # missing values # TODO(pygeos) does not support empty strings, np.nan, or pd.NA missing_values = [None] - if not compat.USE_PYGEOS: + if not (compat.USE_SHAPELY_20 or compat.USE_PYGEOS): missing_values.extend([f(""), np.nan]) missing_values.append(pd.NA) diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py index 3cff510..db63944 100644 --- a/geopandas/tests/test_geodataframe.py +++ b/geopandas/tests/test_geodataframe.py @@ -23,7 +23,7 @@ from pandas.testing import assert_frame_equal, assert_index_equal, assert_series import pytest -TEST_NEAREST = compat.PYGEOS_GE_010 and compat.USE_PYGEOS +TEST_NEAREST = compat.USE_SHAPELY_20 or (compat.PYGEOS_GE_010 and compat.USE_PYGEOS) pandas_133 = Version(pd.__version__) == Version("1.3.3") @@ -895,7 +895,7 @@ class TestDataFrame: @pytest.mark.parametrize("how", ["left", "inner", "right"]) @pytest.mark.parametrize("predicate", ["intersects", "within", "contains"]) @pytest.mark.skipif( - not compat.USE_PYGEOS and not compat.HAS_RTREE, + not (compat.USE_PYGEOS and compat.HAS_RTREE and compat.USE_SHAPELY_20), reason="sjoin needs `rtree` or `pygeos` dependency", ) def test_sjoin(self, how, predicate): diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py index fdd0e3e..dbefd4b 100644 --- a/geopandas/tests/test_geom_methods.py +++ b/geopandas/tests/test_geom_methods.py @@ -591,7 +591,7 @@ class TestGeomMethods: assert_series_equal(res, exp) @pytest.mark.skipif( - not compat.USE_PYGEOS, + not (compat.USE_PYGEOS or compat.USE_SHAPELY_20), reason="covered_by is only implemented for pygeos, not shapely", ) def test_covered_by(self): diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py index c36c77d..eaa2122 100644 --- a/geopandas/tests/test_sindex.py +++ b/geopandas/tests/test_sindex.py @@ -17,8 +17,10 @@ from geopandas import GeoDataFrame, GeoSeries, read_file, datasets import pytest import numpy as np -if compat.USE_PYGEOS: - import pygeos +if compat.USE_SHAPELY_20: + import shapely as mod +elif compat.USE_PYGEOS: + import pygeos as mod @pytest.mark.skip_no_sindex @@ -698,7 +700,7 @@ class TestPygeosInterface: # ------------------------- `nearest` tests ------------------------- # @pytest.mark.skipif( - compat.USE_PYGEOS, + compat.USE_PYGEOS or compat.USE_SHAPELY_20, reason=("RTree supports sindex.nearest with different behaviour"), ) def test_rtree_nearest_warns(self): @@ -709,7 +711,7 @@ class TestPygeosInterface: df.sindex.nearest((0, 0, 1, 1), num_results=2) @pytest.mark.skipif( - not (compat.USE_PYGEOS and not compat.PYGEOS_GE_010), + compat.USE_SHAPELY_20 or not (compat.USE_PYGEOS and not compat.PYGEOS_GE_010), reason=("PyGEOS < 0.10 does not support sindex.nearest"), ) def test_pygeos_error(self): @@ -718,7 +720,7 @@ class TestPygeosInterface: df.sindex.nearest(None) @pytest.mark.skipif( - not (compat.USE_PYGEOS and compat.PYGEOS_GE_010), + not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)), reason=("PyGEOS >= 0.10 is required to test sindex.nearest"), ) @pytest.mark.parametrize("return_all", [True, False]) @@ -730,19 +732,19 @@ class TestPygeosInterface: ], ) def test_nearest_single(self, geometry, expected, return_all): - geoms = pygeos.points(np.arange(10), np.arange(10)) + geoms = mod.points(np.arange(10), np.arange(10)) df = geopandas.GeoDataFrame({"geometry": geoms}) p = Point(geometry) res = df.sindex.nearest(p, return_all=return_all) assert_array_equal(res, expected) - p = pygeos.points(geometry) + p = mod.points(geometry) res = df.sindex.nearest(p, return_all=return_all) assert_array_equal(res, expected) @pytest.mark.skipif( - not compat.USE_PYGEOS or not compat.PYGEOS_GE_010, + not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)), reason=("PyGEOS >= 0.10 is required to test sindex.nearest"), ) @pytest.mark.parametrize("return_all", [True, False]) @@ -754,14 +756,14 @@ class TestPygeosInterface: ], ) def test_nearest_multi(self, geometry, expected, return_all): - geoms = pygeos.points(np.arange(10), np.arange(10)) + geoms = mod.points(np.arange(10), np.arange(10)) df = geopandas.GeoDataFrame({"geometry": geoms}) ps = [Point(p) for p in geometry] res = df.sindex.nearest(ps, return_all=return_all) assert_array_equal(res, expected) - ps = pygeos.points(geometry) + ps = mod.points(geometry) res = df.sindex.nearest(ps, return_all=return_all) assert_array_equal(res, expected) @@ -775,7 +777,7 @@ class TestPygeosInterface: assert_array_equal(res, expected) @pytest.mark.skipif( - not compat.USE_PYGEOS or not compat.PYGEOS_GE_010, + not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)), reason=("PyGEOS >= 0.10 is required to test sindex.nearest"), ) @pytest.mark.parametrize("return_all", [True, False]) @@ -787,14 +789,14 @@ class TestPygeosInterface: ], ) def test_nearest_none(self, geometry, expected, return_all): - geoms = pygeos.points(np.arange(10), np.arange(10)) + geoms = mod.points(np.arange(10), np.arange(10)) df = geopandas.GeoDataFrame({"geometry": geoms}) res = df.sindex.nearest(geometry, return_all=return_all) assert_array_equal(res, expected) @pytest.mark.skipif( - not compat.USE_PYGEOS or not compat.PYGEOS_GE_010, + not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)), reason=("PyGEOS >= 0.10 is required to test sindex.nearest"), ) @pytest.mark.parametrize("return_distance", [True, False]) @@ -810,7 +812,7 @@ class TestPygeosInterface: def test_nearest_max_distance( self, expected, max_distance, return_all, return_distance ): - geoms = pygeos.points(np.arange(10), np.arange(10)) + geoms = mod.points(np.arange(10), np.arange(10)) df = geopandas.GeoDataFrame({"geometry": geoms}) ps = [Point(0.5, 0.5), Point(0, 10)] diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index c693ca1..19c7edf 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -360,10 +360,10 @@ def _nearest_query( how: str, return_distance: bool, ): - if not (compat.PYGEOS_GE_010 and compat.USE_PYGEOS): + if not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)): raise NotImplementedError( - "Currently, only PyGEOS >= 0.10.0 supports `nearest_all`. " - + compat.INSTALL_PYGEOS_ERROR + "Currently, only PyGEOS >= 0.10.0 or Shapely >= 2.0 supports " + "`nearest_all`. " + compat.INSTALL_PYGEOS_ERROR ) # use the opposite of the join direction for the index use_left_as_sindex = how == "right" diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py index 50bdfac..953a1fa 100644 --- a/geopandas/tools/tests/test_sjoin.py +++ b/geopandas/tools/tests/test_sjoin.py @@ -17,7 +17,7 @@ from pandas.testing import assert_frame_equal import pytest -TEST_NEAREST = compat.PYGEOS_GE_010 and compat.USE_PYGEOS +TEST_NEAREST = compat.USE_SHAPELY_20 or (compat.PYGEOS_GE_010 and compat.USE_PYGEOS) pytestmark = pytest.mark.skip_no_sindex @@ -583,7 +583,7 @@ def test_no_nearest_all(): df2 = geopandas.GeoDataFrame({"geometry": []}) with pytest.raises( NotImplementedError, - match="Currently, only PyGEOS >= 0.10.0 supports `nearest_all`", + match="Currently, only PyGEOS >= 0.10.0 or Shapely >= 2.0 supports", ): sjoin_nearest(df1, df2)