Repository: geopandas/geopandas
Pull request #3163 (2024-01-27 12:05:39)

MAINT: Add optional-dependencies to pyproject.toml
## What

The specification for `optional-dependencies` is part of PEP 621
<https://peps.python.org/pep-0621/>. See
<https://packaging.python.org/en/latest/specifications/pyproject-toml/#pyproject-toml-spec>
for latest specification.

I added the dependencies to `optional-dependencies` based on
`./requirements-dev.txt`, `./environment.yml` and #2990.

## Why?

This would be useful for defining the optional dependencies when packaging
`geopandas` for e.g. `nixpkgs`
(<https://github.com/NixOS/nixpkgs/blob/master/pkgs/development/python-modules/geopandas/default.nix>)
as it is better for them to be defined here in the main repository than
downstream. As there is already code in `./geopandas/` that use these optional
dependencies, I believe they should be documented somewhere.

## Misc

I did not update `requirements-dev.txt` as I do not know if it is currently
up-to-date (<https://github.com/geopandas/geopandas/pull/2438#issuecomment-1205268708>)
or being maintained.

Marked as draft as no feedback yet given from maintainers. Hopefully this sounds good!


nialov: Took a look at some other packages and I am actually wondering if it would be simpler to just define two optional dependency groups: `all` and `dev`. `all` includes all dependencies except the development ones. `dev` contains `pytest` and others needed for testing and development.`all` is used by well-known packages, such as `pandas` and `typer`, to indicate all optional dependencies. 

If `all` could be defined by pointing to the already defined groups (`postgis`, `plot`, etc.), that would be the best way to define `all`. But that is not possible without introducing a possible circular dependency issue. E.g. `pandas` defines `all` by redefining all the dependencies of all other groups in `all` (<https://github.com/pandas-dev/pandas/blob/56b5979f136e72ce78e5221392739481c025d5e7/pyproject.toml#L88>). This is not nice to maintain.

The grouping feature is not apparent to users and writing e.g. `pip install geopandas[postgis,postgis-write,plot]` is demanding a lot from any casual installer of `geopandas`. I would also wager that either a) an user wants to minimize dependencies, in which case they can add the optional dependencies as they run into errors as they run `geopandas` code or b) just want all functionality of `geopandas` quickly. For a) the groupings are not much better than just manually installing the required packages. For b) the `all` group gives a simple way to install all optionals.
martinfleis: I just want to say that I believe it is nice to have optional dependencies defined in some way in pyproject. However, I don't think I've ever used it a as a conda/mamba user myself so don't really care if it is split in groups of just `all` and `dev`. The latter seems more maintainable.
m-richards: I also don't have a strong leaning on all & dev vs split out. I use pip rather than conda for work, but I don't tend to use extras, I'll just install the base package. Even if I know and extra exists I usually have to double check what it's called  - and the most reliable way of doing that is looking at the pyproject.toml, at which point I no longer need the extra because I can see the packages underneath it. 
nialov: Updated the pull request to use the two groups discussed above, `all` and `dev`. Also updated documentation to inform users and developers of these options. One change, that I do not know if is appropriate, is the use of `psycopg2-binary` instead of `psycopg2`. Using this package is easier as it contains binary components of `psycopg2` directly (<https://github.com/psycopg/psycopg2#installation>).

Patch:
diff --git a/.github/workflows/versioneer.yml b/.github/workflows/versioneer.yml
index 1605edd6..38283c44 100644
--- a/.github/workflows/versioneer.yml
+++ b/.github/workflows/versioneer.yml
@@ -32,7 +32,7 @@ jobs:
           git reset -- geopandas/__init__.py
           git checkout -- geopandas/__init__.py
       - name: Create PR
-        uses: peter-evans/create-pull-request@v5
+        uses: peter-evans/create-pull-request@v6
         with:
           title: "Update Versioneer"
           branch: update-versioneer
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 9d3e00f2..cf293cb5 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -5,7 +5,7 @@ ci:
 files: 'geopandas\/'
 repos:
     - repo: https://github.com/psf/black
-      rev: 23.1.0
+      rev: 24.2.0
       hooks:
           - id: black
             language_version: python3
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ef544f6c..b0a18651 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,9 +11,13 @@ Notes on dependencies:
 API changes:
 
 - `unary_union` is now deprecated and replaced by the `union_all` method (#3007).
+- The `sjoin` method will now preserve the name of the index of the right
+  GeoDataFrame, if it has one, instead of always using `"index_right"` as the
+  name for the resulting column in the return value (#846, #2144).
 
 New methods:
 
+- Added `set_precision` and `get_precision` methods from shapely to GeoSeries/GeoDataframe (#3175).
 - Added `count_coordinates` method from shapely to GeoSeries/GeoDataframe (#3026).
 - Added `minimum_clearance` method from shapely to GeoSeries/GeoDataframe (#2989).
 - Added `is_ccw` method from shapely to GeoSeries/GeoDataframe (#3027).
@@ -22,6 +26,8 @@ New methods:
 - Added `contains_properly` method from shapely to GeoSeries/GeoDataframe (#3105).
 - Added `snap` method from shapely to GeoSeries/GeoDataframe (#3086).
 - Added `transform` method from shapely to GeoSeries/GeoDataFrame (#3075).
+- Added `dwithin` method to check for a "distance within" predicate on
+  GeoSeries/GeoDataFrame (#3153).
 - Added `to_geo_dict` method to generate GeoJSON-like dictionary from a GeoDataFrame (#3132).
 
 New features and improvements:
@@ -35,6 +41,13 @@ New features and improvements:
   `buffer(0)` (#3113).
 - Passing `"geometry"` as `dtype` to `pd.read_csv` will now return a GeoSeries for
   the specified columns (#3101).
+- 
+API changes:
+- Added support for ``mask`` keyword for pyogrio engine for pyogrio >= 0.7.0 (#3062).
+
+Backwards incompatible API changes:
+- The deprecated default value of GeoDataFrame/ GeoSeries `explode(.., index_parts=True)` is now
+  set to false for consistency with pandas (#3174)
 
 Potentially breaking changes:
 
@@ -42,9 +55,6 @@ Potentially breaking changes:
   now returns a Pandas DataFrame instead of a GeoDataFrame with an empty
   ``geometry`` column.
 
-API changes:
-- Added support for ``mask`` keyword for pyogrio engine for pyogrio >= 0.7.0 (#3062).
-
 Bug fixes:
 
 - Fix `GeoDataFrame.merge()` incorrectly returning a `DataFrame` instead of a
@@ -56,8 +66,20 @@ Deprecations and compatibility notes:
 - The deprecation of `geopandas.datasets` has been enforced and the module has been
   removed. New sample datasets are now available in the
   [geodatasets](https://geodatasets.readthedocs.io/en/latest/) package (#3084).
-
-## Version 0.14.3 (Jan ??, 2024)
+- Many longstanding deprecated functions, methods and properties have been removed (#3174)
+  - Removed deprecated functions  
+    `geopandas.io.read_file`, `geopandas.io.to_file` and `geopandas.io.sql.read_postgis`. 
+    `geopandas.read_file`, `geopandas.read_postgis` and the GeoDataFrame/GeoSeries `to_file(..)` 
+    method should be used instead.
+  - Removed deprecated `GeometryArray.data` property, `np.asarray(..)` or the `to_numpy()`
+    method should be used instead.
+  - Removed deprecated `sindex.query_bulk` method, using `sindex.query` instead.
+  - Removed deprecated `sjoin` parameter `op`, `predicate` should be supplied instead.
+  - Removed deprecated GeoSeries/ GeoDataFrame methods `__xor__`, `__or__`, `__and__` and 
+    `__sub__`. Instead use methods `symmetric_difference`, `union`, `intersection` and 
+    `difference` respectively.
+
+## Version 0.14.3 (Jan 31, 2024)
 
 - Several fixes for compatibility with the latest pandas 2.2 release.
 - Fix bug in `pandas.concat` CRS consistency checking where CRS differing by WKT
diff --git a/ci/envs/310-pd20-conda-forge.yaml b/ci/envs/310-pd20-conda-forge.yaml
index 91e55bc9..4854fdb8 100644
--- a/ci/envs/310-pd20-conda-forge.yaml
+++ b/ci/envs/310-pd20-conda-forge.yaml
@@ -11,7 +11,8 @@ dependencies:
   - pyproj
   - packaging
   # testing
-  - pytest
+  # TEMP pin, 8.1 breaks doctestplus, see https://github.com/scientific-python/pytest-doctestplus/issues/239
+  - pytest==8.0.*
   - pytest-cov
   - pytest-xdist
   - fsspec
diff --git a/ci/envs/311-dev.yaml b/ci/envs/311-dev.yaml
index bc7964a3..93c8c534 100644
--- a/ci/envs/311-dev.yaml
+++ b/ci/envs/311-dev.yaml
@@ -24,7 +24,8 @@ dependencies:
     - mapclassify>=2.4.0
     # dev versions of packages
     - --pre --index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple --extra-index-url https://pypi.fury.io/arrow-nightlies/ --extra-index-url https://pypi.org/simple
-    - fiona
+    # TEMP pin, see https://github.com/geopandas/geopandas/issues/3207
+    - fiona==1.9.*
     - pandas
     - matplotlib
     - pyarrow
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 9beb9109..f8ddd866 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -107,7 +107,7 @@ master_doc = "index"
 
 # General information about the project.
 project = "GeoPandas"
-copyright = "2013–2022, GeoPandas developers"
+copyright = "2013–, GeoPandas developers"
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index d0f08e9f..6a8fbbee 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -35,6 +35,8 @@ General methods and attributes
    GeoSeries.z
    GeoSeries.get_coordinates
    GeoSeries.count_coordinates
+   GeoSeries.set_precision
+   GeoSeries.get_precision
 
 Unary predicates
 ----------------
@@ -61,6 +63,7 @@ Binary predicates
    GeoSeries.contains_properly
    GeoSeries.crosses
    GeoSeries.disjoint
+   GeoSeries.dwithin
    GeoSeries.geom_equals
    GeoSeries.geom_almost_equals
    GeoSeries.geom_equals_exact
diff --git a/doc/source/docs/user_guide/mergingdata.rst b/doc/source/docs/user_guide/mergingdata.rst
index 85627e3c..b2be6314 100644
--- a/doc/source/docs/user_guide/mergingdata.rst
+++ b/doc/source/docs/user_guide/mergingdata.rst
@@ -39,7 +39,7 @@ In the following examples, these datasets are used:
 Appending
 ---------
 
-Appending :class:`GeoDataFrame` and :class:`GeoSeries` uses pandas :meth:`~pandas.DataFrame.append` methods.
+Appending :class:`GeoDataFrame` and :class:`GeoSeries` uses pandas :func:`~pandas.concat` function.
 Keep in mind, that appended geometry columns needs to have the same CRS.
 
 .. ipython:: python
diff --git a/doc/source/getting_started/install.rst b/doc/source/getting_started/install.rst
index 286564a1..2abb210b 100644
--- a/doc/source/getting_started/install.rst
+++ b/doc/source/getting_started/install.rst
@@ -101,6 +101,10 @@ as well::
     installing those. Using conda (see above) avoids the need to compile the
     dependencies yourself.
 
+Optional runtime dependencies can also be installed all at once::
+
+    pip install 'geopandas[all]'
+
 Installing from source
 ----------------------
 
@@ -111,6 +115,11 @@ You may install the latest development version by cloning the
     cd geopandas
     pip install .
 
+Development dependencies can be installed using the dev optional
+dependency group::
+
+    pip install '.[dev]'
+
 It is also possible to install the latest development version
 directly from the GitHub repository with::
 
diff --git a/doc/source/index.rst b/doc/source/index.rst
index a18559f4..d36504e8 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -1,3 +1,5 @@
+:html_theme.sidebar_secondary.remove:
+
 GeoPandas |version|
 ===================
 
@@ -72,10 +74,3 @@ travel, workshops, and a variety of other needs.
     :alt: numfocus
     :width: 400
     :target: https://numfocus.org/project/geopandas
-
-Indices and tables
-------------------
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index f254bd64..034d28fb 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -14,7 +14,8 @@ PANDAS_GE_14 = Version(pd.__version__) >= Version("1.4.0rc0")
 PANDAS_GE_15 = Version(pd.__version__) >= Version("1.5.0")
 PANDAS_GE_20 = Version(pd.__version__) >= Version("2.0.0")
 PANDAS_GE_21 = Version(pd.__version__) >= Version("2.1.0")
-PANDAS_GE_22 = Version(pd.__version__) >= Version("2.2.0.dev0")
+PANDAS_GE_22 = Version(pd.__version__) >= Version("2.2.0")
+PANDAS_GE_30 = Version(pd.__version__) >= Version("3.0.0.dev0")
 
 
 # -----------------------------------------------------------------------------
diff --git a/geopandas/_config.py b/geopandas/_config.py
index 7090395d..7815b606 100644
--- a/geopandas/_config.py
+++ b/geopandas/_config.py
@@ -5,6 +5,7 @@ Based on https://github.com/topper-123/optioneer, but simplified (don't deal
 with nested options, deprecated options, ..), just the attribute-style dict
 like holding the options and giving a nice repr.
 """
+
 from collections import namedtuple
 import textwrap
 
diff --git a/geopandas/_decorator.py b/geopandas/_decorator.py
index 6c263c75..6c3b2fe9 100644
--- a/geopandas/_decorator.py
+++ b/geopandas/_decorator.py
@@ -39,9 +39,11 @@ def doc(*docstrings: Union[str, Callable], **params) -> Callable:
 
         # formatting templates and concatenating docstring
         decorated.__doc__ = "".join(
-            component.format(**params)
-            if isinstance(component, str)
-            else dedent(component.__doc__ or "")
+            (
+                component.format(**params)
+                if isinstance(component, str)
+                else dedent(component.__doc__ or "")
+            )
             for component in docstring_components
         )
 
diff --git a/geopandas/array.py b/geopandas/array.py
index 7486b85d..ed71bb24 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -332,17 +332,6 @@ class GeometryArray(ExtensionArray):
         self.crs = crs
         self._sindex = None
 
-    @property
-    def data(self):
-        warnings.warn(
-            "Accessing the underlying geometries through the `.data` attribute is "
-            "deprecated and will be removed in GeoPandas 1.0. You can use "
-            "`np.asarray(..)` or the `to_numpy()` method instead.",
-            DeprecationWarning,
-            stacklevel=2,
-        )
-        return self._data
-
     @property
     def sindex(self):
         if self._sindex is None:
@@ -536,6 +525,9 @@ class GeometryArray(ExtensionArray):
         out[:] = [shapely.count_coordinates(s) for s in self._data]
         return out
 
+    def get_precision(self):
+        return shapely.get_precision(self._data)
+
     #
     # Unary operations that return new geometries
     #
@@ -657,6 +649,12 @@ class GeometryArray(ExtensionArray):
             shapely.transform(self._data, transformation, include_z), crs=self.crs
         )
 
+    def set_precision(self, grid_size, mode="valid_output"):
+        return GeometryArray(
+            shapely.set_precision(self._data, grid_size=grid_size, mode=mode),
+            crs=self.crs,
+        )
+
     #
     # Binary predicates
     #
@@ -708,6 +706,10 @@ class GeometryArray(ExtensionArray):
     def within(self, other):
         return self._binary_method("within", self, other)
 
+    def dwithin(self, other, distance):
+        self.check_geographic_crs(stacklevel=6)
+        return self._binary_method("dwithin", self, other, distance=distance)
+
     def geom_equals_exact(self, other, tolerance):
         return self._binary_method("equals_exact", self, other, tolerance=tolerance)
 
diff --git a/geopandas/base.py b/geopandas/base.py
index 053fda62..dfa59bc4 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -540,6 +540,56 @@ GeometryCollection
         """
         return _delegate_property("has_z", self)
 
+    def get_precision(self):
+        """Returns a ``Series`` of the precision of each geometry.
+
+        If a precision has not been previously set, it will be 0, indicating regular
+        double precision coordinates are in use. Otherwise, it will return the precision
+        grid size that was set on a geometry.
+
+        Returns NaN for not-a-geometry values.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Point(0, 1),
+        ...         Point(0, 1, 2),
+        ...         Point(0, 1.5, 2),
+        ...     ]
+        ... )
+        >>> s
+        0          POINT (0 1)
+        1      POINT Z (0 1 2)
+        2    POINT Z (0 1.5 2)
+        dtype: geometry
+
+        >>> s.get_precision()
+        0    0.0
+        1    0.0
+        2    0.0
+        dtype: float64
+
+        >>> s1 = s.set_precision(1)
+        >>> s1
+        0        POINT (0 1)
+        1    POINT Z (0 1 2)
+        2    POINT Z (0 2 2)
+        dtype: geometry
+
+        >>> s1.get_precision()
+        0    1.0
+        1    1.0
+        2    1.0
+        dtype: float64
+
+        See also
+        --------
+        GeoSeries.set_precision : set precision grid size
+        """
+        return Series(self.geometry.values.get_precision(), index=self.index)
+
     #
     # Unary operations that return a GeoSeries
     #
@@ -1056,6 +1106,95 @@ GeometryCollection
         """
         return _delegate_geo_method("remove_repeated_points", self, tolerance=tolerance)
 
+    def set_precision(self, grid_size, mode="valid_output"):
+        """Returns a ``GeoSeries`` with the precision set to a precision grid size.
+
+        By default, geometries use double precision coordinates (``grid_size=0``).
+
+        Coordinates will be rounded if a precision grid is less precise than the input
+        geometry. Duplicated vertices will be dropped from lines and polygons for grid
+        sizes greater than 0. Line and polygon geometries may collapse to empty
+        geometries if all vertices are closer together than ``grid_size``. Spikes or
+        sections in Polygons narrower than ``grid_size`` after rounding the vertices
+        will be removed, which can lead to MultiPolygons or empty geometries. Z values,
+        if present, will not be modified.
+
+        Parameters
+        ----------
+        grid_size : float
+            Precision grid size. If 0, will use double precision (will not modify
+            geometry if precision grid size was not previously set). If this value is
+            more precise than input geometry, the input geometry will not be modified.
+        mode : {'valid_output', 'pointwise', 'keep_collapsed'}, default 'valid_output'
+            This parameter determines the way a precision reduction is applied on the
+            geometry. There are three modes:
+
+            * ``'valid_output'`` (default): The output is always valid. Collapsed
+              geometry elements (including both polygons and lines) are removed.
+              Duplicate vertices are removed.
+            * ``'pointwise'``: Precision reduction is performed pointwise. Output
+              geometry may be invalid due to collapse or self-intersection. Duplicate
+              vertices are not removed.
+            * ``'keep_collapsed'``: Like the default mode, except that collapsed linear
+              geometry elements are preserved. Collapsed polygonal input elements are
+              removed. Duplicate vertices are removed.
+
+        Examples
+        --------
+
+        >>> from shapely import LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...        Point(0.9, 0.9),
+        ...        Point(0.9, 0.9, 0.9),
+        ...        LineString([(0, 0), (0, 0.1), (0, 1), (1, 1)]),
+        ...        LineString([(0, 0), (0, 0.1), (0.1, 0.1)])
+        ...     ],
+        ... )
+        >>> s
+        0                      POINT (0.9 0.9)
+        1                POINT Z (0.9 0.9 0.9)
+        2    LINESTRING (0 0, 0 0.1, 0 1, 1 1)
+        3     LINESTRING (0 0, 0 0.1, 0.1 0.1)
+        dtype: geometry
+
+        >>> s.set_precision(1)
+        0                   POINT (1 1)
+        1             POINT Z (1 1 0.9)
+        2    LINESTRING (0 0, 0 1, 1 1)
+        3            LINESTRING Z EMPTY
+        dtype: geometry
+
+        >>> s.set_precision(1, mode="pointwise")
+        0                        POINT (1 1)
+        1                  POINT Z (1 1 0.9)
+        2    LINESTRING (0 0, 0 0, 0 1, 1 1)
+        3         LINESTRING (0 0, 0 0, 0 0)
+        dtype: geometry
+
+        >>> s.set_precision(1, mode="keep_collapsed")
+        0                   POINT (1 1)
+        1             POINT Z (1 1 0.9)
+        2    LINESTRING (0 0, 0 1, 1 1)
+        3         LINESTRING (0 0, 0 0)
+        dtype: geometry
+
+        Notes
+        -----
+        Subsequent operations will always be performed in the precision of the geometry
+        with higher precision (smaller ``grid_size``). That same precision will be
+        attached to the operation outputs.
+
+        Input geometries should be geometrically valid; unexpected results may occur if
+        input geometries are not. You can check the validity with
+        :meth:`~GeoSeries.is_valid` and fix invalid geometries with
+        :meth:`~GeoSeries.make_valid` methods.
+
+        """
+        return _delegate_geo_method(
+            "set_precision", self, grid_size=grid_size, mode=mode
+        )
+
     def representative_point(self):
         """Returns a ``GeoSeries`` of (cheaply computed) points that are
         guaranteed to be within each geometry.
@@ -1768,6 +1907,117 @@ GeometryCollection
         """
         return _binary_op("contains_properly", self, other, align)
 
+    def dwithin(self, other, distance, align=True):
+        """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
+        each aligned geometry that is within a set distance from ``other``.
+
+        The operation works on a 1-to-1 row-wise manner:
+
+        .. image:: ../../../_static/binary_op-01.svg
+           :align: center
+
+        Parameters
+        ----------
+        other : GeoSeries or geometric object
+            The GeoSeries (elementwise) or geometric object to test for
+            equality.
+        distance : float, np.array, pd.Series
+            Distance(s) to test if each geometry is within. A scalar distance will be
+            applied to all geometries. An array or Series will be applied elementwise.
+            If np.array or pd.Series are used then it must have same length as the
+            GeoSeries.
+        align : bool (default True)
+            If True, automatically aligns GeoSeries based on their indices.
+            If False, the order of elements is preserved.
+
+        Returns
+        -------
+        Series (bool)
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1), (0, 1)]),
+        ...         LineString([(0, 0), (0, 2)]),
+        ...         LineString([(0, 0), (0, 1)]),
+        ...         Point(0, 1),
+        ...     ],
+        ...     index=range(0, 4),
+        ... )
+        >>> s2 = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(1, 0), (4, 2), (2, 2)]),
+        ...         Polygon([(2, 0), (3, 2), (2, 2)]),
+        ...         LineString([(2, 0), (2, 2)]),
+        ...         Point(1, 1),
+        ...     ],
+        ...     index=range(1, 5),
+        ... )
+
+        >>> s
+        0    POLYGON ((0 0, 1 1, 0 1, 0 0))
+        1             LINESTRING (0 0, 0 2)
+        2             LINESTRING (0 0, 0 1)
+        3                       POINT (0 1)
+        dtype: geometry
+
+        >>> s2
+        1    POLYGON ((1 0, 4 2, 2 2, 1 0))
+        2    POLYGON ((2 0, 3 2, 2 2, 2 0))
+        3             LINESTRING (2 0, 2 2)
+        4                       POINT (1 1)
+        dtype: geometry
+
+        We can check if each geometry of GeoSeries contains a single
+        geometry:
+
+        .. image:: ../../../_static/binary_op-03.svg
+           :align: center
+
+        >>> point = Point(0, 1)
+        >>> s2.dwithin(point, 1.8)
+        1     True
+        2    False
+        3    False
+        4     True
+        dtype: bool
+
+        We can also check two GeoSeries against each other, row by row.
+        The GeoSeries above have different indices. We can either align both GeoSeries
+        based on index values and compare elements with the same index using
+        ``align=True`` or ignore index and compare elements based on their matching
+        order using ``align=False``:
+
+        .. image:: ../../../_static/binary_op-02.svg
+
+        >>> s.dwithin(s2, distance=1, align=True)
+        0    False
+        1     True
+        2    False
+        3    False
+        4    False
+        dtype: bool
+
+        >>> s.dwithin(s2, distance=1, align=False)
+        0     True
+        1    False
+        2    False
+        3     True
+        dtype: bool
+
+        Notes
+        -----
+        This method works in a row-wise manner. It does not check if an element
+        of one GeoSeries is within the set distance of *any* element of the other one.
+
+        See also
+        --------
+        GeoSeries.within
+        """
+        return _binary_op("dwithin", self, other, distance=distance, align=align)
+
     def geom_equals(self, other, align=True):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry equal to `other`.
@@ -4994,11 +5244,13 @@ GeometryCollection
                 )
             sample_function = getattr(pointpats.random, method)
             result = self.geometry.apply(
-                lambda x: points_from_xy(
-                    *sample_function(x, size=size, **kwargs).T
-                ).union_all()
-                if not (x.is_empty or x is None or "Polygon" not in x.geom_type)
-                else MultiPoint(),
+                lambda x: (
+                    points_from_xy(
+                        *sample_function(x, size=size, **kwargs).T
+                    ).union_all()
+                    if not (x.is_empty or x is None or "Polygon" not in x.geom_type)
+                    else MultiPoint()
+                ),
             )
 
         return GeoSeries(result, name="sampled_points", crs=self.crs, index=self.index)
@@ -5050,7 +5302,7 @@ def _get_index_for_parts(orig_idx, outer_idx, ignore_index, index_parts):
             index_arrays.append(inner_index)
 
             index = pd.MultiIndex.from_arrays(
-                index_arrays, names=orig_idx.names + [None]
+                index_arrays, names=list(orig_idx.names) + [None]
             )
 
         else:
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 42bc8e68..63051078 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1850,7 +1850,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         return aggregated
 
     # overrides the pandas native explode method to break up features geometrically
-    def explode(self, column=None, ignore_index=False, index_parts=None, **kwargs):
+    def explode(self, column=None, ignore_index=False, index_parts=False, **kwargs):
         """
         Explode multi-part geometries into multiple single geometries.
 
@@ -1867,7 +1867,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         ignore_index : bool, default False
             If True, the resulting index will be labelled 0, 1, …, n - 1,
             ignoring `index_parts`.
-        index_parts : boolean, default True
+        index_parts : boolean, default False
             If True, the resulting index will be a multi-index (original
             index with an additional level indicating the multiple
             geometries: a new zero-based index for each single part geometry
@@ -1933,18 +1933,6 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         if not isinstance(self[column].dtype, GeometryDtype):
             return super().explode(column, ignore_index=ignore_index, **kwargs)
 
-        if index_parts is None:
-            if not ignore_index:
-                warnings.warn(
-                    "Currently, index_parts defaults to True, but in the future, "
-                    "it will default to False to be consistent with Pandas. "
-                    "Use `index_parts=True` to keep the current behavior and "
-                    "True/False to silence the warning.",
-                    FutureWarning,
-                    stacklevel=2,
-                )
-            index_parts = True
-
         exploded_geom = self.geometry.reset_index(drop=True).explode(index_parts=True)
 
         df = self.drop(self._geometry_column_name, axis=1).take(
@@ -2059,47 +2047,6 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             self, name, con, schema, if_exists, index, index_label, chunksize, dtype
         )
 
-        #
-        # Implement standard operators for GeoSeries
-        #
-
-    def __xor__(self, other):
-        """Implement ^ operator as for builtin set type"""
-        warnings.warn(
-            "'^' operator will be deprecated. Use the 'symmetric_difference' "
-            "method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.geometry.symmetric_difference(other)
-
-    def __or__(self, other):
-        """Implement | operator as for builtin set type"""
-        warnings.warn(
-            "'|' operator will be deprecated. Use the 'union' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.geometry.union(other)
-
-    def __and__(self, other):
-        """Implement & operator as for builtin set type"""
-        warnings.warn(
-            "'&' operator will be deprecated. Use the 'intersection' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.geometry.intersection(other)
-
-    def __sub__(self, other):
-        """Implement - operator as for builtin set type"""
-        warnings.warn(
-            "'-' operator will be deprecated. Use the 'difference' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.geometry.difference(other)
-
     plot = CachedAccessor("plot", geopandas.plotting.GeoplotAccessor)
 
     @doc(_explore)
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index c4634882..2a5a7480 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -208,9 +208,16 @@ class GeoSeries(GeoPandasBase, Series):
                         "Non geometry data passed to GeoSeries constructor, "
                         f"received data of dtype '{s.dtype}'"
                     )
-            # try to convert to GeometryArray, if fails return plain Series
+            # extract object-dtype numpy array from pandas Series; with CoW this
+            # gives a read-only array, so we try to set the flag back to writeable
+            data = s.to_numpy()
             try:
-                data = from_shapely(s.values, crs)
+                data.flags.writeable = True
+            except ValueError:
+                pass
+            # try to convert to GeometryArray
+            try:
+                data = from_shapely(data, crs)
             except TypeError:
                 raise TypeError(
                     "Non geometry data passed to GeoSeries constructor, "
@@ -899,7 +906,7 @@ class GeoSeries(GeoPandasBase, Series):
         """Interactive map based on folium/leaflet.js"""
         return _explore_geoseries(self, *args, **kwargs)
 
-    def explode(self, ignore_index=False, index_parts=None) -> GeoSeries:
+    def explode(self, ignore_index=False, index_parts=False) -> GeoSeries:
         """
         Explode multi-part geometries into multiple single geometries.
 
@@ -912,7 +919,7 @@ class GeoSeries(GeoPandasBase, Series):
         ignore_index : bool, default False
             If True, the resulting index will be labelled 0, 1, …, n - 1,
             ignoring `index_parts`.
-        index_parts : boolean, default True
+        index_parts : boolean, default False
             If True, the resulting index will be a multi-index (original
             index with an additional level indicating the multiple
             geometries: a new zero-based index for each single part geometry
@@ -950,17 +957,6 @@ class GeoSeries(GeoPandasBase, Series):
         """
         from .base import _get_index_for_parts
 
-        if index_parts is None and not ignore_index:
-            warnings.warn(
-                "Currently, index_parts defaults to True, but in the future, "
-                "it will default to False to be consistent with Pandas. "
-                "Use `index_parts=True` to keep the current behavior and True/False "
-                "to silence the warning.",
-                FutureWarning,
-                stacklevel=2,
-            )
-            index_parts = True
-
         geometries, outer_idx = shapely.get_parts(self.values._data, return_index=True)
 
         index = _get_index_for_parts(
@@ -1290,47 +1286,6 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
         """
         return Series(to_wkt(self.array, **kwargs), index=self.index)
 
-    #
-    # Implement standard operators for GeoSeries
-    #
-
-    def __xor__(self, other):
-        """Implement ^ operator as for builtin set type"""
-        warnings.warn(
-            "'^' operator will be deprecated. Use the 'symmetric_difference' "
-            "method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.symmetric_difference(other)
-
-    def __or__(self, other):
-        """Implement | operator as for builtin set type"""
-        warnings.warn(
-            "'|' operator will be deprecated. Use the 'union' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.union(other)
-
-    def __and__(self, other):
-        """Implement & operator as for builtin set type"""
-        warnings.warn(
-            "'&' operator will be deprecated. Use the 'intersection' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.intersection(other)
-
-    def __sub__(self, other):
-        """Implement - operator as for builtin set type"""
-        warnings.warn(
-            "'-' operator will be deprecated. Use the 'difference' method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.difference(other)
-
     def clip(self, mask, keep_geom_type: bool = False) -> GeoSeries:
         """Clip points, lines, or polygon geometries to the mask extent.
 
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index eaf44ce1..36485bc8 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -507,29 +507,6 @@ def _read_file_pyogrio(path_or_bytes, bbox=None, mask=None, rows=None, **kwargs)
     return pyogrio.read_dataframe(path_or_bytes, bbox=bbox, **kwargs)
 
 
-def read_file(*args, **kwargs):
-    warnings.warn(
-        "geopandas.io.file.read_file() is intended for internal "
-        "use only, and will be deprecated. Use geopandas.read_file() instead.",
-        FutureWarning,
-        stacklevel=2,
-    )
-
-    return _read_file(*args, **kwargs)
-
-
-def to_file(*args, **kwargs):
-    warnings.warn(
-        "geopandas.io.file.to_file() is intended for internal "
-        "use only, and will be deprecated. Use GeoDataFrame.to_file() "
-        "or GeoSeries.to_file() instead.",
-        FutureWarning,
-        stacklevel=2,
-    )
-
-    return _to_file(*args, **kwargs)
-
-
 def _detect_driver(path):
     """
     Attempt to auto-detect driver based on the extension
diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py
index b321bb86..7bcb9c55 100644
--- a/geopandas/io/sql.py
+++ b/geopandas/io/sql.py
@@ -183,19 +183,6 @@ def _read_postgis(
         return (_df_to_geodf(df, geom_col=geom_col, crs=crs) for df in df_generator)
 
 
-def read_postgis(*args, **kwargs):
-    import warnings
-
-    warnings.warn(
-        "geopandas.io.sql.read_postgis() is intended for internal "
-        "use only, and will be deprecated. Use geopandas.read_postgis() instead.",
-        FutureWarning,
-        stacklevel=2,
-    )
-
-    return _read_postgis(*args, **kwargs)
-
-
 def _get_geometry_type(gdf):
     """
     Get basic geometry type of a GeoDataFrame. See more info from:
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 34ba0650..b006f61d 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -1,5 +1,3 @@
-import warnings
-
 import numpy as np
 import shapely
 from shapely.geometry.base import BaseGeometry
@@ -249,82 +247,6 @@ class SpatialIndex:
         else:
             return np.asarray(geometry)
 
-    def query_bulk(self, geometry, predicate=None, sort=False, distance=None):
-        """
-        DEPRECATED: use `query` instead.
-
-        Returns all combinations of each input geometry and geometries in
-        the tree where the envelope of each input geometry intersects with
-        the envelope of a tree geometry.
-
-        In the context of a spatial join, input geometries are the “left”
-        geometries that determine the order of the results, and tree geometries
-        are “right” geometries that are joined against the left geometries.
-        This effectively performs an inner join, where only those combinations
-        of geometries that can be joined based on envelope overlap or optional
-        predicate are returned.
-
-        Parameters
-        ----------
-        geometry : {GeoSeries, GeometryArray, numpy.array of Shapely geometries}
-            Accepts GeoPandas geometry iterables (GeoSeries, GeometryArray)
-            or a numpy array of Shapely geometries.
-        predicate : {None, "contains", "contains_properly", "covered_by", "covers", \
-"crosses", "intersects", "overlaps", "touches", "within"}, optional
-            If predicate is provided, the input geometries are tested using
-            the predicate function against each item in the tree whose extent
-            intersects the envelope of the each input geometry:
-            predicate(input_geometry, tree_geometry).  If possible, prepared
-            geometries are used to help speed up the predicate operation.
-        sort : bool, default False
-            If True, results sorted lexicographically using
-            geometry's indexes as the primary key and the sindex's indexes as the
-            secondary key. If False, no additional sorting is applied.
-
-        Returns
-        -------
-        ndarray with shape (2, n)
-            The first subarray contains input geometry integer indexes.
-            The second subarray contains tree geometry integer indexes.
-
-        Examples
-        --------
-        >>> from shapely.geometry import Point, box
-        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-        >>> s
-        0    POINT (0 0)
-        1    POINT (1 1)
-        2    POINT (2 2)
-        3    POINT (3 3)
-        4    POINT (4 4)
-        5    POINT (5 5)
-        6    POINT (6 6)
-        7    POINT (7 7)
-        8    POINT (8 8)
-        9    POINT (9 9)
-        dtype: geometry
-        >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
-        >>> s2
-        0    POLYGON ((4 2, 4 4, 2 4, 2 2, 4 2))
-        1    POLYGON ((6 5, 6 6, 5 6, 5 5, 6 5))
-        dtype: geometry
-
-        >>> s.sindex.query_bulk(s2)
-        array([[0, 0, 0, 1, 1],
-                [2, 3, 4, 5, 6]])
-
-        >>> s.sindex.query_bulk(s2, predicate="contains")
-        array([[0],
-                [3]])
-        """
-        warnings.warn(
-            "The `query_bulk()` method is deprecated and will be removed in "
-            "GeoPandas 1.0. You can use the `query()` method instead.",
-            FutureWarning,
-            stacklevel=2,
-        )
-        return self.query(geometry, predicate=predicate, sort=sort, distance=distance)
-
     def nearest(
         self,
         geometry,
diff --git a/geopandas/tools/clip.py b/geopandas/tools/clip.py
index a47d02f6..728ff0d7 100644
--- a/geopandas/tools/clip.py
+++ b/geopandas/tools/clip.py
@@ -5,6 +5,7 @@ geopandas.clip
 A module to clip vector data using GeoPandas.
 
 """
+
 import warnings
 
 import numpy as np
@@ -60,13 +61,13 @@ def _clip_gdf_with_mask(gdf, mask):
     if isinstance(gdf_sub, GeoDataFrame):
         clipped = gdf_sub.copy()
         if clipping_by_rectangle:
-            clipped.loc[
-                non_point_mask, clipped._geometry_column_name
-            ] = gdf_sub.geometry.values[non_point_mask].clip_by_rect(*mask)
+            clipped.loc[non_point_mask, clipped._geometry_column_name] = (
+                gdf_sub.geometry.values[non_point_mask].clip_by_rect(*mask)
+            )
         else:
-            clipped.loc[
-                non_point_mask, clipped._geometry_column_name
-            ] = gdf_sub.geometry.values[non_point_mask].intersection(mask)
+            clipped.loc[non_point_mask, clipped._geometry_column_name] = (
+                gdf_sub.geometry.values[non_point_mask].intersection(mask)
+            )
     else:
         # GeoSeries
         clipped = gdf_sub.copy()
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 71c9f495..97425da3 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -1,3 +1,4 @@
+from functools import partial
 from typing import Optional
 import warnings
 
@@ -97,26 +98,6 @@ def sjoin(
     Every operation in GeoPandas is planar, i.e. the potential third
     dimension is not taken into account.
     """
-    if "op" in kwargs:
-        op = kwargs.pop("op")
-        deprecation_message = (
-            "The `op` parameter is deprecated and will be removed"
-            " in a future release. Please use the `predicate` parameter"
-            " instead."
-        )
-        if predicate not in ("intersects", op):
-            override_message = (
-                "A non-default value for `predicate` was passed"
-                f' (got `predicate="{predicate}"`'
-                f' in combination with `op="{op}"`).'
-                " The value of `predicate` will be overridden by the value of `op`,"
-                " , which may result in unexpected behavior."
-                f"\n{deprecation_message}"
-            )
-            warnings.warn(override_message, UserWarning, stacklevel=4)
-        else:
-            warnings.warn(deprecation_message, FutureWarning, stacklevel=4)
-        predicate = op
     if kwargs:
         first = next(iter(kwargs.keys()))
         raise TypeError(f"sjoin() got an unexpected keyword argument '{first}'")
@@ -167,18 +148,6 @@ def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
     if not _check_crs(left_df, right_df):
         _crs_mismatch_warn(left_df, right_df, stacklevel=4)
 
-    index_left = "index_{}".format(lsuffix)
-    index_right = "index_{}".format(rsuffix)
-
-    # due to GH 352
-    if any(left_df.columns.isin([index_left, index_right])) or any(
-        right_df.columns.isin([index_left, index_right])
-    ):
-        raise ValueError(
-            "'{0}' and '{1}' cannot be names in the frames being"
-            " joined".format(index_left, index_right)
-        )
-
 
 def _geom_predicate_query(left_df, right_df, predicate, distance):
     """Compute geometric comparisons and get matching indices.
@@ -196,27 +165,21 @@ def _geom_predicate_query(left_df, right_df, predicate, distance):
         DataFrame with matching indices in
         columns named `_key_left` and `_key_right`.
     """
-    with warnings.catch_warnings():
-        # We don't need to show our own warning here
-        # TODO remove this once the deprecation has been enforced
-        warnings.filterwarnings(
-            "ignore", "Generated spatial index is empty", FutureWarning
-        )
 
-        original_predicate = predicate
+    original_predicate = predicate
 
-        if predicate == "within":
-            # within is implemented as the inverse of contains
-            # contains is a faster predicate
-            # see discussion at https://github.com/geopandas/geopandas/pull/1421
-            predicate = "contains"
-            sindex = left_df.sindex
-            input_geoms = right_df.geometry
-        else:
-            # all other predicates are symmetric
-            # keep them the same
-            sindex = right_df.sindex
-            input_geoms = left_df.geometry
+    if predicate == "within":
+        # within is implemented as the inverse of contains
+        # contains is a faster predicate
+        # see discussion at https://github.com/geopandas/geopandas/pull/1421
+        predicate = "contains"
+        sindex = left_df.sindex
+        input_geoms = right_df.geometry
+    else:
+        # all other predicates are symmetric
+        # keep them the same
+        sindex = right_df.sindex
+        input_geoms = left_df.geometry
 
     if sindex:
         l_idx, r_idx = sindex.query(
@@ -237,6 +200,101 @@ def _geom_predicate_query(left_df, right_df, predicate, distance):
     return indices
 
 
+def _reset_index_with_suffix(df, suffix, other):
+    """
+    Equivalent of df.reset_index(), but with adding 'suffix' to auto-generated
+    column names.
+    """
+    index_original = df.index.names
+    df_reset = df.reset_index()
+    column_names = df_reset.columns.to_numpy(copy=True)
+    for i, label in enumerate(index_original):
+        # if the original label was None, add suffix to auto-generated name
+        if label is None:
+            new_label = column_names[i]
+            if "level" in new_label:
+                # reset_index of MultiIndex gives "level_i" names, preserve the "i"
+                lev = new_label.split("_")[1]
+                new_label = f"index_{suffix}{lev}"
+            else:
+                new_label = f"index_{suffix}"
+            # check new label will not be in other dataframe
+            if new_label in df.columns or new_label in other.columns:
+                raise ValueError(
+                    "'{0}' cannot be a column name in the frames being"
+                    " joined".format(new_label)
+                )
+            column_names[i] = new_label
+    return df_reset, pd.Index(column_names)
+
+
+def _process_column_names_with_suffix(
+    left: pd.Index, right: pd.Index, suffixes, left_df, right_df
+):
+    """
+    Add suffixes to overlapping labels (ignoring the geometry column).
+
+    This is based on pandas' merge logic at https://github.com/pandas-dev/pandas/blob/
+    a0779adb183345a8eb4be58b3ad00c223da58768/pandas/core/reshape/merge.py#L2300-L2370
+    """
+    to_rename = left.intersection(right)
+    if len(to_rename) == 0:
+        return left, right
+
+    lsuffix, rsuffix = suffixes
+
+    if not lsuffix and not rsuffix:
+        raise ValueError(f"columns overlap but no suffix specified: {to_rename}")
+
+    def renamer(x, suffix, geometry):
+        if x in to_rename and x != geometry and suffix is not None:
+            return f"{x}_{suffix}"
+        return x
+
+    lrenamer = partial(renamer, suffix=lsuffix, geometry=left_df._geometry_column_name)
+    rrenamer = partial(renamer, suffix=rsuffix, geometry=right_df._geometry_column_name)
+
+    # TODO retain index name?
+    left_renamed = pd.Index([lrenamer(lab) for lab in left])
+    right_renamed = pd.Index([rrenamer(lab) for lab in right])
+
+    dups = []
+    if not left_renamed.is_unique:
+        # Only warn when duplicates are caused because of suffixes, already duplicated
+        # columns in origin should not warn
+        dups = left_renamed[(left_renamed.duplicated()) & (~left.duplicated())].tolist()
+    if not right_renamed.is_unique:
+        dups.extend(
+            right_renamed[(right_renamed.duplicated()) & (~right.duplicated())].tolist()
+        )
+    # TODO turn this into an error (pandas has done so as well)
+    if dups:
+        warnings.warn(
+            f"Passing 'suffixes' which cause duplicate columns {set(dups)} in the "
+            f"result is deprecated and will raise a MergeError in a future version.",
+            FutureWarning,
+            stacklevel=4,
+        )
+
+    return left_renamed, right_renamed
+
+
+def _restore_index(joined, index_names, index_names_original):
+    """
+    Set back the the original index columns, and restoring their name as `None`
+    if they didn't have a name originally.
+    """
+    joined = joined.set_index(list(index_names))
+
+    # restore the fact that the index didn't have a name
+    joined_index_names = list(joined.index.names)
+    for i, label in enumerate(index_names_original):
+        if label is None:
+            joined_index_names[i] = None
+    joined.index.names = joined_index_names
+    return joined
+
+
 def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
     """Join the GeoDataFrames at the DataFrame level.
 
@@ -267,33 +325,29 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
     # index in geopandas may be any arbitrary dtype. so reset both indices now
     # and store references to the original indices, to be reaffixed later.
     # GH 352
-    index_left = "index_{}".format(lsuffix)
     left_df = left_df.copy(deep=True)
-    try:
-        left_index_name = left_df.index.name
-        left_df.index = left_df.index.rename(index_left)
-    except TypeError:
-        index_left = [
-            "index_{}".format(lsuffix + str(pos))
-            for pos, ix in enumerate(left_df.index.names)
-        ]
-        left_index_name = left_df.index.names
-        left_df.index = left_df.index.rename(index_left)
-    left_df = left_df.reset_index()
-
-    index_right = "index_{}".format(rsuffix)
+    left_nlevels = left_df.index.nlevels
+    left_index_original = left_df.index.names
+    left_df, left_column_names = _reset_index_with_suffix(left_df, lsuffix, right_df)
+
     right_df = right_df.copy(deep=True)
-    try:
-        right_index_name = right_df.index.name
-        right_df.index = right_df.index.rename(index_right)
-    except TypeError:
-        index_right = [
-            "index_{}".format(rsuffix + str(pos))
-            for pos, ix in enumerate(right_df.index.names)
-        ]
-        right_index_name = right_df.index.names
-        right_df.index = right_df.index.rename(index_right)
-    right_df = right_df.reset_index()
+    right_nlevels = right_df.index.nlevels
+    right_index_original = right_df.index.names
+    right_df, right_column_names = _reset_index_with_suffix(right_df, rsuffix, left_df)
+
+    # if conflicting names in left and right, add suffix
+    left_column_names, right_column_names = _process_column_names_with_suffix(
+        left_column_names,
+        right_column_names,
+        (lsuffix, rsuffix),
+        left_df,
+        right_df,
+    )
+    left_df.columns = left_column_names
+    right_df.columns = right_column_names
+
+    left_index = left_df.columns[:left_nlevels]
+    right_index = right_df.columns[:right_nlevels]
 
     # perform join on the dataframes
     if how == "inner":
@@ -304,15 +358,11 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
                 right_df.drop(right_df.geometry.name, axis=1),
                 left_on="_key_right",
                 right_index=True,
-                suffixes=("_{}".format(lsuffix), "_{}".format(rsuffix)),
+                suffixes=(None, None),
             )
-            .set_index(index_left)
             .drop(["_key_right"], axis=1)
         )
-        if isinstance(index_left, list):
-            joined.index.names = left_index_name
-        else:
-            joined.index.name = left_index_name
+        joined = _restore_index(joined, left_index, left_index_original)
 
     elif how == "left":
         join_df = join_df.set_index("_key_left")
@@ -323,15 +373,11 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
                 how="left",
                 left_on="_key_right",
                 right_index=True,
-                suffixes=("_{}".format(lsuffix), "_{}".format(rsuffix)),
+                suffixes=(None, None),
             )
-            .set_index(index_left)
             .drop(["_key_right"], axis=1)
         )
-        if isinstance(index_left, list):
-            joined.index.names = left_index_name
-        else:
-            joined.index.name = left_index_name
+        joined = _restore_index(joined, left_index, left_index_original)
 
     else:  # how == 'right':
         joined = (
@@ -343,16 +389,12 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
                 left_index=True,
                 right_on="_key_left",
                 how="right",
-                suffixes=("_{}".format(lsuffix), "_{}".format(rsuffix)),
+                suffixes=(None, None),
             )
-            .set_index(index_right)
             .drop(["_key_left", "_key_right"], axis=1)
             .set_geometry(right_df.geometry.name)
         )
-        if isinstance(index_right, list):
-            joined.index.names = right_index_name
-        else:
-            joined.index.name = right_index_name
+        joined = _restore_index(joined, right_index, right_index_original)
 
     return joined
 
diff --git a/pyproject.toml b/pyproject.toml
index 4a46fe0e..44153b6e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,6 +28,28 @@ dependencies = [
     "shapely >= 2.0.0",
 ]
 
+[project.optional-dependencies]
+all = [
+    "psycopg2-binary>=2.8.0",
+    "SQLAlchemy>=1.3",
+    "geopy",
+    "matplotlib>=3.5.0",
+    "mapclassify",
+    "xyzservices",
+    "folium",
+    "GeoAlchemy2",
+    "pyarrow>=8.0.0",
+]
+dev = [
+    "pytest>=3.1.0",
+    "pytest-cov",
+    "pytest-xdist",
+    "codecov",
+    "black",
+    "pre-commit",
+]
+
+
 [project.readme]
 text = """\
 GeoPandas is a project to add support for geographic data to
diff --git a/readthedocs.yml b/readthedocs.yml
index af979a2f..a92d29fb 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,8 +1,13 @@
 version: 2
 build:
-  os: ubuntu-20.04
+  os: ubuntu-22.04
   tools:
-    python: mambaforge-4.10
+    python: mambaforge-latest
+  jobs:
+    post_checkout:
+      # we need the tags for versioneer to work
+      - git fetch origin --depth 150
+      - git fetch --tags
 python:
   install:
   - method: pip
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 02ab08bc..911f439c 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,7 +6,7 @@ shapely>=2.0.0
 packaging
 
 # geodatabase access
-psycopg2>=2.8.0
+psycopg2-binary>=2.8.0
 SQLAlchemy>=1.3
 
 # geocoding
@@ -15,10 +15,13 @@ geopy
 # plotting
 matplotlib>=3.5.0
 mapclassify
+xyzservices
+folium
 
 # testing
 pytest>=3.1.0
 pytest-cov
+pytest-xdist
 codecov
 
 # styling
@@ -27,3 +30,6 @@ pre-commit
 
 # PostGIS writing
 GeoAlchemy2
+
+# parquet
+pyarrow>=8.0.0


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 3f89af96..bbeaadad 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -19,7 +19,7 @@ jobs:
     steps:
       - uses: actions/checkout@v4
       - uses: actions/setup-python@v5
-      - uses: pre-commit/action@v3.0.0
+      - uses: pre-commit/action@v3.0.1
 
   Test:
     needs: Linting
@@ -86,14 +86,14 @@ jobs:
         env:
           PANDAS_COPY_ON_WRITE: ${{ matrix.pandas_copy_on_write || '0' }}
         run: |
-          pytest -v -r a -n auto --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/
+          pytest -v -r fEs -n auto --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/
 
       - name: Test with Copy-on-Write
         if: ${{ always() && matrix.dev }}
         env:
           PANDAS_COPY_ON_WRITE: '1'
         run: |
-          pytest -v -r a -n auto --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/
+          pytest -v -r fEs -n auto --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/
 
       - name: Test with PostGIS
         if: (contains(matrix.env, '310-pd20-conda-forge.yaml') || contains(matrix.env, '311-latest-conda-forge.yaml'))
@@ -106,11 +106,11 @@ jobs:
         run: |
           conda install postgis -c conda-forge
           sh ci/scripts/setup_postgres.sh
-          pytest -v -r a --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/io/tests/test_sql.py | tee /dev/stderr | if grep SKIPPED >/dev/null;then echo "TESTS SKIPPED, FAILING" && exit 1;fi
+          pytest -v -r fEs --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/io/tests/test_sql.py | tee /dev/stderr | if grep SKIPPED >/dev/null;then echo "TESTS SKIPPED, FAILING" && exit 1;fi
 
       - name: Test docstrings
         if: contains(matrix.env, '310-pd20-conda-forge.yaml') && contains(matrix.os, 'ubuntu')
         run: |
           pytest -v --color=yes --doctest-only geopandas --ignore=geopandas/datasets
 
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
diff --git a/geopandas/io/tests/generate_legacy_storage_files.py b/geopandas/io/tests/generate_legacy_storage_files.py
index 453c5c64..7dd1dc48 100644
--- a/geopandas/io/tests/generate_legacy_storage_files.py
+++ b/geopandas/io/tests/generate_legacy_storage_files.py
@@ -19,6 +19,7 @@ pickles and test versus the current data that is generated
 (with master). These are then compared.
 
 """
+
 import os
 import pickle
 import platform
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index 9ca44ae4..ff543369 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -408,12 +408,6 @@ def test_to_file_empty(tmpdir, engine):
         input_empty_df.to_file(tempfilename, engine=engine)
 
 
-def test_to_file_privacy(tmpdir, df_nybb):
-    tempfilename = os.path.join(str(tmpdir), "test.shp")
-    with pytest.warns(FutureWarning):
-        geopandas.io.file.to_file(df_nybb, tempfilename)
-
-
 def test_to_file_schema(tmpdir, df_nybb, engine):
     """
     Ensure that the file is written according to the schema
@@ -1035,11 +1029,6 @@ def test_read_file_empty_shapefile(tmpdir, engine):
     assert all(empty.columns == ["A", "Z", "geometry"])
 
 
-def test_read_file_privacy(tmpdir, df_nybb, nybb_filename):
-    with pytest.warns(FutureWarning):
-        geopandas.io.file.read_file(nybb_filename)
-
-
 class FileNumber(object):
     def __init__(self, tmpdir, base, ext):
         self.tmpdir = str(tmpdir)
diff --git a/geopandas/io/tests/test_pickle.py b/geopandas/io/tests/test_pickle.py
index 07fea2cc..6d962807 100644
--- a/geopandas/io/tests/test_pickle.py
+++ b/geopandas/io/tests/test_pickle.py
@@ -2,6 +2,7 @@
 See generate_legacy_storage_files.py for the creation of the legacy files.
 
 """
+
 import glob
 import os
 import pathlib
diff --git a/geopandas/io/tests/test_sql.py b/geopandas/io/tests/test_sql.py
index 0d48f5f1..286d9b74 100644
--- a/geopandas/io/tests/test_sql.py
+++ b/geopandas/io/tests/test_sql.py
@@ -4,6 +4,7 @@ The spatial database tests may not work without additional system
 configuration. postGIS tests require a test database to have been setup;
 see geopandas.tests.util for more information.
 """
+
 import os
 import warnings
 
@@ -336,14 +337,6 @@ class TestIO:
         # by user; should not be set to 0, as from get_srid failure
         assert df.crs is None
 
-    def test_read_postgis_privacy(self, connection_postgis, df_nybb):
-        con = connection_postgis
-        create_postgis(con, df_nybb)
-
-        sql = "SELECT * FROM nybb;"
-        with pytest.warns(FutureWarning):
-            geopandas.io.sql.read_postgis(sql, con)
-
     def test_write_postgis_default(self, engine_postgis, df_nybb):
         """Tests that GeoDataFrame can be written to PostGIS with defaults."""
         engine = engine_postgis
diff --git a/geopandas/testing.py b/geopandas/testing.py
index f3cb7db9..62328e45 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -1,6 +1,7 @@
 """
 Testing functionality for geopandas objects.
 """
+
 import warnings
 
 import pandas as pd
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 34ae3f74..0e1f7bbc 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -267,16 +267,6 @@ def test_to_wkt():
     assert res[0] is None
 
 
-def test_data():
-    arr = from_shapely(points_no_missing)
-    with pytest.warns(DeprecationWarning):
-        np_arr = arr.data
-
-    assert isinstance(np_arr, np.ndarray)
-    assert arr.to_numpy() is np_arr
-    assert np.asarray(arr) is np_arr
-
-
 def test_as_array():
     arr = from_shapely(points_no_missing)
     np_arr1 = np.asarray(arr)
@@ -316,9 +306,11 @@ def test_predicates_vector_scalar(attr, args):
         assert result.dtype == bool
 
         expected = [
-            getattr(tri, attr if "geom" not in attr else attr[5:])(other, *args)
-            if tri is not None
-            else na_value
+            (
+                getattr(tri, attr if "geom" not in attr else attr[5:])(other, *args)
+                if tri is not None
+                else na_value
+            )
             for tri in triangles
         ]
 
@@ -555,9 +547,11 @@ def test_binary_distance():
     # vector - vector
     result = P[: len(T)].distance(T[::-1])
     expected = [
-        getattr(p, attr)(t)
-        if not ((t is None or t.is_empty) or (p is None or p.is_empty))
-        else na_value
+        (
+            getattr(p, attr)(t)
+            if not ((t is None or t.is_empty) or (p is None or p.is_empty))
+            else na_value
+        )
         for t, p in zip(triangles[::-1], points)
     ]
     np.testing.assert_allclose(result, expected)
@@ -614,9 +608,11 @@ def test_binary_project(normalized):
 
     result = L.project(P, normalized=normalized)
     expected = [
-        line.project(p, normalized=normalized)
-        if line is not None and p is not None
-        else na_value
+        (
+            line.project(p, normalized=normalized)
+            if line is not None and p is not None
+            else na_value
+        )
         for p, line in zip(points, lines)
     ]
     np.testing.assert_allclose(result, expected)
@@ -628,9 +624,13 @@ def test_binary_project(normalized):
 def test_buffer(resolution, cap_style, join_style):
     na_value = None
     expected = [
-        p.buffer(0.1, resolution=resolution, cap_style=cap_style, join_style=join_style)
-        if p is not None
-        else na_value
+        (
+            p.buffer(
+                0.1, resolution=resolution, cap_style=cap_style, join_style=join_style
+            )
+            if p is not None
+            else na_value
+        )
         for p in points
     ]
     result = P.buffer(
@@ -836,7 +836,7 @@ def test_equality_ops():
 
 def test_dir():
     assert "contains" in dir(P)
-    assert "data" in dir(P)
+    assert "to_numpy" in dir(P)
 
 
 def test_chaining():
diff --git a/geopandas/tests/test_decorator.py b/geopandas/tests/test_decorator.py
index 180ba922..efc3d37e 100644
--- a/geopandas/tests/test_decorator.py
+++ b/geopandas/tests/test_decorator.py
@@ -27,18 +27,15 @@ def cumsum(whatever):
     method="cumavg",
     operation="average",
 )
-def cumavg(whatever):
-    ...
+def cumavg(whatever): ...
 
 
 @doc(cumsum, method="cummax", operation="maximum")
-def cummax(whatever):
-    ...
+def cummax(whatever): ...
 
 
 @doc(cummax, method="cummin", operation="minimum")
-def cummin(whatever):
-    ...
+def cummin(whatever): ...
 
 
 def test_docstring_formatting():
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index 68412248..dcc961aa 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -9,7 +9,7 @@ from geopandas import GeoDataFrame, read_file
 from pandas.testing import assert_frame_equal
 import pytest
 
-from geopandas._compat import PANDAS_GE_15, PANDAS_GE_20
+from geopandas._compat import PANDAS_GE_15, PANDAS_GE_20, PANDAS_GE_30
 from geopandas.testing import assert_geodataframe_equal, geom_almost_equals
 
 
@@ -260,6 +260,7 @@ def test_dissolve_categorical():
 
     # when observed=False we get an additional observation
     # that wasn't in the original data
+    none_val = "GEOMETRYCOLLECTION EMPTY" if PANDAS_GE_30 else None
     expected_gdf_observed_false = geopandas.GeoDataFrame(
         {
             "cat": pd.Categorical(["a", "a", "b", "b"]),
@@ -267,7 +268,7 @@ def test_dissolve_categorical():
             "geometry": geopandas.array.from_wkt(
                 [
                     "MULTIPOINT (0 0, 1 1)",
-                    None,
+                    none_val,
                     "POINT (2 2)",
                     "POINT (3 3)",
                 ]
diff --git a/geopandas/tests/test_explore.py b/geopandas/tests/test_explore.py
index 111c736f..f39c470d 100644
--- a/geopandas/tests/test_explore.py
+++ b/geopandas/tests/test_explore.py
@@ -361,12 +361,12 @@ class TestExplore:
         m = self.world.explore(
             style_kwds={
                 "style_function": lambda x: {
-                    "fillColor": "red"
-                    if x["properties"]["gdp_md_est"] < 10**6
-                    else "green",
-                    "color": "black"
-                    if x["properties"]["gdp_md_est"] < 10**6
-                    else "white",
+                    "fillColor": (
+                        "red" if x["properties"]["gdp_md_est"] < 10**6 else "green"
+                    ),
+                    "color": (
+                        "black" if x["properties"]["gdp_md_est"] < 10**6 else "white"
+                    ),
                 }
             }
         )
diff --git a/geopandas/tests/test_extension_array.py b/geopandas/tests/test_extension_array.py
index 7b91d64f..5861ff14 100644
--- a/geopandas/tests/test_extension_array.py
+++ b/geopandas/tests/test_extension_array.py
@@ -13,6 +13,7 @@ A set of fixtures are defined to provide data for the tests (the fixtures
 expected to be available to pytest by the inherited pandas tests).
 
 """
+
 import operator
 
 import numpy as np
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index d65fe5aa..5dbc2e57 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -213,27 +213,6 @@ class TestGeomMethods:
         fcmp = assert_series_equal
         self._binary_op_test(op, expected, a, b, fcmp, True, False, *args, **kwargs)
 
-    def _test_binary_operator(self, op, expected, a, b):
-        """
-        The operators only have GeoSeries on the left, but can have
-        GeoSeries or GeoDataFrame on the right.
-        If GeoDataFrame is on the left, geometry column is used.
-
-        """
-        if isinstance(expected, GeoPandasBase):
-            fcmp = assert_geoseries_equal
-        else:
-
-            def fcmp(a, b):
-                assert geom_equals(a, b)
-
-        if isinstance(b, GeoPandasBase):
-            right_df = True
-        else:
-            right_df = False
-
-        self._binary_op_test(op, expected, a, b, fcmp, False, right_df)
-
     def _binary_op_test(
         self, op, expected, left, right, fcmp, left_df, right_df, *args, **kwargs
     ):
@@ -520,6 +499,17 @@ class TestGeomMethods:
             expected, self.g0.contains_properly(self.g9, align=False)
         )
 
+    @pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="requires GEOS>=3.10")
+    def test_dwithin(self):
+        expected = [True, True, True, False, True, True, False]
+        assert_array_dtype_equal(expected, self.g0.dwithin(self.p0, 6))
+
+        expected = [False, True, True, True, True, True, False, False]
+        with pytest.warns(UserWarning, match="The indices .+ different"):
+            assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=True))
+        expected = [True, True, True, True, False, False, False]
+        assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=False))
+
     def test_length(self):
         expected = Series(np.array([2 + np.sqrt(2), 4]), index=self.g1.index)
         self._test_unary_real("length", expected, self.g1)
@@ -635,14 +625,10 @@ class TestGeomMethods:
 
     def test_hausdorff_distance(self):
         # closest point is (0, 0) in self.p1
-        expected = Series(
-            np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index
-        )
+        expected = Series(np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index)
         assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
 
-        expected = Series(
-            np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index
-        )
+        expected = Series(np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index)
         assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
@@ -671,9 +657,7 @@ class TestGeomMethods:
     )
     def test_frechet_distance(self):
         # closest point is (0, 0) in self.p1
-        expected = Series(
-            np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index
-        )
+        expected = Series(np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index)
         assert_array_dtype_equal(expected, self.na_none.frechet_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
@@ -694,9 +678,7 @@ class TestGeomMethods:
             expected, self.g0.frechet_distance(self.g9, align=False)
         )
 
-        expected = Series(
-            np.array([np.sqrt(100**2 + (100 - 5) ** 2)]), self.g12.index
-        )
+        expected = Series(np.array([np.sqrt(100**2 + (100 - 5) ** 2)]), self.g12.index)
         assert_array_dtype_equal(
             expected, self.g12.frechet_distance(self.g13, densify=0.25)
         )
@@ -1291,8 +1273,7 @@ class TestGeomMethods:
             index=MultiIndex.from_tuples(index, names=expected_index_name),
             crs=4326,
         )
-        with pytest.warns(FutureWarning, match="Currently, index_parts defaults"):
-            assert_geoseries_equal(expected, s.explode())
+        assert_geoseries_equal(expected, s.explode(index_parts=True))
 
     @pytest.mark.parametrize("index_name", [None, "test"])
     def test_explode_geodataframe(self, index_name):
@@ -1300,8 +1281,7 @@ class TestGeomMethods:
         df = GeoDataFrame({"col": [1, 2], "geometry": s})
         df.index.name = index_name
 
-        with pytest.warns(FutureWarning, match="Currently, index_parts defaults"):
-            test_df = df.explode()
+        test_df = df.explode(index_parts=True)
 
         expected_s = GeoSeries([Point(1, 2), Point(2, 3), Point(5, 5)])
         expected_df = GeoDataFrame({"col": [1, 1, 2], "geometry": expected_s})
@@ -1631,47 +1611,6 @@ class TestGeomMethods:
         assert test_df.geometry.name == test_df._geometry_column_name
         assert "geometry" in test_df.columns
 
-    #
-    # Test '&', '|', '^', and '-'
-    #
-    def test_intersection_operator(self):
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__and__", self.t1, self.g1, self.g2)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__and__", self.t1, self.gdf1, self.g2)
-
-    def test_union_operator(self):
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__or__", self.sq, self.g1, self.g2)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__or__", self.sq, self.gdf1, self.g2)
-
-    def test_union_operator_polygon(self):
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__or__", self.sq, self.g1, self.t2)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__or__", self.sq, self.gdf1, self.t2)
-
-    def test_symmetric_difference_operator(self):
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__xor__", self.sq, self.g3, self.g4)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__xor__", self.sq, self.gdf3, self.g4)
-
-    def test_difference_series2(self):
-        expected = GeoSeries([GeometryCollection(), self.t2])
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__sub__", expected, self.g1, self.g2)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__sub__", expected, self.gdf1, self.g2)
-
-    def test_difference_poly2(self):
-        expected = GeoSeries([self.t1, self.t1])
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__sub__", expected, self.g1, self.t2)
-        with pytest.warns(FutureWarning):
-            self._test_binary_operator("__sub__", expected, self.gdf1, self.t2)
-
     def test_get_coordinates(self):
         expected = DataFrame(
             data=self.expected_2d,
@@ -1855,3 +1794,62 @@ class TestGeomMethods:
             ],
         )
         assert_geoseries_equal(expected, self.g1.force_3d([1, 2]))
+
+    @pytest.mark.skipif(
+        shapely.geos_version < (3, 9, 5), reason="Empty geom bug in GEOS<3.9.5"
+    )
+    def test_set_precision(self):
+        expected = GeoSeries(
+            [
+                Point(-74, 41, 30.3244),
+                Point(-74, 41, 31.2344),
+                Point(-74, 41),
+                self.pt_empty,
+            ],
+            crs=4326,
+        )
+        assert_geoseries_equal(expected, self.landmarks_mixed_empty.set_precision(1))
+
+        s = GeoSeries(
+            [
+                LineString([(0, 0), (0, 0.1), (0, 1), (1, 1)]),
+                LineString([(0, 0), (0, 0.1), (0.1, 0.1)]),
+            ],
+        )
+        expected = GeoSeries(
+            [
+                LineString([(0, 0), (0, 1), (1, 1)]),
+                LineString(),
+            ],
+        )
+        assert_geoseries_equal(expected, s.set_precision(1))
+
+        expected = GeoSeries(
+            [
+                LineString([(0, 0), (0, 0), (0, 1), (1, 1)]),
+                LineString([(0, 0), (0, 0), (0, 0)]),
+            ]
+        )
+        assert_series_equal(
+            expected.to_wkt(), s.set_precision(1, mode="pointwise").to_wkt()
+        )
+
+        expected = GeoSeries(
+            [
+                LineString([(0, 0), (0, 1), (1, 1)]),
+                LineString([(0, 0), (0, 0)]),
+            ]
+        )
+        assert_series_equal(
+            expected.to_wkt(), s.set_precision(1, mode="keep_collapsed").to_wkt()
+        )
+
+    def test_get_precision(self):
+        expected = Series([0.0, 0.0, 0.0, 0.0], index=self.landmarks_mixed_empty.index)
+        assert_series_equal(expected, self.landmarks_mixed_empty.get_precision())
+        with_precision = self.landmarks_mixed_empty.set_precision(1)
+        expected = Series([1.0, 1.0, 1.0, 1.0], index=with_precision.index)
+        assert_series_equal(expected, with_precision.get_precision())
+        mixed = concat([self.landmarks_mixed_empty, with_precision])
+        expected = Series([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], index=mixed.index)
+        assert_series_equal(expected, mixed.get_precision())
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 4ebf2463..acdd287d 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -24,6 +24,7 @@ from shapely.geometry.base import BaseGeometry
 
 from geopandas import GeoSeries, GeoDataFrame, read_file, clip
 from geopandas.array import GeometryArray, GeometryDtype
+import geopandas._compat as compat
 from geopandas.testing import assert_geoseries_equal, geom_almost_equals
 
 from geopandas.tests.util import geom_equals
@@ -563,6 +564,25 @@ class TestConstructor:
         assert s.name == g.name
         assert s.index is g.index
 
+    def test_copy(self):
+        # default is to copy with CoW / pandas 3+
+        arr = np.array([Point(x, x) for x in range(3)], dtype=object)
+        result = GeoSeries(arr)
+        # modifying result doesn't change original array
+        result.loc[0] = Point(10, 10)
+        if compat.PANDAS_GE_30 or getattr(pd.options.mode, "copy_on_write", False):
+            assert arr[0] == Point(0, 0)
+        else:
+            assert arr[0] == Point(10, 10)
+
+        # avoid copy with copy=False
+        arr = np.array([Point(x, x) for x in range(3)], dtype=object)
+        result = GeoSeries(arr, copy=False)
+        assert result.array._data.flags.writeable
+        # now modifying result also updates original array
+        result.loc[0] = Point(10, 10)
+        assert arr[0] == Point(10, 10)
+
     # GH 1216
     @pytest.mark.parametrize("name", [None, "geometry", "Points"])
     @pytest.mark.parametrize("crs", [None, "epsg:4326"])
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 7014b761..85c0c73c 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -749,6 +749,7 @@ def test_apply_loc_len1(df):
     np.testing.assert_allclose(result, expected)
 
 
+@pytest.mark.skipif(compat.PANDAS_GE_30, reason="convert_dtype is removed in pandas 3")
 def test_apply_convert_dtypes_keyword(s):
     # ensure the convert_dtypes keyword is accepted
     if not compat.PANDAS_GE_21:
@@ -874,3 +875,17 @@ def test_preserve_flags(df):
 
     with pytest.raises(ValueError):
         pd.concat([df, df])
+
+
+def test_ufunc():
+    # this is calling a shapely ufunc, but we currently rely on pandas' implementation
+    # of `__array_ufunc__` to wrap the result back into a GeoSeries
+    ser = GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)])
+    result = shapely.buffer(ser, 2)
+    assert isinstance(result, GeoSeries)
+
+    # ensure the result is still writeable
+    # (https://github.com/geopandas/geopandas/issues/3178)
+    assert result.array._data.flags.writeable
+    result.loc[0] = Point(10, 10)
+    assert result.iloc[0] == Point(10, 10)
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index 51fb5dc0..43be87e2 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -647,7 +647,7 @@ class TestShapelyInterface:
         ),
     )
     def test_query_bulk(self, predicate, test_geom, expected):
-        """Tests the `query_bulk` method with valid
+        """Tests the `query` method with valid
         inputs and valid predicates.
         """
         res = self.df.sindex.query(
@@ -667,12 +667,12 @@ class TestShapelyInterface:
         ],
     )
     def test_query_bulk_empty_geometry(self, test_geoms, expected_value):
-        """Tests the `query_bulk` method with an empty geometry."""
+        """Tests the `query` method with an empty geometries."""
         res = self.df.sindex.query(test_geoms)
         assert_array_equal(res, expected_value)
 
     def test_query_bulk_empty_input_array(self):
-        """Tests the `query_bulk` method with an empty input array."""
+        """Tests the `query` method with an empty input array."""
         test_array = np.array([], dtype=object)
         expected_value = [[], []]
         res = self.df.sindex.query(test_array)
@@ -680,14 +680,14 @@ class TestShapelyInterface:
 
     def test_query_bulk_invalid_input_geometry(self):
         """
-        Tests the `query_bulk` method with invalid input for the `geometry` parameter.
+        Tests the `query` method with invalid input for the `geometry` parameter.
         """
         test_array = "notanarray"
         with pytest.raises(TypeError):
             self.df.sindex.query(test_array)
 
     def test_query_bulk_invalid_predicate(self):
-        """Tests the `query_bulk` method with invalid predicates."""
+        """Tests the `query` method with invalid predicates."""
         test_geom_bounds = (-1, -1, -0.5, -0.5)
         test_predicate = "test"
 
@@ -703,7 +703,7 @@ class TestShapelyInterface:
         ),
     )
     def test_query_bulk_input_type(self, predicate, test_geom, expected):
-        """Tests that query_bulk can accept a GeoSeries, GeometryArray or
+        """Tests that query can accept a GeoSeries, GeometryArray or
         numpy array.
         """
         # pass through GeoSeries to test input type
@@ -738,7 +738,7 @@ class TestShapelyInterface:
         ),
     )
     def test_query_bulk_sorting(self, sort, expected):
-        """Check that results from `query_bulk` don't depend
+        """Check that results from `query` don't depend
         on the order of geometries.
         """
         # these geometries come from a reported issue:
diff --git a/geopandas/tools/tests/test_clip.py b/geopandas/tools/tests/test_clip.py
index 8cc51c9d..cac8e09b 100644
--- a/geopandas/tools/tests/test_clip.py
+++ b/geopandas/tools/tests/test_clip.py
@@ -1,6 +1,5 @@
 """Tests for the clip module."""
 
-
 import numpy as np
 
 import shapely
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index 1c3f79a0..beeecb02 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -114,31 +114,6 @@ class TestSpatialJoin:
         with pytest.warns(UserWarning, match="CRS mismatch between the CRS"):
             sjoin(df1, df2)
 
-    @pytest.mark.parametrize("dfs", ["default-index"], indirect=True)
-    @pytest.mark.parametrize("op", ["intersects", "contains", "within"])
-    def test_deprecated_op_param(self, dfs, op):
-        _, df1, df2, _ = dfs
-        with pytest.warns(FutureWarning, match="`op` parameter is deprecated"):
-            sjoin(df1, df2, op=op)
-
-    @pytest.mark.parametrize("dfs", ["default-index"], indirect=True)
-    @pytest.mark.parametrize("op", ["intersects", "contains", "within"])
-    @pytest.mark.parametrize("predicate", ["contains", "within"])
-    def test_deprecated_op_param_nondefault_predicate(self, dfs, op, predicate):
-        _, df1, df2, _ = dfs
-        match = "use the `predicate` parameter instead"
-        if op != predicate:
-            warntype = UserWarning
-            match = (
-                "`predicate` will be overridden by the value of `op`"  # noqa: ISC003
-                + r"(.|\s)*"
-                + match
-            )
-        else:
-            warntype = FutureWarning
-        with pytest.warns(warntype, match=match):
-            sjoin(df1, df2, predicate=predicate, op=op)
-
     @pytest.mark.parametrize("dfs", ["default-index"], indirect=True)
     def test_unknown_kwargs(self, dfs):
         _, df1, df2, _ = dfs
@@ -148,7 +123,6 @@ class TestSpatialJoin:
         ):
             sjoin(df1, df2, extra_param="test")
 
-    @pytest.mark.filterwarnings("ignore:The `op` parameter:FutureWarning")
     @pytest.mark.parametrize(
         "dfs",
         [
@@ -161,11 +135,10 @@ class TestSpatialJoin:
         indirect=True,
     )
     @pytest.mark.parametrize("predicate", ["intersects", "contains", "within"])
-    @pytest.mark.parametrize("predicate_kw", ["predicate", "op"])
-    def test_inner(self, predicate, predicate_kw, dfs):
+    def test_inner(self, predicate, dfs):
         index, df1, df2, expected = dfs
 
-        res = sjoin(df1, df2, how="inner", **{predicate_kw: predicate})
+        res = sjoin(df1, df2, how="inner", predicate=predicate)
 
         exp = expected[predicate].dropna().copy()
         exp = exp.drop("geometry_y", axis=1).rename(columns={"geometry_x": "geometry"})
@@ -176,7 +149,7 @@ class TestSpatialJoin:
             ].astype("int64")
         if index == "named-index":
             exp[["df1_ix", "df2_ix"]] = exp[["df1_ix", "df2_ix"]].astype("int64")
-            exp = exp.set_index("df1_ix").rename(columns={"df2_ix": "index_right"})
+            exp = exp.set_index("df1_ix")
         if index in ["default-index", "string-index"]:
             exp = exp.set_index("index_left")
             exp.index.name = None
@@ -186,10 +159,7 @@ class TestSpatialJoin:
             )
             exp.index.names = df1.index.names
         if index == "named-multi-index":
-            exp = exp.set_index(["df1_ix1", "df1_ix2"]).rename(
-                columns={"df2_ix1": "index_right0", "df2_ix2": "index_right1"}
-            )
-            exp.index.names = df1.index.names
+            exp = exp.set_index(["df1_ix1", "df1_ix2"])
 
         assert_frame_equal(res, exp)
 
@@ -226,7 +196,7 @@ class TestSpatialJoin:
             res["index_right"] = res["index_right"].astype(float)
         elif index == "named-index":
             exp[["df1_ix"]] = exp[["df1_ix"]].astype("int64")
-            exp = exp.set_index("df1_ix").rename(columns={"df2_ix": "index_right"})
+            exp = exp.set_index("df1_ix")
         if index in ["default-index", "string-index"]:
             exp = exp.set_index("index_left")
             exp.index.name = None
@@ -236,10 +206,7 @@ class TestSpatialJoin:
             )
             exp.index.names = df1.index.names
         if index == "named-multi-index":
-            exp = exp.set_index(["df1_ix1", "df1_ix2"]).rename(
-                columns={"df2_ix1": "index_right0", "df2_ix2": "index_right1"}
-            )
-            exp.index.names = df1.index.names
+            exp = exp.set_index(["df1_ix1", "df1_ix2"])
 
         assert_frame_equal(res, exp)
 
@@ -342,7 +309,7 @@ class TestSpatialJoin:
             res["index_left"] = res["index_left"].astype(float)
         elif index == "named-index":
             exp[["df2_ix"]] = exp[["df2_ix"]].astype("int64")
-            exp = exp.set_index("df2_ix").rename(columns={"df1_ix": "index_left"})
+            exp = exp.set_index("df2_ix")
         if index in ["default-index", "string-index"]:
             exp = exp.set_index("index_right")
             exp = exp.reindex(columns=res.columns)
@@ -353,10 +320,8 @@ class TestSpatialJoin:
             )
             exp.index.names = df2.index.names
         if index == "named-multi-index":
-            exp = exp.set_index(["df2_ix1", "df2_ix2"]).rename(
-                columns={"df1_ix1": "index_left0", "df1_ix2": "index_left1"}
-            )
-            exp.index.names = df2.index.names
+            exp = exp.set_index(["df2_ix1", "df2_ix2"])
+
         if predicate == "within":
             exp = exp.sort_index()
 
@@ -419,6 +384,232 @@ class TestSpatialJoin:
         assert_frame_equal(expected_gdf.sort_index(), joined.sort_index())
 
 
+class TestIndexNames:
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_preserve_index_names(self, how):
+        # preserve names of both left and right index
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame({"geometry": geoms}, index=pd.Index([1, 2], name="myidx1"))
+        df2 = GeoDataFrame(
+            {"geometry": geoms}, index=pd.Index(["a", "b"], name="myidx2")
+        )
+        result = sjoin(df1, df2, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {"myidx1": [1, 2], "geometry": geoms, "myidx2": ["a", "b"]}
+            ).set_index("myidx1")
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {"myidx2": ["a", "b"], "myidx1": [1, 2], "geometry": geoms},
+            ).set_index("myidx2")
+        assert_geodataframe_equal(result, expected)
+
+        # but also add suffixes if both left and right have the same index
+        df1.index.name = "myidx"
+        df2.index.name = "myidx"
+        result = sjoin(df1, df2, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {"myidx_left": [1, 2], "geometry": geoms, "myidx_right": ["a", "b"]}
+            ).set_index("myidx_left")
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {"myidx_right": ["a", "b"], "myidx_left": [1, 2], "geometry": geoms},
+            ).set_index("myidx_right")
+        assert_geodataframe_equal(result, expected)
+
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_preserve_index_names_multiindex(self, how):
+        # preserve names of both left and right index
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame(
+            {"geometry": geoms},
+            index=pd.MultiIndex.from_tuples(
+                [("a", 1), ("b", 2)], names=["myidx1", "level2"]
+            ),
+        )
+        df2 = GeoDataFrame(
+            {"geometry": geoms},
+            index=pd.MultiIndex.from_tuples(
+                [("c", 3), ("d", 4)], names=["myidx2", None]
+            ),
+        )
+        result = sjoin(df1, df2, how=how)
+        expected_base = GeoDataFrame(
+            {
+                "myidx1": ["a", "b"],
+                "level2": [1, 2],
+                "geometry": geoms,
+                "myidx2": ["c", "d"],
+                "index_right1": [3, 4],
+            }
+        )
+        if how in ("inner", "left"):
+            expected = expected_base.set_index(["myidx1", "level2"])
+        else:
+            # right join
+            expected = expected_base.set_index(["myidx2", "index_right1"])
+            # if it was originally None, that is preserved
+            expected.index.names = ["myidx2", None]
+        assert_geodataframe_equal(result, expected)
+
+        # but also add suffixes if both left and right have the same index
+        df1.index.names = ["myidx", "level2"]
+        df2.index.names = ["myidx", None]
+        result = sjoin(df1, df2, how=how)
+        expected_base = GeoDataFrame(
+            {
+                "myidx_left": ["a", "b"],
+                "level2": [1, 2],
+                "geometry": geoms,
+                "myidx_right": ["c", "d"],
+                "index_right1": [3, 4],
+            }
+        )
+        if how in ("inner", "left"):
+            expected = expected_base.set_index(["myidx_left", "level2"])
+        else:
+            # right join
+            expected = expected_base.set_index(["myidx_right", "index_right1"])
+            # if it was originally None, that is preserved
+            expected.index.names = ["myidx_right", None]
+        assert_geodataframe_equal(result, expected)
+
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_duplicate_column_index_name(self, how):
+        # case where a left column and the right index have the same name or the
+        # other way around -> correctly add suffix or preserve index name
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame({"myidx": [1, 2], "geometry": geoms})
+        df2 = GeoDataFrame(
+            {"geometry": geoms}, index=pd.Index(["a", "b"], name="myidx")
+        )
+        result = sjoin(df1, df2, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {"myidx_left": [1, 2], "geometry": geoms, "myidx_right": ["a", "b"]}
+            )
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {"index_left": [0, 1], "myidx_left": [1, 2], "geometry": geoms},
+                index=pd.Index(["a", "b"], name="myidx_right"),
+            )
+        assert_geodataframe_equal(result, expected)
+
+        result = sjoin(df2, df1, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {"geometry": geoms, "index_right": [0, 1], "myidx_right": [1, 2]},
+                index=pd.Index(["a", "b"], name="myidx_left"),
+            )
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {"myidx_left": ["a", "b"], "myidx_right": [1, 2], "geometry": geoms},
+            )
+        assert_geodataframe_equal(result, expected)
+
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_duplicate_column_index_name_multiindex(self, how):
+        # case where a left column and the right index have the same name or the
+        # other way around -> correctly add suffix or preserve index name
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame({"myidx": [1, 2], "geometry": geoms})
+        df2 = GeoDataFrame(
+            {"geometry": geoms},
+            index=pd.MultiIndex.from_tuples(
+                [("a", 1), ("b", 2)], names=["myidx", "level2"]
+            ),
+        )
+        result = sjoin(df1, df2, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {
+                    "myidx_left": [1, 2],
+                    "geometry": geoms,
+                    "myidx_right": ["a", "b"],
+                    "level2": [1, 2],
+                }
+            )
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {"index_left": [0, 1], "myidx_left": [1, 2], "geometry": geoms},
+                index=pd.MultiIndex.from_tuples(
+                    [("a", 1), ("b", 2)], names=["myidx_right", "level2"]
+                ),
+            )
+        assert_geodataframe_equal(result, expected)
+
+        result = sjoin(df2, df1, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {"geometry": geoms, "index_right": [0, 1], "myidx_right": [1, 2]},
+                index=pd.MultiIndex.from_tuples(
+                    [("a", 1), ("b", 2)], names=["myidx_left", "level2"]
+                ),
+            )
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {
+                    "myidx_left": ["a", "b"],
+                    "level2": [1, 2],
+                    "myidx_right": [1, 2],
+                    "geometry": geoms,
+                },
+            )
+        assert_geodataframe_equal(result, expected)
+
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_conflicting_column_index_name(self, how):
+        # test case where the auto-generated index name conflicts
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame({"index_right": [1, 2], "geometry": geoms})
+        df2 = GeoDataFrame({"geometry": geoms})
+        with pytest.raises(ValueError, match="'index_right' cannot be a column name"):
+            sjoin(df1, df2, how=how)
+
+    @pytest.mark.xfail
+    @pytest.mark.parametrize("how", ["inner", "left", "right"])
+    def test_conflicting_column_with_suffix(self, how):
+        # test case where the auto-generated index name conflicts
+        geoms = [Point(1, 1), Point(2, 2)]
+        df1 = GeoDataFrame(
+            {"column": [1, 2], "column_right": ["a", "b"], "geometry": geoms}
+        )
+        df2 = GeoDataFrame({"column": [0.1, 0.2], "geometry": geoms})
+
+        result = sjoin(df1, df2, how=how)
+        if how in ("inner", "left"):
+            expected = GeoDataFrame(
+                {1: [1, 2], 2: ["a", "b"], 3: geoms, 4: [0, 1], 5: [0.1, 0.2]}
+            )
+            expected.columns = [
+                "column_left",
+                "column_right",
+                "geometry",
+                "index_right",
+                "column_right",
+            ]
+        else:
+            # right join
+            expected = GeoDataFrame(
+                {1: [0, 1], 2: [1, 2], 3: ["a", "b"], 4: [0.1, 0.2], 5: geoms}
+            )
+            expected.columns = [
+                "index_left",
+                "column_left",
+                "column_right",
+                "column_right",
+                "geometry",
+            ]
+        assert_geodataframe_equal(result, expected)
+
+
 @pytest.mark.usefixtures("_setup_class_nybb_filename")
 class TestSpatialJoinNYBB:
     def setup_method(self):