Repository: geopandas/geopandas

Issue:
ENH: Improves Examples in Gallery to Better Explain Plotting Coordinates With Background Maps
#### Is your feature request related to a problem?

The example "[Creating a GeoDataFrame from a DataFrame with coordinates](https://geopandas.org/en/stable/gallery/create_geopandas_from_pandas.html)" in the examples gallery creates a `GeoDataFrame` object with no CRS, ie. is `None`. A world geometry is then loaded and the two are plotted together.

This example functions correctly because the loaded world geometry has a CRS of "EPSG:4326" which seems compatible with a None CRS when plotting. However the data frame is not compatible with plotting with Contextily or other libraries.

#### Describe the solution you'd like

The created data frame should have a CRS explicitly set and the example should discuss what this CRS means. With this in place the data frame in this example can be plotted using code from the "[Adding a background map to plot](https://geopandas.org/en/stable/gallery/plotting_basemap_background.html)" example. As an example for less experienced users it's important to do things like this explicitly so that the purpose of CRS is known, this is especially important since the errors one sees from failed plots are not terribly helpful in understanding the root problem.

In general as an inexperienced user of Geopandas and GIS data I found it difficult to go from a set of Lat/Lon degree coordinates to points plotted with a background map. Adding a CRS value to this example would help with that in conjunction with the example on plotting with a background map which explains doing CRS conversion when necessary. A single unified example of a set of coordinates being plotted on a background map, with some further elaboration on other plotting possibilities, would help a lot in this regard.

#### API breaking implications

This is purely documentation. It appears that plotting with a None CRS does work, perhaps this isn't desired behaviour and adding a warning or error would help others discover that this was an issue.

#### Describe alternatives you've considered

I had to figure out everything mentioned here by trial and error, specifically looking at what the CRS of `world` is in the below example and using that. This might still be incorrect but appears to produce a correct map.

#### Additional context

Below is the code from the example with the proposed addition to the `GeoDataFrame` constructor and use of the Contextily library (which I wouldn't put into this example but I leave it here for demonstration):

```python
import pandas as pd
import geopandas
import contextily as cx

df = pd.DataFrame(
    {'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'],
     'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'],
     'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48],
     'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]})

gdf = geopandas.GeoDataFrame(
    df, 
    geometry=geopandas.points_from_xy(df.Longitude, df.Latitude), 
    crs="EPSG:4326"  # <-- additional argument, without it gdf.crs is None
)

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

# this will plot South America correctly
ax = world[world.continent == 'South America'].plot(color='white', edgecolor='black')

# if the CRS of gdf is None this still works
gdf.plot(ax=ax)  

cx.add_basemap(ax, crs=gdf.crs)  # need to set crs argument or move gdf and world to EPSG:3857  
```



Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f755ca98..f6ac18b0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,11 +13,16 @@ New features and improvements:
   the distance along the Hilbert curve (#2070).
 - Added ``minimum_bounding_circle()`` method from shapely to GeoSeries/GeoDataframe (#2621).
 - Support specifying ``min_zoom`` and ``max_zoom`` inside the ``map_kwds`` argument for ``.explore()`` (#2599).
+- Added `minimum_bounding_radius()` as GeoSeries method (#2827).
+- Added support for append (``mode="a"`` or ``append=True``) in ``to_file()``
+  using ``engine="pyogrio"`` (#2788).
 
 Deprecations and compatibility notes:
 
 - Added warning that ``unary_union`` will return ``'GEOMETRYCOLLECTION EMPTY'`` instead
   of None for all-None GeoSeries. (#2618)
+- The ``query_bulk()`` method of the spatial index `.sindex` property is deprecated
+  in favor of ``query()`` (#2823).
 
 Bug fixes:
 
@@ -30,6 +35,8 @@ Bug fixes:
 - Fix `to_parquet`/`to_feather` to use correct WKB flavor for 3D geometries (#2654)
 - Fix `read_file` to avoid reading all file bytes prior to calling Fiona or
   Pyogrio if provided a URL as input (#2796)
+- Fix `copy()` downcasting GeoDataFrames without an active geometry column to a 
+  DataFrame (#2775)
 
 Notes on (optional) dependencies:
 
diff --git a/benchmarks/sindex.py b/benchmarks/sindex.py
index 0767bb1f..bbd15a62 100644
--- a/benchmarks/sindex.py
+++ b/benchmarks/sindex.py
@@ -90,7 +90,7 @@ class BenchQuery:
         self.data = generate_test_df()
 
     def time_query_bulk(self, predicate, input_geom_type, tree_geom_type):
-        self.data[tree_geom_type].sindex.query_bulk(
+        self.data[tree_geom_type].sindex.query(
             self.data[input_geom_type].geometry.values,
             predicate=predicate,
         )
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 21c81598..8c52d8d1 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -175,7 +175,7 @@ html_theme_options = {
             "url": "https://twitter.com/geopandas",
             "icon": "fab fa-twitter-square fa-xl",
         },
-   ]
+    ]
 }
 
 # Add any paths that contain custom themes here, relative to this directory.
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 465c040e..253aa810 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -26,6 +26,7 @@ General methods and attributes
    GeoSeries.representative_point
    GeoSeries.exterior
    GeoSeries.interiors
+   GeoSeries.minimum_bounding_radius
    GeoSeries.x
    GeoSeries.y
    GeoSeries.z
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
index bb8c0d93..5536d6d9 100644
--- a/doc/source/docs/reference/sindex.rst
+++ b/doc/source/docs/reference/sindex.rst
@@ -32,7 +32,6 @@ methods:
     is_empty
     nearest
     query
-    query_bulk
     size
     valid_query_predicates
 
diff --git a/doc/source/gallery/create_geopandas_from_pandas.ipynb b/doc/source/gallery/create_geopandas_from_pandas.ipynb
index 4907e276..a8f9f86c 100644
--- a/doc/source/gallery/create_geopandas_from_pandas.ipynb
+++ b/doc/source/gallery/create_geopandas_from_pandas.ipynb
@@ -62,8 +62,10 @@
     "``points_from_xy()`` to transform **Longitude** and **Latitude** into a list\n",
     "of ``shapely.Point`` objects and set it as a ``geometry`` while creating the\n",
     "``GeoDataFrame``. (note that ``points_from_xy()`` is an enhanced wrapper for\n",
-    "``[Point(x, y) for x, y in zip(df.Longitude, df.Latitude)]``)\n",
-    "\n"
+    "``[Point(x, y) for x, y in zip(df.Longitude, df.Latitude)]``). The ``crs``\n",
+    "value is also set to explicitly state the geometry data defines latitude/\n",
+    "longitude world geodetic degree values. This is important for the correct \n",
+    "interpretation of the data, such as when plotting with data in other formats."
    ]
   },
   {
@@ -73,7 +75,7 @@
    "outputs": [],
    "source": [
     "gdf = geopandas.GeoDataFrame(\n",
-    "    df, geometry=geopandas.points_from_xy(df.Longitude, df.Latitude)\n",
+    "    df, geometry=geopandas.points_from_xy(df.Longitude, df.Latitude), crs=\"EPSG:4326\"\n",
     ")"
    ]
   },
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 19e51fa4..97534d52 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -688,6 +688,18 @@ def minimum_bounding_circle(data):
         )
 
 
+def minimum_bounding_radius(data):
+    if compat.USE_SHAPELY_20:
+        return shapely.minimum_bounding_radius(data)
+    elif compat.USE_PYGEOS:
+        return pygeos.minimum_bounding_radius(data)
+    else:
+        raise NotImplementedError(
+            f"shapely >= 2.0 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+
+
 #
 # Binary predicates
 #
diff --git a/geopandas/array.py b/geopandas/array.py
index a5a3b72d..1474ddf6 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -538,6 +538,9 @@ class GeometryArray(ExtensionArray):
             vectorized.minimum_bounding_circle(self._data), crs=self.crs
         )
 
+    def minimum_bounding_radius(self):
+        return vectorized.minimum_bounding_radius(self._data)
+
     def normalize(self):
         return GeometryArray(vectorized.normalize(self._data), crs=self.crs)
 
diff --git a/geopandas/base.py b/geopandas/base.py
index 6955219e..52fb17fa 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -726,6 +726,38 @@ GeometryCollection
         """
         return _delegate_geo_method("minimum_bounding_circle", self)
 
+    def minimum_bounding_radius(self):
+        """Returns a `Series` of the radii of the minimum bounding circles
+        that enclose each geometry.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point, LineString, Polygon
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1), (0, 1), (0, 0)]),
+        ...         LineString([(0, 0), (1, 1), (1, 0)]),
+        ...         Point(0,0),
+        ...     ]
+        ... )
+        >>> s
+        0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
+        1    LINESTRING (0.00000 0.00000, 1.00000 1.00000, ...
+        2                              POINT (0.00000 0.00000)
+        dtype: geometry
+        >>> s.minimum_bounding_radius()
+        0    0.707107
+        1    0.707107
+        2    0.000000
+        dtype: float64
+
+        See also
+        --------
+        GeoSeries.minumum_bounding_circle : minimum bounding circle (geometry)
+
+        """
+        return Series(self.geometry.values.minimum_bounding_radius(), index=self.index)
+
     def normalize(self):
         """Returns a ``GeoSeries`` of normalized
         geometries to normal form (or canonical form).
@@ -2810,13 +2842,13 @@ GeometryCollection
         1    POLYGON ((5.00000 4.00000, 5.00000 5.00000, 4....
         dtype: geometry
 
-        >>> s.sindex.query_bulk(s2)
+        >>> s.sindex.query(s2)
         array([[0, 0, 0, 1],
                [1, 2, 3, 4]])
 
         Query the spatial index with an array of geometries based on the predicate:
 
-        >>> s.sindex.query_bulk(s2, predicate="contains")
+        >>> s.sindex.query(s2, predicate="contains")
         array([[0],
                [2]])
         """
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 9b816974..887506d3 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1459,6 +1459,13 @@ individually so that features may have different properties
     #
     # Implement pandas methods
     #
+    @doc(pd.DataFrame)
+    def copy(self, deep=True):
+        copied = super().copy(deep=deep)
+        if type(copied) is pd.DataFrame:
+            copied.__class__ = GeoDataFrame
+            copied._geometry_column_name = self._geometry_column_name
+        return copied
 
     def merge(self, *args, **kwargs):
         r"""Merge two ``GeoDataFrame`` objects with a database-style join.
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index 5ab40a22..ef7f9645 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -477,10 +477,15 @@ def _to_file(
         .. versionadded:: 0.7
             Previously the index was not written.
     mode : string, default 'w'
-        The write mode, 'w' to overwrite the existing file and 'a' to append.
-        Not all drivers support appending. The drivers that support appending
-        are listed in fiona.supported_drivers or
-        https://github.com/Toblerity/Fiona/blob/master/fiona/drvsupport.py
+        The write mode, 'w' to overwrite the existing file and 'a' to append;
+        when using the pyogrio engine, you can also pass ``append=True``.
+        Not all drivers support appending. For the fiona engine, the drivers
+        that support appending are listed in fiona.supported_drivers or
+        https://github.com/Toblerity/Fiona/blob/master/fiona/drvsupport.py.
+        For the pyogrio engine, you should be able to use any driver that
+        is available in your installation of GDAL that supports append
+        capability; see the specific driver entry at
+        https://gdal.org/drivers/vector/index.html for more information.
     crs : pyproj.CRS, default None
         If specified, the CRS is passed to Fiona to
         better control how the file is written. If None, GeoPandas
@@ -527,6 +532,9 @@ def _to_file(
             stacklevel=3,
         )
 
+    if mode not in ("w", "a"):
+        raise ValueError("'mode' should be one of 'w' or 'a', got '{mode}' instead")
+
     if engine == "fiona":
         _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs)
     elif engine == "pyogrio":
@@ -568,10 +576,8 @@ def _to_file_pyogrio(df, filename, driver, schema, crs, mode, **kwargs):
             "The 'schema' argument is not supported with the 'pyogrio' engine."
         )
 
-    if mode != "w":
-        raise ValueError(
-            "Only mode='w' is supported for now with the 'pyogrio' engine."
-        )
+    if mode == "a":
+        kwargs["append"] = True
 
     if crs is not None:
         raise ValueError("Passing 'crs' it not supported with the 'pyogrio' engine.")
diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py
index bbbc6222..3ba1125f 100644
--- a/geopandas/io/sql.py
+++ b/geopandas/io/sql.py
@@ -31,7 +31,10 @@ def _get_conn(conn_or_engine):
     from sqlalchemy.engine.base import Engine, Connection
 
     if isinstance(conn_or_engine, Connection):
-        with conn_or_engine.begin():
+        if not conn_or_engine.in_transaction():
+            with conn_or_engine.begin():
+                yield conn_or_engine
+        else:
             yield conn_or_engine
     elif isinstance(conn_or_engine, Engine):
         with conn_or_engine.begin() as conn:
@@ -395,6 +398,7 @@ def _write_postgis(
     """
     try:
         from geoalchemy2 import Geometry
+        from sqlalchemy import text
     except ImportError:
         raise ImportError("'to_postgis()' requires geoalchemy2 package.")
 
@@ -431,8 +435,10 @@ def _write_postgis(
             # Only check SRID if table exists
             if connection.dialect.has_table(connection, name, schema):
                 target_srid = connection.execute(
-                    "SELECT Find_SRID('{schema}', '{table}', '{geom_col}');".format(
-                        schema=schema_name, table=name, geom_col=geom_name
+                    text(
+                        "SELECT Find_SRID('{schema}', '{table}', '{geom_col}');".format(
+                            schema=schema_name, table=name, geom_col=geom_name
+                        )
                     )
                 ).fetchone()[0]
 
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 6f055416..7334f4c3 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -45,32 +45,66 @@ class BaseSpatialIndex:
         raise NotImplementedError
 
     def query(self, geometry, predicate=None, sort=False):
-        """Return the index of all geometries in the tree with extents that
-        intersect the envelope of the input geometry.
-
-        When using the ``rtree`` package, this is not a vectorized function.
-        If speed is important, please use PyGEOS.
+        """
+        Return the integer indices of all combinations of each input geometry
+        and tree geometries where the bounding box of each input geometry
+        intersects the bounding box of a tree geometry.
+
+        If the input geometry is a scalar, this returns an array of shape (n, ) with
+        the indices of the matching tree geometries.  If the input geometry is an
+        array_like, this returns an array with shape (2,n) where the subarrays
+        correspond to the indices of the input geometries and indices of the
+        tree geometries associated with each.  To generate an array of pairs of
+        input geometry index and tree geometry index, simply transpose the
+        result.
+
+        If a predicate is provided, the tree geometries are first queried based
+        on the bounding box of the input geometry and then are further filtered
+        to those that meet the predicate when comparing the input geometry to
+        the tree geometry: ``predicate(geometry, tree_geometry)``.
+
+        Bounding boxes are limited to two dimensions and are axis-aligned
+        (equivalent to the ``bounds`` property of a geometry); any Z values
+        present in input geometries are ignored when querying the tree.
+
+        Any input geometry that is None or empty will never match geometries in
+        the tree.
 
         Parameters
         ----------
-        geometry : shapely geometry
-            A single shapely geometry to query against the spatial index.
+        geometry : shapely.Geometry or array-like of geometries \
+(numpy.ndarray, GeoSeries, GeometryArray)
+            A single shapely geometry or array of geometries to query against
+            the spatial index. For array-like, accepts both GeoPandas geometry
+            iterables (GeoSeries, GeometryArray) or a numpy array of Shapely
+            or PyGEOS geometries.
         predicate : {None, "contains", "contains_properly", "covered_by", "covers", \
 "crosses", "intersects", "overlaps", "touches", "within"}, optional
-            If predicate is provided, the input geometry is
-            tested using the predicate function against each item
-            in the tree whose extent intersects the envelope of the
-            input geometry: predicate(input_geometry, tree_geometry).
-            If possible, prepared geometries are used to help
-            speed up the predicate operation.
+            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 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, the results will be sorted in ascending order.
-            If False, results are often sorted but there is no guarantee.
+            If True, the results will be sorted in ascending order. In case
+            of 2D array, the result is sorted lexicographically using the
+            geometries' indexes as the primary key and the sindex's indexes
+            as the secondary key.
+            If False, no additional sorting is applied (results are often
+            sorted but there is no guarantee).
 
         Returns
         -------
-        matches : ndarray of shape (n_results, )
-            Integer indices for matching geometries from the spatial index.
+        ndarray with shape (n,) if geometry is a scalar
+            Integer indices for matching geometries from the spatial index
+            tree geometries.
+
+        OR
+
+        ndarray with shape (2, n) if geometry is an array_like
+            The first subarray contains input geometry integer indices.
+            The second subarray contains tree geometry integer indices.
 
         Examples
         --------
@@ -89,16 +123,45 @@ class BaseSpatialIndex:
         9    POINT (9.00000 9.00000)
         dtype: geometry
 
+        Querying the tree with a scalar geometry:
+
         >>> s.sindex.query(box(1, 1, 3, 3))
         array([1, 2, 3])
 
         >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
         array([2])
+
+        Querying the tree with an array of geometries:
+
+        >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
+        >>> s2
+        0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
+        1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
+        dtype: geometry
+
+        >>> s.sindex.query(s2)
+        array([[0, 0, 0, 1, 1],
+               [2, 3, 4, 5, 6]])
+
+        >>> s.sindex.query(s2, predicate="contains")
+        array([[0],
+               [3]])
+
+        Notes
+        -----
+        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 overlapping bounding boxes or
+        optional predicate are returned.
         """
         raise NotImplementedError
 
     def query_bulk(self, geometry, predicate=None, sort=False):
         """
+        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.
@@ -451,6 +514,19 @@ if compat.HAS_RTREE:
                     )
                 )
 
+            if hasattr(geometry, "__array__") and not isinstance(
+                geometry, BaseGeometry
+            ):
+                # Iterates over geometry, applying func.
+                tree_index = []
+                input_geometry_index = []
+
+                for i, geo in enumerate(geometry):
+                    res = self.query(geo, predicate=predicate, sort=sort)
+                    tree_index.extend(res)
+                    input_geometry_index.extend([i] * len(res))
+                return np.vstack([input_geometry_index, tree_index])
+
             # handle empty / invalid geometries
             if geometry is None:
                 # return an empty integer array, similar to pygeos.STRtree.query.
@@ -524,15 +600,13 @@ if compat.HAS_RTREE:
 
         @doc(BaseSpatialIndex.query_bulk)
         def query_bulk(self, geometry, predicate=None, sort=False):
-            # Iterates over geometry, applying func.
-            tree_index = []
-            input_geometry_index = []
-
-            for i, geo in enumerate(geometry):
-                res = self.query(geo, predicate=predicate, sort=sort)
-                tree_index.extend(res)
-                input_geometry_index.extend([i] * len(res))
-            return np.vstack([input_geometry_index, tree_index])
+            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)
 
         def nearest(self, coordinates, num_results=1, objects=False):
             """
@@ -690,15 +764,26 @@ if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS:
                     )
                 )
 
-            if isinstance(geometry, BaseGeometry):
-                geometry = array._shapely_to_geom(geometry)
+            geometry = self._as_geometry_array(geometry)
 
-            matches = self._tree.query(geometry=geometry, predicate=predicate)
+            if compat.USE_SHAPELY_20:
+                indices = self._tree.query(geometry, predicate=predicate)
+            else:
+                if isinstance(geometry, np.ndarray):
+                    indices = self._tree.query_bulk(geometry, predicate=predicate)
+                else:
+                    indices = self._tree.query(geometry, predicate=predicate)
 
             if sort:
-                return np.sort(matches)
+                if indices.ndim == 1:
+                    return np.sort(indices)
+                else:
+                    # sort by first array (geometry) and then second (tree)
+                    geo_idx, tree_idx = indices
+                    sort_indexer = np.lexsort((tree_idx, geo_idx))
+                    return np.vstack((geo_idx[sort_indexer], tree_idx[sort_indexer]))
 
-            return matches
+            return indices
 
         @staticmethod
         def _as_geometry_array(geometry):
@@ -715,6 +800,11 @@ if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS:
             np.ndarray
                 A numpy array of pygeos geometries.
             """
+            # to ensure pygeos.Geometry as input is treated the same as shapely
+            # geometrie. TODO can be removed when we remove pygeos support
+            if isinstance(geometry, mod.Geometry):
+                geometry = array._geom_to_shapely(geometry)
+
             if isinstance(geometry, np.ndarray):
                 return array.from_shapely(geometry)._data
             elif isinstance(geometry, geoseries.GeoSeries):
@@ -739,27 +829,13 @@ if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS:
 
         @doc(BaseSpatialIndex.query_bulk)
         def query_bulk(self, geometry, predicate=None, sort=False):
-            if predicate not in self.valid_query_predicates:
-                raise ValueError(
-                    "Got `predicate` = `{}`, `predicate` must be one of {}".format(
-                        predicate, self.valid_query_predicates
-                    )
-                )
-
-            geometry = self._as_geometry_array(geometry)
-
-            if compat.USE_SHAPELY_20:
-                res = self._tree.query(geometry, predicate)
-            else:
-                res = self._tree.query_bulk(geometry, predicate)
-
-            if sort:
-                # sort by first array (geometry) and then second (tree)
-                geo_res, tree_res = res
-                indexing = np.lexsort((tree_res, geo_res))
-                return np.vstack((geo_res[indexing], tree_res[indexing]))
-
-            return res
+            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)
 
         @doc(BaseSpatialIndex.nearest)
         def nearest(
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index 674072b8..5c1dec9a 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -27,7 +27,7 @@ def _overlay_intersection(df1, df2):
     Overlay Intersection operation used in overlay function
     """
     # Spatial Index to create intersections
-    idx1, idx2 = df2.sindex.query_bulk(df1.geometry, predicate="intersects", sort=True)
+    idx1, idx2 = df2.sindex.query(df1.geometry, predicate="intersects", sort=True)
     # Create pairs of geometries in both dataframes to be intersected
     if idx1.size > 0 and idx2.size > 0:
         left = df1.geometry.take(idx1)
@@ -77,7 +77,7 @@ def _overlay_difference(df1, df2):
     Overlay Difference operation used in overlay function
     """
     # spatial index query to find intersections
-    idx1, idx2 = df2.sindex.query_bulk(df1.geometry, predicate="intersects", sort=True)
+    idx1, idx2 = df2.sindex.query(df1.geometry, predicate="intersects", sort=True)
     idx1_unique, idx1_unique_indices = np.unique(idx1, return_index=True)
     idx2_split = np.split(idx2, idx1_unique_indices[1:])
     sidx = [
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index fcb679a1..892444f3 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -217,7 +217,7 @@ def _geom_predicate_query(left_df, right_df, predicate):
             input_geoms = left_df.geometry
 
     if sindex:
-        l_idx, r_idx = sindex.query_bulk(input_geoms, predicate=predicate, sort=False)
+        l_idx, r_idx = sindex.query(input_geoms, predicate=predicate, sort=False)
         indices = pd.DataFrame({"_key_left": l_idx, "_key_right": r_idx})
     else:
         # when sindex is empty / has no valid geometries


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index f51c004e..29d7a4b0 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -107,7 +107,7 @@ jobs:
           pytest -v -r a -n auto --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/
 
       - name: Test with PostGIS
-        if: contains(matrix.env, '39-pd13-conda-forge.yaml') && contains(matrix.os, 'ubuntu')
+        if: (contains(matrix.env, '39-pd13-conda-forge.yaml') || contains(matrix.env, '311-latest-conda-forge.yaml')) && contains(matrix.os, 'ubuntu')
         env:
           PGUSER: postgres
           PGPASSWORD: postgres
diff --git a/ci/envs/311-latest-conda-forge.yaml b/ci/envs/311-latest-conda-forge.yaml
index dd3e0a95..748afa76 100644
--- a/ci/envs/311-latest-conda-forge.yaml
+++ b/ci/envs/311-latest-conda-forge.yaml
@@ -4,7 +4,7 @@ channels:
 dependencies:
   - python=3.11
   # required
-  - pandas
+  # - pandas
   - shapely
   - fiona
   - pyproj
@@ -26,10 +26,13 @@ dependencies:
   - geopy
   # installed in tests.yaml, because not available on windows
   # - postgis
-  - SQLalchemy<2
+  - SQLalchemy>=2
   - psycopg2
   - libspatialite
   - geoalchemy2
   - pyarrow
   # doctest testing
   - pytest-doctestplus
+  - pip
+  - pip:
+    - pandas==2.0.0rc1
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index bf2647df..f4f08855 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -221,6 +221,8 @@ def test_to_file_datetime(tmpdir, driver, ext, time, engine):
             df["b"].dt.tz_convert(pytz.utc), df_read["b"].dt.tz_convert(pytz.utc)
         )
     else:
+        if engine == "pyogrio" and PANDAS_GE_20:
+            df["b"] = df["b"].astype("datetime64[ms]")
         assert_series_equal(df["b"], df_read["b"])
 
 
@@ -475,13 +477,8 @@ def test_to_file_with_duplicate_columns(tmpdir, engine):
 @pytest.mark.parametrize("driver,ext", driver_ext_pairs)
 def test_append_file(tmpdir, df_nybb, df_null, driver, ext, engine):
     """Test to_file with append mode and from_file"""
-    skip_pyogrio_not_supported(engine)
-    from fiona import supported_drivers
-
     tempfilename = os.path.join(str(tmpdir), "boros" + ext)
     driver = driver if driver else _detect_driver(tempfilename)
-    if "a" not in supported_drivers[driver]:
-        return None
 
     df_nybb.to_file(tempfilename, driver=driver, engine=engine)
     df_nybb.to_file(tempfilename, mode="a", driver=driver, engine=engine)
@@ -492,6 +489,15 @@ def test_append_file(tmpdir, df_nybb, df_null, driver, ext, engine):
     expected = pd.concat([df_nybb] * 2, ignore_index=True)
     assert_geodataframe_equal(df, expected, check_less_precise=True)
 
+    if engine == "pyogrio":
+        # for pyogrio also ensure append=True works
+        tempfilename = os.path.join(str(tmpdir), "boros2" + ext)
+        df_nybb.to_file(tempfilename, driver=driver, engine=engine)
+        df_nybb.to_file(tempfilename, append=True, driver=driver, engine=engine)
+        # Read layer back in
+        df = GeoDataFrame.from_file(tempfilename, engine=engine)
+        assert len(df) == (len(df_nybb) * 2)
+
     # Write layer with null geometry out to file
     tempfilename = os.path.join(str(tmpdir), "null_geom" + ext)
     df_null.to_file(tempfilename, driver=driver, engine=engine)
@@ -504,6 +510,12 @@ def test_append_file(tmpdir, df_nybb, df_null, driver, ext, engine):
     assert_geodataframe_equal(df, expected, check_less_precise=True)
 
 
+def test_mode_unsupported(tmpdir, df_nybb, engine):
+    tempfilename = os.path.join(str(tmpdir), "data.shp")
+    with pytest.raises(ValueError, match="'mode' should be one of 'w' or 'a'"):
+        df_nybb.to_file(tempfilename, mode="r", engine=engine)
+
+
 @pytest.mark.parametrize("driver,ext", driver_ext_pairs)
 def test_empty_crs(tmpdir, driver, ext, engine):
     """Test handling of undefined CRS with GPKG driver (GH #1975)."""
diff --git a/geopandas/io/tests/test_sql.py b/geopandas/io/tests/test_sql.py
index 360fefa6..25ca56b4 100644
--- a/geopandas/io/tests/test_sql.py
+++ b/geopandas/io/tests/test_sql.py
@@ -5,6 +5,7 @@ configuration. postGIS tests require a test database to have been setup;
 see geopandas.tests.util for more information.
 """
 import os
+import warnings
 
 import pandas as pd
 
@@ -15,6 +16,13 @@ from geopandas.io.sql import _get_conn as get_conn, _write_postgis as write_post
 from geopandas.tests.util import create_postgis, create_spatialite, validate_boro_df
 import pytest
 
+try:
+    from sqlalchemy import text
+except ImportError:
+    # Avoid local imports for text in all sqlalchemy tests
+    # all tests using text use engine_postgis, which ensures sqlalchemy is available
+    text = str
+
 
 @pytest.fixture
 def df_nybb():
@@ -43,8 +51,11 @@ def connection_postgis():
         )
     except OperationalError:
         pytest.skip("Cannot connect with postgresql database")
-
-    yield con
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", message="pandas only supports SQLAlchemy connectable.*"
+        )
+        yield con
     con.close()
 
 
@@ -73,7 +84,7 @@ def engine_postgis():
                 port=port,
             )
         )
-        con.begin()
+        con.connect()
     except Exception:
         pytest.skip("Cannot connect with postgresql database")
 
@@ -115,11 +126,15 @@ def drop_table_if_exists(conn_or_engine, table):
     sqlalchemy = pytest.importorskip("sqlalchemy")
 
     if sqlalchemy.inspect(conn_or_engine).has_table(table):
-        metadata = sqlalchemy.MetaData(conn_or_engine)
-        metadata.reflect()
+        metadata = sqlalchemy.MetaData()
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore", message="Did not recognize type 'geometry' of column.*"
+            )
+            metadata.reflect(conn_or_engine)
         table = metadata.tables.get(table)
         if table is not None:
-            table.drop(checkfirst=True)
+            table.drop(conn_or_engine, checkfirst=True)
 
 
 @pytest.fixture
@@ -340,7 +355,7 @@ class TestIO:
         # Write to db
         write_postgis(df_nybb, con=engine, name=table, if_exists="fail")
         # Validate
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         validate_boro_df(df)
 
@@ -355,7 +370,7 @@ class TestIO:
         # Write to db
         write_postgis(df_nybb, con=engine, name=table, if_exists="fail")
         # Validate
-        sql = 'SELECT * FROM "{table}";'.format(table=table)
+        sql = text('SELECT * FROM "{table}";'.format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         validate_boro_df(df)
 
@@ -370,7 +385,7 @@ class TestIO:
             # Write to db
             write_postgis(df_nybb, con=con, name=table, if_exists="fail")
             # Validate
-            sql = "SELECT * FROM {table};".format(table=table)
+            sql = text("SELECT * FROM {table};".format(table=table))
             df = read_postgis(sql, con, geom_col="geometry")
             validate_boro_df(df)
 
@@ -406,7 +421,7 @@ class TestIO:
         # Overwrite
         write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
         # Validate
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         validate_boro_df(df)
 
@@ -423,7 +438,7 @@ class TestIO:
         write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
         write_postgis(df_nybb, con=engine, name=table, if_exists="append")
         # Validate
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         new_rows, new_cols = df.shape
 
@@ -449,13 +464,16 @@ class TestIO:
         # Write to db
         df_nybb = df_nybb
         df_nybb.crs = None
-        write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
+        with pytest.warns(UserWarning, match="Could not parse CRS from the GeoDataF"):
+            write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
         # Validate that srid is -1
-        target_srid = engine.execute(
+        sql = text(
             "SELECT Find_SRID('{schema}', '{table}', '{geom_col}');".format(
                 schema="public", table=table, geom_col="geometry"
             )
-        ).fetchone()[0]
+        )
+        with engine.connect() as conn:
+            target_srid = conn.execute(sql).fetchone()[0]
         assert target_srid == 0, "SRID should be 0, found %s" % target_srid
 
     def test_write_postgis_with_esri_authority(self, engine_postgis, df_nybb):
@@ -471,11 +489,13 @@ class TestIO:
         df_nybb_esri = df_nybb.to_crs("ESRI:102003")
         write_postgis(df_nybb_esri, con=engine, name=table, if_exists="replace")
         # Validate that srid is 102003
-        target_srid = engine.execute(
+        sql = text(
             "SELECT Find_SRID('{schema}', '{table}', '{geom_col}');".format(
                 schema="public", table=table, geom_col="geometry"
             )
-        ).fetchone()[0]
+        )
+        with engine.connect() as conn:
+            target_srid = conn.execute(sql).fetchone()[0]
         assert target_srid == 102003, "SRID should be 102003, found %s" % target_srid
 
     def test_write_postgis_geometry_collection(
@@ -491,11 +511,14 @@ class TestIO:
         write_postgis(df_geom_collection, con=engine, name=table, if_exists="replace")
 
         # Validate geometry type
-        sql = "SELECT DISTINCT(GeometryType(geometry)) FROM {table} ORDER BY 1;".format(
-            table=table
+        sql = text(
+            "SELECT DISTINCT(GeometryType(geometry)) FROM {table} ORDER BY 1;".format(
+                table=table
+            )
         )
-        geom_type = engine.execute(sql).fetchone()[0]
-        sql = "SELECT * FROM {table};".format(table=table)
+        with engine.connect() as conn:
+            geom_type = conn.execute(sql).fetchone()[0]
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
 
         assert geom_type.upper() == "GEOMETRYCOLLECTION"
@@ -516,10 +539,13 @@ class TestIO:
         )
 
         # Validate geometry type
-        sql = "SELECT DISTINCT GeometryType(geometry) FROM {table} ORDER BY 1;".format(
-            table=table
+        sql = text(
+            "SELECT DISTINCT GeometryType(geometry) FROM {table} ORDER BY 1;".format(
+                table=table
+            )
         )
-        res = engine.execute(sql).fetchall()
+        with engine.connect() as conn:
+            res = conn.execute(sql).fetchall()
         assert res[0][0].upper() == "LINESTRING"
         assert res[1][0].upper() == "MULTILINESTRING"
         assert res[2][0].upper() == "POINT"
@@ -535,10 +561,13 @@ class TestIO:
         write_postgis(df_linear_ring, con=engine, name=table, if_exists="replace")
 
         # Validate geometry type
-        sql = "SELECT DISTINCT(GeometryType(geometry)) FROM {table} ORDER BY 1;".format(
-            table=table
+        sql = text(
+            "SELECT DISTINCT(GeometryType(geometry)) FROM {table} ORDER BY 1;".format(
+                table=table
+            )
         )
-        geom_type = engine.execute(sql).fetchone()[0]
+        with engine.connect() as conn:
+            geom_type = conn.execute(sql).fetchone()[0]
 
         assert geom_type.upper() == "LINESTRING"
 
@@ -558,15 +587,19 @@ class TestIO:
             chunksize=1,
         )
         # Validate row count
-        sql = "SELECT COUNT(geometry) FROM {table};".format(table=table)
-        row_cnt = engine.execute(sql).fetchone()[0]
+        sql = text("SELECT COUNT(geometry) FROM {table};".format(table=table))
+        with engine.connect() as conn:
+            row_cnt = conn.execute(sql).fetchone()[0]
         assert row_cnt == 3
 
         # Validate geometry type
-        sql = "SELECT DISTINCT GeometryType(geometry) FROM {table} ORDER BY 1;".format(
-            table=table
+        sql = text(
+            "SELECT DISTINCT GeometryType(geometry) FROM {table} ORDER BY 1;".format(
+                table=table
+            )
         )
-        res = engine.execute(sql).fetchall()
+        with engine.connect() as conn:
+            res = conn.execute(sql).fetchall()
         assert res[0][0].upper() == "LINESTRING"
         assert res[1][0].upper() == "MULTILINESTRING"
         assert res[2][0].upper() == "POINT"
@@ -579,15 +612,16 @@ class TestIO:
 
         table = "nybb"
         schema_to_use = "test"
-        sql = "CREATE SCHEMA IF NOT EXISTS {schema};".format(schema=schema_to_use)
-        engine.execute(sql)
+        sql = text("CREATE SCHEMA IF NOT EXISTS {schema};".format(schema=schema_to_use))
+        with engine.begin() as conn:
+            conn.execute(sql)
 
         write_postgis(
             df_nybb, con=engine, name=table, if_exists="replace", schema=schema_to_use
         )
         # Validate
-        sql = "SELECT * FROM {schema}.{table};".format(
-            schema=schema_to_use, table=table
+        sql = text(
+            "SELECT * FROM {schema}.{table};".format(schema=schema_to_use, table=table)
         )
 
         df = read_postgis(sql, engine, geom_col="geometry")
@@ -603,16 +637,19 @@ class TestIO:
 
         table = "nybb"
         schema_to_use = "test"
-        sql = "CREATE SCHEMA IF NOT EXISTS {schema};".format(schema=schema_to_use)
-        engine.execute(sql)
+        sql = text("CREATE SCHEMA IF NOT EXISTS {schema};".format(schema=schema_to_use))
+        with engine.begin() as conn:
+            conn.execute(sql)
 
         try:
             write_postgis(
                 df_nybb, con=engine, name=table, if_exists="fail", schema=schema_to_use
             )
             # Validate
-            sql = "SELECT * FROM {schema}.{table};".format(
-                schema=schema_to_use, table=table
+            sql = text(
+                "SELECT * FROM {schema}.{table};".format(
+                    schema=schema_to_use, table=table
+                )
             )
 
             df = read_postgis(sql, engine, geom_col="geometry")
@@ -627,8 +664,8 @@ class TestIO:
             df_nybb, con=engine, name=table, if_exists="replace", schema=schema_to_use
         )
         # Validate
-        sql = "SELECT * FROM {schema}.{table};".format(
-            schema=schema_to_use, table=table
+        sql = text(
+            "SELECT * FROM {schema}.{table};".format(schema=schema_to_use, table=table)
         )
 
         df = read_postgis(sql, engine, geom_col="geometry")
@@ -645,7 +682,7 @@ class TestIO:
         write_postgis(df_3D_geoms, con=engine, name=table, if_exists="replace")
 
         # Check that all geometries have 3 dimensions
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         assert list(df.geometry.has_z) == [True, True, True]
 
@@ -661,7 +698,7 @@ class TestIO:
         write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
 
         # Check that the row order matches
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         assert df["BoroCode"].tolist() == correct_order
 
@@ -678,7 +715,7 @@ class TestIO:
         write_postgis(df_nybb, con=engine, name=table, if_exists="append")
 
         # Check that the row order matches
-        sql = "SELECT * FROM {table};".format(table=table)
+        sql = text("SELECT * FROM {table};".format(table=table))
         df = read_postgis(sql, engine, geom_col="geometry")
         validate_boro_df(df)
 
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 632db49e..73e8463e 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -542,6 +542,19 @@ class TestDataFrame:
         assert type(df2) is GeoDataFrame
         assert self.df.crs == df2.crs
 
+    def test_empty_copy(self):
+        # https://github.com/geopandas/geopandas/issues/2765
+        df = GeoDataFrame()
+        df2 = df.copy()
+        assert type(df2) is GeoDataFrame
+        df3 = df.copy(deep=True)
+        assert type(df3) is GeoDataFrame
+
+    def test_no_geom_copy(self):
+        df = GeoDataFrame(pd.DataFrame({"a": [1, 2, 3]}))
+        assert type(df) is GeoDataFrame
+        assert type(df.copy()) is GeoDataFrame
+
     def test_bool_index(self):
         # Find boros with 'B' in their name
         df = self.df[self.df["BoroName"].str.contains("B")]
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index bacc908d..7bebc84b 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -1466,3 +1466,32 @@ class TestGeomMethods:
             NotImplementedError, match="shapely >= 2.0 or PyGEOS are required"
         ):
             self.g11.get_coordinates()
+
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="minimum_bounding_radius not implemented for shapely<2",
+    )
+    def test_minimum_bounding_radius(self):
+        mbr_geoms = self.g1.minimum_bounding_radius()
+
+        assert_series_equal(
+            mbr_geoms,
+            Series([0.707106, 0.707106]),
+        )
+
+        mbr_lines = self.g5.minimum_bounding_radius()
+
+        assert_series_equal(
+            mbr_lines,
+            Series([0.707106, 0.707106]),
+        )
+
+    @pytest.mark.skipif(
+        (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="minimum_bounding_radius not implemented for shapely<2",
+    )
+    def test_minimium_bounding_radius_not(self):
+        with pytest.raises(
+            NotImplementedError, match="shapely >= 2.0 or PyGEOS is required"
+        ):
+            self.g1.minimum_bounding_radius()
diff --git a/geopandas/tests/test_op_output_types.py b/geopandas/tests/test_op_output_types.py
index 52d048fe..ea2ed500 100644
--- a/geopandas/tests/test_op_output_types.py
+++ b/geopandas/tests/test_op_output_types.py
@@ -62,9 +62,7 @@ def _check_metadata_gs(gs, name="geometry", crs=crs_wgs):
     assert gs.crs == crs
 
 
-def assert_object(
-    result, expected_type, geo_name="geometry", crs=crs_wgs, check_none_name=False
-):
+def assert_object(result, expected_type, geo_name="geometry", crs=crs_wgs):
     """
     Helper method to make tests easier to read. Checks result is of the expected
     type. If result is a GeoDataFrame or GeoSeries, checks geo_name
@@ -75,36 +73,55 @@ def assert_object(
     assert type(result) is expected_type
 
     if expected_type == GeoDataFrame:
-        if geo_name is not None:
-            _check_metadata_gdf(result, geo_name=geo_name, crs=crs)
-        else:
-            if check_none_name:  # TODO this is awkward
-                assert result._geometry_column_name is None
-
-            if result._geometry_column_name is None:
-                msg = (
-                    "You are calling a geospatial method on the GeoDataFrame, "
-                    "but the active"
-                )
-            else:
-                msg = (
-                    "You are calling a geospatial method on the GeoDataFrame, but "
-                    r"the active geometry column \("
-                    rf"'{result._geometry_column_name}'\) is not present"
-                )
-            with pytest.raises(AttributeError, match=msg):
-                result.geometry.name  # be explicit that geometry is invalid here
+        assert geo_name is not None
+        _check_metadata_gdf(result, geo_name=geo_name, crs=crs)
+
     elif expected_type == GeoSeries:
         _check_metadata_gs(result, name=geo_name, crs=crs)
 
 
+def assert_obj_no_active_geo_col(result, expected_type, geo_colname=None):
+    """
+    Helper method to make tests easier to read. Checks result is of the expected
+    type. Asserts that accessing result.geometry.name raises, corresponding to
+    _geometry_column_name being in an invalid state
+    (either None, or a column no longer present)
+    This amounts to testing the assertion raised (geometry column is unset, vs
+    old geometry column is missing)
+
+    We assert that _geometry_column_name = int_geo_colname
+
+    """
+    if expected_type == GeoDataFrame:
+        if geo_colname is None:
+            assert result._geometry_column_name is None
+        else:
+            assert geo_colname == result._geometry_column_name
+
+        if result._geometry_column_name is None:
+            msg = (
+                "You are calling a geospatial method on the GeoDataFrame, "
+                "but the active"
+            )
+        else:
+            msg = (
+                "You are calling a geospatial method on the GeoDataFrame, but "
+                r"the active geometry column \("
+                rf"'{result._geometry_column_name}'\) is not present"
+            )
+        with pytest.raises(AttributeError, match=msg):
+            result.geometry.name  # be explicit that geometry is invalid here
+    else:
+        raise NotImplementedError()
+
+
 def test_getitem(df):
     geo_name = df.geometry.name
     assert_object(df[["value1", "value2"]], pd.DataFrame)
     assert_object(df[[geo_name, "geometry2"]], GeoDataFrame, geo_name)
     assert_object(df[[geo_name]], GeoDataFrame, geo_name)
-    assert_object(df[["geometry2", "value1"]], GeoDataFrame, None)
-    assert_object(df[["geometry2"]], GeoDataFrame, None)
+    assert_obj_no_active_geo_col(df[["geometry2", "value1"]], GeoDataFrame)
+    assert_obj_no_active_geo_col(df[["geometry2"]], GeoDataFrame)
     assert_object(df[["value1"]], pd.DataFrame)
     # Series
     assert_object(df[geo_name], GeoSeries, geo_name)
@@ -117,8 +134,13 @@ def test_loc(df):
     assert_object(df.loc[:, ["value1", "value2"]], pd.DataFrame)
     assert_object(df.loc[:, [geo_name, "geometry2"]], GeoDataFrame, geo_name)
     assert_object(df.loc[:, [geo_name]], GeoDataFrame, geo_name)
-    assert_object(df.loc[:, ["geometry2", "value1"]], GeoDataFrame, None)
-    assert_object(df.loc[:, ["geometry2"]], GeoDataFrame, None)
+    # These two are inconsistent with getitem, active geom col dropped,
+    # but other geometry columns present
+    assert_obj_no_active_geo_col(
+        df.loc[:, ["geometry2", "value1"]], GeoDataFrame, geo_name
+    )
+    assert_obj_no_active_geo_col(df.loc[:, ["geometry2"]], GeoDataFrame, geo_name)
+    # #####
     assert_object(df.loc[:, ["value1"]], pd.DataFrame)
     # Series
     assert_object(df.loc[:, geo_name], GeoSeries, geo_name)
@@ -131,8 +153,11 @@ def test_iloc(df):
     assert_object(df.iloc[:, 0:2], pd.DataFrame)
     assert_object(df.iloc[:, 2:4], GeoDataFrame, geo_name)
     assert_object(df.iloc[:, [2]], GeoDataFrame, geo_name)
-    assert_object(df.iloc[:, [3, 0]], GeoDataFrame, None)
-    assert_object(df.iloc[:, [3]], GeoDataFrame, None)
+    # These two are inconsistent with getitem, active geom col dropped,
+    # but other geometry columns present
+    assert_obj_no_active_geo_col(df.iloc[:, [3, 0]], GeoDataFrame, geo_name)
+    assert_obj_no_active_geo_col(df.iloc[:, [3]], GeoDataFrame, geo_name)
+    # #####
     assert_object(df.iloc[:, [0]], pd.DataFrame)
     # Series
     assert_object(df.iloc[:, 2], GeoSeries, geo_name)
@@ -164,8 +189,15 @@ def test_reindex(df):
     assert_object(df.reindex(columns=[geo_name, "geometry2"]), GeoDataFrame, geo_name)
     assert_object(df.reindex(columns=[geo_name]), GeoDataFrame, geo_name)
     assert_object(df.reindex(columns=["new_col", geo_name]), GeoDataFrame, geo_name)
-    assert_object(df.reindex(columns=["geometry2", "value1"]), GeoDataFrame, None)
-    assert_object(df.reindex(columns=["geometry2"]), GeoDataFrame, None)
+    # These two are inconsistent with getitem, active geom col dropped,
+    # but other geometry columns present
+    assert_obj_no_active_geo_col(
+        df.reindex(columns=["geometry2", "value1"]), GeoDataFrame, geo_name
+    )
+    assert_obj_no_active_geo_col(
+        df.reindex(columns=["geometry2"]), GeoDataFrame, geo_name
+    )
+    # #####
     assert_object(df.reindex(columns=["value1"]), pd.DataFrame)
 
     # reindexing the rows always preserves the GeoDataFrame
@@ -184,8 +216,15 @@ def test_drop(df):
     assert_object(df.drop(columns=["value1", "value2"]), GeoDataFrame, geo_name)
     cols = ["value1", "value2", "geometry2"]
     assert_object(df.drop(columns=cols), GeoDataFrame, geo_name)
-    assert_object(df.drop(columns=[geo_name, "value2"]), GeoDataFrame, None)
-    assert_object(df.drop(columns=["value1", "value2", geo_name]), GeoDataFrame, None)
+    # These two are inconsistent with getitem, active geom col dropped,
+    # but other geometry columns present
+    assert_obj_no_active_geo_col(
+        df.drop(columns=[geo_name, "value2"]), GeoDataFrame, geo_name
+    )
+    assert_obj_no_active_geo_col(
+        df.drop(columns=["value1", "value2", geo_name]), GeoDataFrame, geo_name
+    )
+    # #####
     assert_object(df.drop(columns=["geometry2", "value2", geo_name]), pd.DataFrame)
 
 
@@ -199,8 +238,13 @@ def test_apply(df):
     assert_object(df[["value1", "value2"]].apply(identity), pd.DataFrame)
     assert_object(df[[geo_name, "geometry2"]].apply(identity), GeoDataFrame, geo_name)
     assert_object(df[[geo_name]].apply(identity), GeoDataFrame, geo_name)
-    assert_object(df[["geometry2", "value1"]].apply(identity), GeoDataFrame, None, None)
-    assert_object(df[["geometry2"]].apply(identity), GeoDataFrame, None, None)
+    expected_geo_col_name = None if compat.PANDAS_GE_14 else "geometry"
+    assert_obj_no_active_geo_col(
+        df[["geometry2", "value1"]].apply(identity), GeoDataFrame, expected_geo_col_name
+    )
+    assert_obj_no_active_geo_col(
+        df[["geometry2"]].apply(identity), GeoDataFrame, expected_geo_col_name
+    )
     assert_object(df[["value1"]].apply(identity), pd.DataFrame)
 
     # axis = 0, Series
@@ -231,7 +275,10 @@ def test_apply_axis1_secondary_geo_cols(df):
     def identity(x):
         return x
 
-    assert_object(df[["geometry2"]].apply(identity, axis=1), GeoDataFrame, None, None)
+    expected_geo_col_name = None if compat.PANDAS_GE_14 else "geometry"
+    assert_obj_no_active_geo_col(
+        df[["geometry2"]].apply(identity, axis=1), GeoDataFrame, expected_geo_col_name
+    )
 
 
 def test_expanddim_in_apply():
@@ -257,9 +304,9 @@ def test_expandim_in_groupby_aggregate_multiple_funcs():
 
     grouped = s.groupby([0, 1, 0])
     agg = grouped.agg([total_area, union])
-    assert_object(agg, GeoDataFrame, None, None, check_none_name=True)
+    assert_obj_no_active_geo_col(agg, GeoDataFrame, geo_colname=None)
     result = grouped.agg([union, total_area])
-    assert_object(result, GeoDataFrame, None, None, check_none_name=True)
+    assert_obj_no_active_geo_col(result, GeoDataFrame, geo_colname=None)
     assert_object(grouped.agg([total_area, total_area]), pd.DataFrame)
     assert_object(grouped.agg([total_area]), pd.DataFrame)
 
@@ -272,17 +319,15 @@ def test_expanddim_in_unstack():
         index=pd.MultiIndex.from_tuples([("A", "a"), ("A", "b"), ("B", "a")]),
     )
     unstack = s.unstack()
-    assert_object(unstack, GeoDataFrame, None, None, False)
-
-    if compat.PANDAS_GE_12:
-        assert unstack._geometry_column_name is None
-    else:  # pandas GH37369, unstack doesn't call finalize
-        assert unstack._geometry_column_name == "geometry"
+    expected_geo_col_name = None if compat.PANDAS_GE_12 else "geometry"
+    assert_obj_no_active_geo_col(
+        unstack, GeoDataFrame, geo_colname=expected_geo_col_name
+    )
 
     # https://github.com/geopandas/geopandas/issues/2486
     s.name = "geometry"
     unstack = s.unstack()
-    assert_object(unstack, GeoDataFrame, None, None)
+    assert_obj_no_active_geo_col(unstack, GeoDataFrame, expected_geo_col_name)
 
 
 # indexing /  constructor_sliced tests
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index 0cc2e42b..e60127cb 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -571,7 +571,7 @@ class TestPygeosInterface:
         test_geom = geopandas.GeoSeries(
             [box(*geom) for geom in test_geom], index=range(len(test_geom))
         )
-        res = self.df.sindex.query_bulk(test_geom, predicate=predicate)
+        res = self.df.sindex.query(test_geom, predicate=predicate)
         assert_array_equal(res, expected)
 
     @pytest.mark.parametrize(
@@ -591,14 +591,14 @@ class TestPygeosInterface:
         # determine if it should use shapely or pygeos geometry objects
         # note: for this test, test_geoms (note plural) is a list already
         test_geoms = geopandas.GeoSeries(test_geoms, index=range(len(test_geoms)))
-        res = self.df.sindex.query_bulk(test_geoms)
+        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."""
         test_array = np.array([], dtype=object)
         expected_value = [[], []]
-        res = self.df.sindex.query_bulk(test_array)
+        res = self.df.sindex.query(test_array)
         assert_array_equal(res, expected_value)
 
     def test_query_bulk_invalid_input_geometry(self):
@@ -607,7 +607,7 @@ class TestPygeosInterface:
         """
         test_array = "notanarray"
         with pytest.raises(TypeError):
-            self.df.sindex.query_bulk(test_array)
+            self.df.sindex.query(test_array)
 
     def test_query_bulk_invalid_predicate(self):
         """Tests the `query_bulk` method with invalid predicates."""
@@ -619,7 +619,7 @@ class TestPygeosInterface:
         test_geom = geopandas.GeoSeries([box(*test_geom_bounds)], index=["0"])
 
         with pytest.raises(ValueError):
-            self.df.sindex.query_bulk(test_geom.geometry, predicate=test_predicate)
+            self.df.sindex.query(test_geom.geometry, predicate=test_predicate)
 
     @pytest.mark.parametrize(
         "predicate, test_geom, expected",
@@ -638,21 +638,21 @@ class TestPygeosInterface:
         test_geom = geopandas.GeoSeries([box(*test_geom)], index=["0"])
 
         # test GeoSeries
-        res = self.df.sindex.query_bulk(test_geom, predicate=predicate)
+        res = self.df.sindex.query(test_geom, predicate=predicate)
         assert_array_equal(res, expected)
 
         # test GeometryArray
-        res = self.df.sindex.query_bulk(test_geom.geometry, predicate=predicate)
+        res = self.df.sindex.query(test_geom.geometry, predicate=predicate)
         assert_array_equal(res, expected)
-        res = self.df.sindex.query_bulk(test_geom.geometry.values, predicate=predicate)
+        res = self.df.sindex.query(test_geom.geometry.values, predicate=predicate)
         assert_array_equal(res, expected)
 
         # test numpy array
-        res = self.df.sindex.query_bulk(
+        res = self.df.sindex.query(
             test_geom.geometry.values.to_numpy(), predicate=predicate
         )
         assert_array_equal(res, expected)
-        res = self.df.sindex.query_bulk(
+        res = self.df.sindex.query(
             test_geom.geometry.values.to_numpy(), predicate=predicate
         )
         assert_array_equal(res, expected)
@@ -686,7 +686,7 @@ class TestPygeosInterface:
         tree_df = geopandas.GeoDataFrame(geometry=tree_polys)
         test_df = geopandas.GeoDataFrame(geometry=test_polys)
 
-        res = tree_df.sindex.query_bulk(test_df.geometry, sort=sort)
+        res = tree_df.sindex.query(test_df.geometry, sort=sort)
 
         # asserting the same elements
         assert sorted(res[0]) == sorted(expected[0])
@@ -880,7 +880,7 @@ class TestPygeosInterface:
         world = read_file(datasets.get_path("naturalearth_lowres"))
         capitals = read_file(datasets.get_path("naturalearth_cities"))
 
-        res = world.sindex.query_bulk(capitals.geometry, predicate)
+        res = world.sindex.query(capitals.geometry, predicate)
         assert res.shape == expected_shape