Repository: geopandas/geopandas

Issue:
BUG: GeoSeries.to_frame when name is not 'geometry' doesn't set geometry column
- [x] I have checked that this issue has not already been reported.

- [x] I have confirmed this bug exists on the latest version of geopandas.

- [x] (optional) I have confirmed this bug exists on the master branch of geopandas.

---

**Note**: Please read [this guide](https://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports) detailing how to provide the necessary information for us to reproduce your bug.

#### Code Sample, a copy-pastable example

```python
In [2]: import geopandas

In [3]:  df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_cities"))

In [4]: srs = df['geometry']

In [5]: type(srs)
Out[5]: geopandas.geoseries.GeoSeries

In [6]: srs.name = 'foo'

In [7]: srs.to_frame().geometry
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-7-ee808c91f778> in <module>
----> 1 srs.to_frame().geometry

C:\MyPrograms\Miniconda3\envs\gpd_latest_conda_forge\lib\site-packages\pandas\core\generic.py in __getattr__(self, name)
   5476         ):
   5477             return self[name]
-> 5478         return object.__getattribute__(self, name)
   5479
   5480     def __setattr__(self, name: str, value) -> None:

c:\data\mydocuments\code\geopandas_test\geopandas-m-richards\geopandas\geodataframe.py in _get_geometry(self)
    196     def _get_geometry(self):
    197         if self._geometry_column_name not in self:
--> 198             raise AttributeError(
    199                 "No geometry data set yet (expected in"
    200                 " column '%s'.)" % self._geometry_column_name

AttributeError: No geometry data set yet (expected in column 'geometry'.)

In [8]: srs.name = 'geometry'

In [9]: srs.to_frame().geometry
Out[9]:
0        POINT (12.45339 41.90328)
1        POINT (12.44177 43.93610)
2         POINT (9.51667 47.13372)
3         POINT (6.13000 49.61166)
4        POINT (158.14997 6.91664)
                  ...
197      POINT (31.24802 30.05191)
198     POINT (139.74946 35.68696)
199       POINT (2.33139 48.86864)
200    POINT (-70.66899 -33.44807)
201      POINT (103.85387 1.29498)
Name: geometry, Length: 202, dtype: geometry
```

#### Problem description
`GeoSeries.to_frame()` should produce valid geometry, regardless of whether the series name is geometry.

#### Expected Output
As in [9]. (Edit: also `srs.to_frame().crs` should be preserved and not None)

I also think that if `GeoSeries.name = None` it might make more sense to produce a GeoDataFrame with a single column "geometry" rather than '0' that is default. 

Yet to look into properly, but I think this probably could be fixed by expanding `GeoSeries._constructor_expanddim`, similar to some of the other stuff I've been looking at.

#### Output of ``geopandas.show_versions()``

<details>

SYSTEM INFO
-----------
python     : 3.8.10 | packaged by conda-forge | (default, May 11 2021, 06:25:23) [MSC v.1916 64 bit (AMD64)]
executable : C:\MyPrograms\Miniconda3\envs\gpd_latest_conda_forge\python.exe
machine    : Windows-10-10.0.19043-SP0

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : None
GEOS lib   : None
GDAL       : 3.3.1
GDAL data dir: None
PROJ       : 8.0.1
PROJ data dir: C:\MyPrograms\Miniconda3\envs\gpd_latest_conda_forge\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.9.0+59.g4536b92
pandas     : 1.3.1
fiona      : 1.8.20
numpy      : 1.21.1
shapely    : 1.7.1
rtree      : 0.9.7
pyproj     : 3.1.0
matplotlib : 3.4.2
mapclassify: 2.4.3
geopy      : 2.2.0
psycopg2   : 2.9.1 (dt dec pq3 ext lo64)
geoalchemy2: None
pyarrow    : 5.0.0
pygeos     : 0.10.1

</details>


Discussion:
jorisvandenbossche: 
Some other cases where `_constructor_expanddim` is used (in addition to `to_frame`):


```
In [28]: s = geopandas.GeoSeries.from_xy([0, 1], [0, 1])

In [29]: s.apply(lambda x: pd.Series([x.x, x.y]))
Out[29]: 
     0    1
0  0.0  0.0
1  1.0  1.0

In [30]: type(_)
Out[30]: geopandas.geodataframe.GeoDataFrame   # <--- wrong
```

or 

```
In [37]: s = geopandas.GeoSeries.from_xy([0, 1, 2], [0, 1, 3])

In [38]: def union(s):
    ...:     return s.unary_union

In [39]: def total_area(s):
    ...:     return s.area.sum()

In [40]: s.groupby([0, 1, 0]).agg([total_area, union])
Out[40]: 
   total_area                                          union
0         0.0  MULTIPOINT (0.00000 0.00000, 2.00000 3.00000)
1         0.0                        POINT (1.00000 1.00000)

In [41]: s.groupby([0, 1, 0]).agg([union, total_area])
Out[41]: 
                                           union  total_area
0  MULTIPOINT (0.00000 0.00000, 2.00000 3.00000)         0.0
1                        POINT (1.00000 1.00000)         0.0

In [42]: type(_)
Out[42]: geopandas.geodataframe.GeoDataFrame
```

or 

```
In [47]: s = geopandas.GeoSeries.from_xy(
    ...:     [0, 1, 2],
    ...:     [0, 1, 3],
    ...:     index=pd.MultiIndex.from_tuples([("A", "a"), ("A", "b"), ("B", "a")]),
    ...: )

In [48]: s
Out[48]: 
A  a    POINT (0.00000 0.00000)
   b    POINT (1.00000 1.00000)
B  a    POINT (2.00000 3.00000)
dtype: geometry

In [49]: s.unstack()
Out[49]: 
                         a                        b
A  POINT (0.00000 0.00000)  POINT (1.00000 1.00000)
B  POINT (2.00000 3.00000)                     None

In [50]: type(_)
Out[50]: geopandas.geodataframe.GeoDataFrame
```

And in many cases the `_geometry_column_name` will also be wrong. 

It are a bit contrived examples, but probably could to check to ensure that the new logic in `_constructor_expandim` doesn't break on those.
jorisvandenbossche: Can you also add some tests based on the examples at https://github.com/geopandas/geopandas/pull/2296#issuecomment-1021966443? (unless those are already covered)
m-richards: > Can you also add some tests based on the examples at [#2296 (comment)](https://github.com/geopandas/geopandas/pull/2296#issuecomment-1021966443)? (unless those are already covered)

Of course! I saw these and then forgot about them going through the review comments, that's my bad. 
```python
s = geopandas.GeoSeries.from_xy([0, 1], [0, 1])
s.apply(lambda x: pd.Series([x.x, x.y]))
```
reveals a big missed case where expanddim can recieve something that's not going to be a geodataframe. I'm now going to double check other usages in pandas. It's a bit fraught trying to alter these really powerful internal methods but I hope it will be worth it in the end.

I've just gone through and done this, the only other two places are in NDFrame.align(a_geo_series) which I can't work out a meaningful test case - I don't use align, but the GeoSeries is interpreted as row rather than a column, so I can't see how that could be used in a meaningful way when applied to a GeoSeries. The other place is in the StringMethods accessor, and .str only works on string dtype columns anyway, so isn't applicable (This is perhaps why you didn't mention these in the list of other cases above).

There's a bit of ugly conditional logic now, so I think I will see if it's possible to clean that up.


m-richards: Was hoping to get this all green but 
https://github.com/pandas-dev/pandas/blob/04e720f201e53ec1a0e9cf474df7d111fd386a31/pandas/core/groupby/generic.py#L314-L315
as a pattern is a problem - there's similar code in concat as well. 

This is only really a problem because I decided that if expanddim recieves a dict, it should set the active geometry column to be the leftmost named column, to be consistent with concat of dataframes. Given the dicussion around getitem, I think it actually makes more sense to leave the geometry column unset in this case and make the user set it (and I think that logic should hold for the dataframe case as well?). 

(This is the case I'm referring to)
```python
s = gpd.GeoSeries.from_xy([1,2], [3,4])
df =pd.concat([s,s], axis=1)
In [132]: df
Out[132]:
                         0                        1
0  POINT (1.00000 3.00000)  POINT (1.00000 3.00000)
1  POINT (2.00000 4.00000)  POINT (2.00000 4.00000)
In [131]: df._geometry_column_name # what should this be, for frames, it is the leftmost item in concat's geom col

```




m-richards: I think this is now in a reviewable state. The problem in my comment above is now resolved.

There is a separate question about whether `concat([a,b])` should actually preserve the left active geometry column or not (This PR mimics the current gdf case, but that could be removed / deferred by removing most of the changes in `__finalize__` - that's what I was going to do so that this could progress, and in stumbling through that I found a fix here.
jorisvandenbossche: > This is only really a problem because I decided that if expanddim recieves a dict, it should set the active geometry column to be the leftmost named column, to be consistent with concat of dataframes. Given the dicussion around getitem, I think it actually makes more sense to leave the geometry column unset in this case and make the user set it (and I think that logic should hold for the dataframe case as well?).

Yes, I think there is indeed an argument for keeping the geometry column name unset (certainly now we have better error messages around that). In any case, in interest of keeping this PR nicely scoped, let's maybe leave that for a separate issue / PR? (from your last comment I understand that this is mostly leaving out the changes in `__finalize__`?)
m-richards: 
> Yes, I think there is indeed an argument for keeping the geometry column name unset (certainly now we have better error messages around that). In any case, in interest of keeping this PR nicely scoped, let's maybe leave that for a separate issue / PR? (from your last comment I understand that this is mostly leaving out the changes in `__finalize__`?)

Yep that's right. Happy to leave it, probably makes sense in a broader discussion alongside the geodataframe case. 


jorisvandenbossche: I did another read through the code and tests and discussion, and added some comments. And thanks again for doing this! :)
jorisvandenbossche: > Was hoping to get this all green but https://github.com/pandas-dev/pandas/blob/04e720f201e53ec1a0e9cf474df7d111fd386a31/pandas/core/groupby/generic.py#L314-L315 as a pattern is a problem - there's similar code in concat as well.

So that pattern is like

```python
# creating a default DataFrame
result = obj._constructor_expanddim(result, ..)
# and then in a second step overriding the columns names
result.columns = Index([...])
```

That's indeed difficult to work with for us (because it means we can't necessarily rely on what is passed to `_constructor_expanddim` to infer the geometry column name). In theory we could override the `columns` setter to map old to new geometry column name, but not sure that's something we should get into.

Currently the above mainly happens for cases where there might be multiple geometry columns? (like the concat) If we for now decide to leave the geometry column name unset (`None`) in those cases, this is OK for now?
m-richards: > Currently the above mainly happens for cases where there might be multiple geometry columns? (like the concat) If we for now decide to leave the geometry column name unset (`None`) in those cases, this is OK for now?

This is no longer an issue here (this was the reason I gave up on #2080 and replaced it with this PR (I was trying to have "geometry" propagate as the first column name / default name in `to_frame` in a less than smart way).

The last sticking point in this PR was in `unstack` but that is now fixed because it calls finalize and I've fixed that there.

overriding the setitem of columns did occur to me as well, really don't want to have to do that if we don't have to - and I think we don't, so that's good.

Unless there's a version specific failure in the CI I think I'm done updating this in response to those comments. So hopefully that's the end of notification spam on this thread for you.

>I did another read through the code and tests and discussion, and added some comments. And thanks again for doing this! :)
Thanks for reviewing and the feedback, would have been a bit of a mess to review.  




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..db351c67 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
@@ -53,6 +53,35 @@ def _geoseries_constructor_with_fallback(data=None, index=None, crs=None, **kwar
         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.
@@ -587,9 +616,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..450fcc09 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 (
@@ -538,13 +538,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 +585,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"