Repository: geopandas/geopandas
Pull request #2314 (2022-01-21 11:54:11)

DEP: remove deprecated fallback to return pd.Series from GeoSeries constructor

martinfleis: Thanks!

I guess you should also remove `_SERIES_WARNING_MSG` on l.31 and that `catch_warnings` context in _geoseries_constructor_with_fallback.
m-richards: The other thing I wasn't sure about was if `__new__` should be moved to `__init__`, since there is no early returning with a pd.Series anymore.
martinfleis: > The other thing I wasn't sure about was if __new__ should be moved to __init__, since there is no early returning with a pd.Series anymore.

Probably yes, guessing from that comment in there. Can you try that to see if it causes some issues? It shouldn't. 

Patch:
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 0173a8cd..1965d115 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -69,6 +69,7 @@ Set-theoretic Methods
 .. autosummary::
    :toctree: api/
 
+   GeoSeries.clip_by_rect
    GeoSeries.difference
    GeoSeries.intersection
    GeoSeries.symmetric_difference
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 457a81fc..159e5781 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -11,6 +11,7 @@ import pandas as pd
 
 import shapely.geometry
 import shapely.geos
+import shapely.ops
 import shapely.wkb
 import shapely.wkt
 
@@ -728,6 +729,20 @@ def almost_equals(self, other, decimal):
 #
 
 
+def clip_by_rect(data, xmin, ymin, xmax, ymax):
+    if compat.USE_PYGEOS:
+        return pygeos.clip_by_rect(data, xmin, ymin, xmax, ymax)
+    else:
+        clipped_geometries = np.empty(len(data), dtype=object)
+        clipped_geometries[:] = [
+            shapely.ops.clip_by_rect(s, xmin, ymin, xmax, ymax)
+            if s is not None
+            else None
+            for s in data
+        ]
+        return clipped_geometries
+
+
 def difference(data, other):
     if compat.USE_PYGEOS:
         return _binary_method("difference", data, other)
diff --git a/geopandas/array.py b/geopandas/array.py
index 5f3f24b6..9b47d6a6 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -564,6 +564,11 @@ class GeometryArray(ExtensionArray):
     # Binary operations that return new geometries
     #
 
+    def clip_by_rect(self, xmin, ymin, xmax, ymax):
+        return GeometryArray(
+            vectorized.clip_by_rect(self.data, xmin, ymin, xmax, ymax), crs=self.crs
+        )
+
     def difference(self, other):
         return GeometryArray(
             self._binary_method("difference", self, other), crs=self.crs
diff --git a/geopandas/base.py b/geopandas/base.py
index 2edf6474..dab414ba 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2536,6 +2536,75 @@ GeometryCollection
         """
         return _binary_geo("intersection", self, other, align)
 
+    def clip_by_rect(self, xmin, ymin, xmax, ymax):
+        """Returns a ``GeoSeries`` of the portions of geometry within the given
+        rectangle.
+
+        Note that the results are not exactly equal to
+        :meth:`~GeoSeries.intersection()`. E.g. in edge cases,
+        :meth:`~GeoSeries.clip_by_rect()` will not return a point just touching the
+        rectangle. Check the examples section below for some of these exceptions.
+
+        The geometry is clipped in a fast but possibly dirty way. The output is not
+        guaranteed to be valid. No exceptions will be raised for topological errors.
+
+        Note: empty geometries or geometries that do not overlap with the specified
+        bounds will result in ``GEOMETRYCOLLECTION EMPTY``.
+
+        Parameters
+        ----------
+        xmin: float
+            Minimum x value of the rectangle
+        ymin: float
+            Minimum y value of the rectangle
+        xmax: float
+            Maximum x value of the rectangle
+        ymax: float
+            Maximum y value of the rectangle
+
+        Returns
+        -------
+        GeoSeries
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         LineString([(0, 0), (2, 2)]),
+        ...         LineString([(2, 0), (0, 2)]),
+        ...         Point(0, 1),
+        ...     ],
+        ...     crs=3857,
+        ... )
+        >>> bounds = (0, 0, 1, 1)
+        >>> s
+        0    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        1    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        2                LINESTRING (0.000 0.000, 2.000 2.000)
+        3                LINESTRING (2.000 0.000, 0.000 2.000)
+        4                                  POINT (0.000 1.000)
+        dtype: geometry
+        >>> s.clip_by_rect(*bounds)
+        0    POLYGON ((0.000 0.000, 0.000 1.000, 1.000 1.00...
+        1    POLYGON ((0.000 0.000, 0.000 1.000, 1.000 1.00...
+        2                LINESTRING (0.000 0.000, 1.000 1.000)
+        3                             GEOMETRYCOLLECTION EMPTY
+        4                             GEOMETRYCOLLECTION EMPTY
+        dtype: geometry
+
+        See also
+        --------
+        GeoSeries.intersection
+        """
+        from .geoseries import GeoSeries
+
+        geometry_array = GeometryArray(self.geometry.values)
+        clipped_geometry = geometry_array.clip_by_rect(xmin, ymin, xmax, ymax)
+        return GeoSeries(clipped_geometry.data, index=self.index, crs=self.crs)
+
     #
     # Other operations
     #
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 3b4b5163..bfd298d4 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -189,11 +189,10 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
             self.set_geometry(geometry, inplace=True)
 
         if geometry is None and crs:
-            warnings.warn(
-                "Assigning CRS to a GeoDataFrame without a geometry column is now "
-                "deprecated and will not be supported in the future.",
-                FutureWarning,
-                stacklevel=2,
+            raise ValueError(
+                "Assigning CRS to a GeoDataFrame without a geometry column is not "
+                "supported. Supply geometry using the 'geometry=' keyword argument, "
+                "or by providing a DataFrame with column name 'geometry'",
             )
 
     def __setattr__(self, attr, val):
@@ -312,7 +311,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         geo_column_name = self._geometry_column_name
         if isinstance(col, (Series, list, np.ndarray, GeometryArray)):
             level = col
-        elif hasattr(col, "ndim") and col.ndim != 1:
+        elif hasattr(col, "ndim") and col.ndim > 1:
             raise ValueError("Must pass array with one dimension only.")
         else:
             try:
@@ -442,13 +441,11 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
     def crs(self, value):
         """Sets the value of the crs"""
         if self._geometry_column_name not in self:
-            warnings.warn(
-                "Assigning CRS to a GeoDataFrame without a geometry column is now "
-                "deprecated and will not be supported in the future.",
-                FutureWarning,
-                stacklevel=4,
+            raise ValueError(
+                "Assigning CRS to a GeoDataFrame without a geometry column is not "
+                "supported. Use GeoDataFrame.set_geometry to set the active "
+                "geometry column.",
             )
-            self._crs = None if not value else CRS.from_user_input(value)
         else:
             if hasattr(self.geometry.values, "crs"):
                 self.geometry.values.crs = value
@@ -1470,6 +1467,11 @@ individually so that features may have different properties
                     f"Please ensure this column from the first DataFrame is not "
                     f"repeated."
                 )
+        elif method == "unstack":
+            # unstack adds multiindex columns and reshapes data.
+            # it never makes sense to retain geometry column
+            self._geometry_column_name = None
+            self._crs = None
         return self
 
     def dissolve(
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 2452bedd..395a6a8b 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -3,7 +3,7 @@ import warnings
 
 import numpy as np
 import pandas as pd
-from pandas import Series, MultiIndex
+from pandas import Series, MultiIndex, DataFrame
 from pandas.core.internals import SingleBlockManager
 
 from pyproj import CRS
@@ -28,12 +28,6 @@ from .array import (
 from .base import is_geometry_type
 
 
-_SERIES_WARNING_MSG = """\
-    You are passing non-geometry data to the GeoSeries constructor. Currently,
-    it falls back to returning a pandas Series. But in the future, we will start
-    to raise a TypeError instead."""
-
-
 def _geoseries_constructor_with_fallback(data=None, index=None, crs=None, **kwargs):
     """
     A flexible constructor for GeoSeries._constructor, which needs to be able
@@ -41,18 +35,40 @@ def _geoseries_constructor_with_fallback(data=None, index=None, crs=None, **kwar
     geometries)
     """
     try:
-        with warnings.catch_warnings():
-            warnings.filterwarnings(
-                "ignore",
-                message=_SERIES_WARNING_MSG,
-                category=FutureWarning,
-                module="geopandas[.*]",
-            )
-            return GeoSeries(data=data, index=index, crs=crs, **kwargs)
+        return GeoSeries(data=data, index=index, crs=crs, **kwargs)
     except TypeError:
         return Series(data=data, index=index, **kwargs)
 
 
+def _geoseries_expanddim(data=None, *args, **kwargs):
+    from geopandas import GeoDataFrame
+
+    # pd.Series._constructor_expanddim == pd.DataFrame
+    df = pd.DataFrame(data, *args, **kwargs)
+    geo_col_name = None
+    if isinstance(data, GeoSeries):
+        # pandas default column name is 0, keep convention
+        geo_col_name = data.name if data.name is not None else 0
+
+    if df.shape[1] == 1:
+        geo_col_name = df.columns[0]
+
+    if (df.dtypes == "geometry").sum() > 0:
+        if geo_col_name is None or not is_geometry_type(df[geo_col_name]):
+            df = GeoDataFrame(df)
+            df._geometry_column_name = None
+        else:
+            df = df.set_geometry(geo_col_name)
+
+    return df
+
+
+# pd.concat (pandas/core/reshape/concat.py) requires this for the
+# concatenation of series since pandas 1.1
+# (https://github.com/pandas-dev/pandas/commit/f9e4c8c84bcef987973f2624cc2932394c171c8c)
+_geoseries_expanddim._get_axis_number = DataFrame._get_axis_number
+
+
 class GeoSeries(GeoPandasBase, Series):
     """
     A Series object designed to store shapely geometry objects.
@@ -133,10 +149,7 @@ class GeoSeries(GeoPandasBase, Series):
 
     _metadata = ["name"]
 
-    def __new__(cls, data=None, index=None, crs=None, **kwargs):
-        # we need to use __new__ because we want to return Series instance
-        # instead of GeoSeries instance in case of non-geometry data
-
+    def __init__(self, data=None, index=None, crs=None, **kwargs):
         if hasattr(data, "crs") and crs:
             if not data.crs:
                 # make a copy to avoid setting CRS to passed GeometryArray
@@ -163,12 +176,11 @@ class GeoSeries(GeoPandasBase, Series):
                     values = data.blocks[0].values
                     block = ExtensionBlock(values, slice(0, len(values), 1), ndim=1)
                     data = SingleBlockManager([block], data.axes[0], fastpath=True)
-                self = super(GeoSeries, cls).__new__(cls)
-                super(GeoSeries, self).__init__(data, index=index, **kwargs)
-                self.crs = getattr(self.values, "crs", crs)
-                return self
-            warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
-            return Series(data, index=index, **kwargs)
+            else:
+                raise TypeError(
+                    "Non geometry data passed to GeoSeries constructor, "
+                    f"received data of dtype '{data.blocks[0].dtype}'"
+                )
 
         if isinstance(data, BaseGeometry):
             # fix problem for scalar geometries passed, ensure the list of
@@ -198,28 +210,24 @@ class GeoSeries(GeoPandasBase, Series):
                     # pd.Series with empty data gives float64 for older pandas versions
                     s = s.astype(object)
                 else:
-                    warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
-                    return s
+                    raise TypeError(
+                        "Non geometry data passed to GeoSeries constructor, "
+                        f"received data of dtype '{s.dtype}'"
+                    )
             # try to convert to GeometryArray, if fails return plain Series
             try:
                 data = from_shapely(s.values, crs)
             except TypeError:
-                warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
-                return s
+                raise TypeError(
+                    "Non geometry data passed to GeoSeries constructor, "
+                    f"received data of dtype '{s.dtype}'"
+                )
             index = s.index
             name = s.name
 
-        self = super(GeoSeries, cls).__new__(cls)
-        super(GeoSeries, self).__init__(data, index=index, name=name, **kwargs)
-
+        super().__init__(data, index=index, name=name, **kwargs)
         if not self.crs:
             self.crs = crs
-        return self
-
-    def __init__(self, *args, **kwargs):
-        # need to overwrite Series init to prevent calling it for GeoSeries
-        # (doesn't know crs, all work is already done above)
-        pass
 
     def append(self, *args, **kwargs):
         return self._wrapped_pandas_method("append", *args, **kwargs)
@@ -587,9 +595,7 @@ class GeoSeries(GeoPandasBase, Series):
 
     @property
     def _constructor_expanddim(self):
-        from geopandas import GeoDataFrame
-
-        return GeoDataFrame
+        return _geoseries_expanddim
 
     def _wrapped_pandas_method(self, mtd, *args, **kwargs):
         """Wrap a generic pandas method to ensure it returns a GeoSeries"""
diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py
index 9bd0f1d1..17e7c404 100644
--- a/geopandas/io/sql.py
+++ b/geopandas/io/sql.py
@@ -230,7 +230,7 @@ def _get_geometry_type(gdf):
 
     # Check for 3D-coordinates
     if any(gdf.geometry.has_z):
-        target_geom_type = target_geom_type + "Z"
+        target_geom_type += "Z"
 
     return target_geom_type, has_curve
 
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index babb67a5..15133eff 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -139,7 +139,7 @@ def _overlay_union(df1, df2):
     # keep geometry column last
     columns = list(dfunion.columns)
     columns.remove("geometry")
-    columns = columns + ["geometry"]
+    columns.append("geometry")
     return dfunion.reindex(columns=columns)
 
 
diff --git a/setup.py b/setup.py
index 3ed6d865..ca1d0ec1 100644
--- a/setup.py
+++ b/setup.py
@@ -58,6 +58,9 @@ setup(
     author="GeoPandas contributors",
     author_email="kjordahl@alum.mit.edu",
     url="http://geopandas.org",
+    project_urls={
+        "Source": "https://github.com/geopandas/geopandas",
+    },
     long_description=LONG_DESCRIPTION,
     packages=[
         "geopandas",


Test patch:
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index 2b995e47..04a4270a 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -242,9 +242,13 @@ class TestGeometryArrayCRS:
         assert df.geometry.crs == self.wgs
         assert df.geometry.values.crs == self.wgs
 
-        df = GeoDataFrame(self.geoms, columns=["geom"], crs=27700)
-        assert df.crs == self.osgb
-        df = df.set_geometry("geom")
+        with pytest.raises(ValueError, match="Assigning CRS to a GeoDataFrame without"):
+            GeoDataFrame(self.geoms, columns=["geom"], crs=27700)
+        with pytest.raises(ValueError, match="Assigning CRS to a GeoDataFrame without"):
+            GeoDataFrame(crs=27700)
+
+        df = GeoDataFrame(self.geoms, columns=["geom"])
+        df = df.set_geometry("geom", crs=27700)
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
@@ -256,14 +260,8 @@ class TestGeometryArrayCRS:
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
 
-        df = GeoDataFrame(crs=27700)
-        df = df.set_geometry(self.geoms)
-        assert df.crs == self.osgb
-        assert df.geometry.crs == self.osgb
-        assert df.geometry.values.crs == self.osgb
-
         # new geometry with set CRS has priority over GDF CRS
-        df = GeoDataFrame(crs=27700)
+        df = GeoDataFrame(geometry=self.geoms, crs=27700)
         df = df.set_geometry(self.geoms, crs=4326)
         assert df.crs == self.wgs
         assert df.geometry.crs == self.wgs
@@ -341,14 +339,18 @@ class TestGeometryArrayCRS:
         "scalar", [None, Point(0, 0), LineString([(0, 0), (1, 1)])]
     )
     def test_scalar(self, scalar):
-        with pytest.warns(FutureWarning):
-            df = GeoDataFrame()
-            df.crs = 4326
+        df = GeoDataFrame()
         df["geometry"] = scalar
+        df.crs = 4326
         assert df.crs == self.wgs
         assert df.geometry.crs == self.wgs
         assert df.geometry.values.crs == self.wgs
 
+    def test_crs_with_no_geom_fails(self):
+        with pytest.raises(ValueError, match="Assigning CRS to a GeoDataFrame without"):
+            df = GeoDataFrame()
+            df.crs = 4326
+
     def test_read_file(self):
         nybb_filename = datasets.get_path("nybb")
         df = read_file(nybb_filename)
@@ -579,21 +581,6 @@ class TestGeometryArrayCRS:
         assert merged.geom.values.crs == self.osgb
         assert merged.crs == self.osgb
 
-    # CRS should be assigned to geometry
-    def test_deprecation(self):
-        with pytest.warns(FutureWarning):
-            df = GeoDataFrame([], crs=27700)
-
-        # https://github.com/geopandas/geopandas/issues/1548
-        # ensure we still have converted the crs value to a CRS object
-        assert isinstance(df.crs, pyproj.CRS)
-
-        with pytest.warns(FutureWarning):
-            df = GeoDataFrame([])
-            df.crs = 27700
-
-        assert isinstance(df.crs, pyproj.CRS)
-
     # make sure that geometry column from list has CRS (__setitem__)
     def test_setitem_geometry(self):
         arr = from_shapely(self.geoms, crs=27700)
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 630021d9..a5ec71bb 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -334,6 +334,11 @@ class TestDataFrame:
         assert isinstance(result, GeoDataFrame)
         assert isinstance(result.index, pd.DatetimeIndex)
 
+    def test_set_geometry_np_int(self):
+        self.df.loc[:, 0] = self.df.geometry
+        df = self.df.set_geometry(np.int64(0))
+        assert df.geometry.name == 0
+
     def test_get_geometry_invalid(self):
         df = GeoDataFrame()
         df["geom"] = self.df.geometry
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 05475a96..05e3d5c1 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -85,8 +85,12 @@ class TestGeomMethods:
         self.g8 = GeoSeries([self.t1, self.t5])
         self.empty = GeoSeries([])
         self.all_none = GeoSeries([None, None])
+        self.all_geometry_collection_empty = GeoSeries(
+            [GeometryCollection([]), GeometryCollection([])]
+        )
         self.empty_poly = Polygon()
         self.g9 = GeoSeries(self.g0, index=range(1, 8))
+        self.g10 = GeoSeries([self.t1, self.t4])
 
         # Crossed lines
         self.l3 = LineString([(0, 0), (1, 1)])
@@ -254,6 +258,15 @@ class TestGeomMethods:
             assert len(self.g0.intersection(self.g9, align=True) == 8)
         assert len(self.g0.intersection(self.g9, align=False) == 7)
 
+    def test_clip_by_rect(self):
+        self._test_binary_topological(
+            "clip_by_rect", self.g1, self.g10, *self.sq.bounds
+        )
+        # self.g1 and self.t3.bounds do not intersect
+        self._test_binary_topological(
+            "clip_by_rect", self.all_geometry_collection_empty, self.g1, *self.t3.bounds
+        )
+
     def test_union_series(self):
         self._test_binary_topological("union", self.sq, self.g1, self.g2)
 
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 9182f03d..c511c3f3 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -7,7 +7,7 @@ import tempfile
 import numpy as np
 from numpy.testing import assert_array_equal
 import pandas as pd
-from pandas.util.testing import assert_index_equal
+from pandas.testing import assert_index_equal
 
 from pyproj import CRS
 from shapely.geometry import (
@@ -464,21 +464,15 @@ class TestConstructor:
             for x in gs:
                 assert x.equals(g)
 
-    def test_no_geometries_fallback(self):
-        with pytest.warns(FutureWarning):
-            s = GeoSeries([True, False, True])
-        assert not isinstance(s, GeoSeries)
-        assert type(s) == pd.Series
+    def test_non_geometry_raises(self):
+        with pytest.raises(TypeError, match="Non geometry data passed to GeoSeries"):
+            GeoSeries([True, False, True])
 
-        with pytest.warns(FutureWarning):
-            s = GeoSeries(["a", "b", "c"])
-        assert not isinstance(s, GeoSeries)
-        assert type(s) == pd.Series
+        with pytest.raises(TypeError, match="Non geometry data passed to GeoSeries"):
+            GeoSeries(["a", "b", "c"])
 
-        with pytest.warns(FutureWarning):
-            s = GeoSeries([[1, 2], [3, 4]])
-        assert not isinstance(s, GeoSeries)
-        assert type(s) == pd.Series
+        with pytest.raises(TypeError, match="Non geometry data passed to GeoSeries"):
+            GeoSeries([[1, 2], [3, 4]])
 
     def test_empty(self):
         s = GeoSeries([])
@@ -494,7 +488,6 @@ class TestConstructor:
     def test_empty_array(self):
         # with empty data that have an explicit dtype, we use the fallback or
         # not depending on the dtype
-        arr = np.array([], dtype="bool")
 
         # dtypes that can never hold geometry-like data
         for arr in [
@@ -504,10 +497,10 @@ class TestConstructor:
             # this gets converted to object dtype by pandas
             # np.array([], dtype="str"),
         ]:
-            with pytest.warns(FutureWarning):
-                s = GeoSeries(arr)
-            assert not isinstance(s, GeoSeries)
-            assert type(s) == pd.Series
+            with pytest.raises(
+                TypeError, match="Non geometry data passed to GeoSeries"
+            ):
+                GeoSeries(arr)
 
         # dtypes that can potentially hold geometry-like data (object) or
         # can come from empty data (float64)
@@ -538,13 +531,38 @@ class TestConstructor:
         assert s.index is g.index
 
     # GH 1216
-    def test_expanddim(self):
+    @pytest.mark.parametrize("name", [None, "geometry", "Points"])
+    @pytest.mark.parametrize("crs", [None, "epsg:4326"])
+    def test_reset_index(self, name, crs):
         s = GeoSeries(
-            [MultiPoint([(0, 0), (1, 1)]), MultiPoint([(2, 2), (3, 3), (4, 4)])]
+            [MultiPoint([(0, 0), (1, 1)]), MultiPoint([(2, 2), (3, 3), (4, 4)])],
+            name=name,
+            crs=crs,
         )
         s = s.explode(index_parts=True)
         df = s.reset_index()
         assert type(df) == GeoDataFrame
+        # name None -> 0, otherwise name preserved
+        assert df.geometry.name == (name if name is not None else 0)
+        assert df.crs == s.crs
+
+    @pytest.mark.parametrize("name", [None, "geometry", "Points"])
+    @pytest.mark.parametrize("crs", [None, "epsg:4326"])
+    def test_to_frame(self, name, crs):
+        s = GeoSeries([Point(0, 0), Point(1, 1)], name=name, crs=crs)
+        df = s.to_frame()
+        assert type(df) == GeoDataFrame
+        # name None -> 0, otherwise name preserved
+        expected_name = name if name is not None else 0
+        assert df.geometry.name == expected_name
+        assert df._geometry_column_name == expected_name
+        assert df.crs == s.crs
+
+        # if name is provided to to_frame, it should override
+        df2 = s.to_frame(name="geom")
+        assert type(df) == GeoDataFrame
+        assert df2.geometry.name == "geom"
+        assert df2.crs == s.crs
 
     def test_explode_without_multiindex(self):
         s = GeoSeries(
@@ -560,7 +578,6 @@ class TestConstructor:
         )
         s = s.explode(ignore_index=True)
         expected_index = pd.Index(range(len(s)))
-        print(expected_index)
         assert_index_equal(s.index, expected_index)
 
         # index_parts is ignored if ignore_index=True
diff --git a/geopandas/tests/test_merge.py b/geopandas/tests/test_merge.py
index a0495ab5..a29fbae8 100644
--- a/geopandas/tests/test_merge.py
+++ b/geopandas/tests/test_merge.py
@@ -1,10 +1,12 @@
 import pandas as pd
 import pytest
 from geopandas.testing import assert_geodataframe_equal
+from pandas.testing import assert_index_equal
 
 from shapely.geometry import Point
 
 from geopandas import GeoDataFrame, GeoSeries
+from geopandas import _compat as compat
 
 
 class TestMerging:
@@ -97,3 +99,23 @@ class TestMerging:
         res3 = pd.concat([df2.set_crs("epsg:4326"), self.gdf], axis=1)
         # check metadata comes from first df
         self._check_metadata(res3, geometry_column_name="geom", crs="epsg:4326")
+
+    @pytest.mark.xfail(
+        not compat.PANDAS_GE_11,
+        reason="pandas <=1.0 hard codes concat([GeoSeries, GeoSeries]) -> "
+        "DataFrame or Union[DataFrame, SparseDataFrame] in 0.25",
+    )
+    def test_concat_axis1_geoseries(self):
+        gseries2 = GeoSeries([Point(i, i) for i in range(3, 6)], crs="epsg:4326")
+        result = pd.concat([gseries2, self.gseries], axis=1)
+        # Note this is not consistent with concat([gdf, gdf], axis=1) where the
+        # left metadata is set on the result. This is deliberate for now.
+        assert type(result) is GeoDataFrame
+        self._check_metadata(result, geometry_column_name=None, crs=None)
+        assert_index_equal(pd.Index([0, 1]), result.columns)
+
+        gseries2.name = "foo"
+        result2 = pd.concat([gseries2, self.gseries], axis=1)
+        assert type(result2) is GeoDataFrame
+        self._check_metadata(result2, geometry_column_name=None, crs=None)
+        assert_index_equal(pd.Index(["foo", 0]), result2.columns)
diff --git a/geopandas/tests/test_op_output_types.py b/geopandas/tests/test_op_output_types.py
index 291abfbd..5e9b5d96 100644
--- a/geopandas/tests/test_op_output_types.py
+++ b/geopandas/tests/test_op_output_types.py
@@ -1,6 +1,7 @@
 import pandas as pd
 import pyproj
 import pytest
+import geopandas._compat as compat
 
 from shapely.geometry import Point
 
@@ -46,7 +47,9 @@ 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):
+def assert_object(
+    result, expected_type, geo_name="geometry", crs=crs_wgs, check_none_name=False
+):
     """
     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
@@ -60,6 +63,9 @@ def assert_object(result, expected_type, geo_name="geometry", crs=crs_wgs):
         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, "
@@ -135,18 +141,11 @@ def test_squeeze(df):
 
 def test_to_frame(df):
     geo_name = df.geometry.name
-    # TODO this reflects current behaviour, but we should fix
-    #  GeoSeries._constructor_expanddim so this doesn't happen
     res1 = df[geo_name].to_frame()
-    if geo_name == "geometry":  # -> this should be doable for any geo_name
-        assert_object(res1, GeoDataFrame, geo_name)
-    assert res1._geometry_column_name == "geometry"  # -> should be geo_name
+    assert_object(res1, GeoDataFrame, geo_name, crs=df[geo_name].crs)
 
     res2 = df["geometry2"].to_frame()
-    assert type(res2) is GeoDataFrame
-    assert res2._geometry_column_name == "geometry"  # -> should be geometry2
-    assert res2.crs is None  # -> should be crs_osgb
-    # also res2.geometry should not crash because geometry isn't set
+    assert_object(res2, GeoDataFrame, "geometry2", crs=crs_osgb)
 
     res3 = df["value1"].to_frame()
     assert_object(res3, pd.DataFrame)
@@ -235,3 +234,57 @@ def test_apply_axis1_secondary_geo_cols(df):
         return x
 
     assert_object(df[["geometry2"]].apply(identity, axis=1), GeoDataFrame, None, None)
+
+
+def test_expanddim_in_apply():
+    # https://github.com/geopandas/geopandas/pull/2296#issuecomment-1021966443
+    s = GeoSeries.from_xy([0, 1], [0, 1])
+    result = s.apply(lambda x: pd.Series([x.x, x.y]))
+    assert_object(result, pd.DataFrame)
+
+
+@pytest.mark.xfail(
+    not compat.PANDAS_GE_11,
+    reason="pandas <1.1 don't preserve subclass through groupby ops",  # Pandas GH33884
+)
+def test_expandim_in_groupby_aggregate_multiple_funcs():
+    # https://github.com/geopandas/geopandas/pull/2296#issuecomment-1021966443
+    # There are two calls to _constructor_expanddim here
+    # SeriesGroupBy._aggregate_multiple_funcs() and
+    # SeriesGroupBy._wrap_series_output() len(output) > 1
+
+    s = GeoSeries.from_xy([0, 1, 2], [0, 1, 3])
+
+    def union(s):
+        return s.unary_union
+
+    def total_area(s):
+        return s.area.sum()
+
+    grouped = s.groupby([0, 1, 0])
+    agg = grouped.agg([total_area, union])
+    assert_object(agg, GeoDataFrame, None, None, check_none_name=True)
+    result = grouped.agg([union, total_area])
+    assert_object(result, GeoDataFrame, None, None, check_none_name=True)
+    assert_object(grouped.agg([total_area, total_area]), pd.DataFrame)
+    assert_object(grouped.agg([total_area]), pd.DataFrame)
+
+
+@pytest.mark.xfail(
+    not compat.PANDAS_GE_11,
+    reason="pandas <1.1 uses concat([Series]) in unstack",  # Pandas GH33356
+)
+def test_expanddim_in_unstack():
+    # https://github.com/geopandas/geopandas/pull/2296#issuecomment-1021966443
+    s = GeoSeries.from_xy(
+        [0, 1, 2],
+        [0, 1, 3],
+        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"