Make GeoPandas work with Shapely 2.0 (#2275)

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

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

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::
 

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

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

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

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

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(

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"
     ),
 )

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 "

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)
 

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

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

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

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"

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)