Repository: geopandas/geopandas

Issue:
ENH: return GeoSeries when accessing other (non "active") geometry columns
Now that that we have the geometry dtype and ExtensionArray based GeometryArray, it is more transparent to have multiple geometry columns (which was already possible before, but they were all "object")

```
In [30]: df = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres')) 

In [31]: df['centroids'] = df.geometry.centroid   

In [32]: df.dtypes   
Out[32]: 
pop_est          int64
continent       object
name            object
iso_a3          object
gdp_md_est     float64
geometry      geometry
centroids     geometry
dtype: object
```

But, when accessing the column, we only box it into a GeoSeries if it was the "active" geometry column, and not for other "geometry"-dtype columns:

```
In [33]: type(df['geometry']) 
Out[33]: geopandas.geoseries.GeoSeries

In [34]: type(df['centroids'])  
Out[34]: pandas.core.series.Series
```

Now that we can check the dtype of the column for GeometryDtype, we can actually also return a GeoSeries for other geometry columns, enabling the same functionality (additional functions provided by the GeoSeries subclass) for those columns. 
(the GeoDataFrame level functions of course keep using the "active" geometry column, just as it is the case now)

Discussion:
martinfleis: I have slightly changed `assert_geoseries_equal` as this PR caused coverage decrease there. While making additional test to cover that I figured out that `check_dtype` was essentially ignored as it first checks for GeometryDtype and then it does `check_dtype`. So I fixed that.

Patch:
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 2dc663c0..617d98c6 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -6,7 +6,7 @@ from pandas import DataFrame, Series
 
 from shapely.geometry import mapping, shape
 
-from geopandas.array import GeometryArray, from_shapely
+from geopandas.array import GeometryArray, from_shapely, GeometryDtype
 from geopandas.base import GeoPandasBase, is_geometry_type
 from geopandas.geoseries import GeoSeries
 import geopandas.io
@@ -587,7 +587,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         """
         result = super(GeoDataFrame, self).__getitem__(key)
         geo_col = self._geometry_column_name
-        if isinstance(key, str) and key == geo_col:
+        if isinstance(result, Series) and isinstance(result.dtype, GeometryDtype):
             result.__class__ = GeoSeries
             result.crs = self.crs
             result._invalidate_sindex()


Test patch:
diff --git a/geopandas/testing.py b/geopandas/testing.py
index d3ed09b7..5c20b13e 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -61,7 +61,7 @@ def geom_almost_equals(this, that):
 def assert_geoseries_equal(
     left,
     right,
-    check_dtype=False,
+    check_dtype=True,
     check_index_type=False,
     check_series_type=True,
     check_less_precise=False,
@@ -92,16 +92,14 @@ def assert_geoseries_equal(
     """
     assert len(left) == len(right), "%d != %d" % (len(left), len(right))
 
-    msg = "dtype should be a GeometryDtype, got {0}"
-    assert isinstance(left.dtype, GeometryDtype), msg.format(left.dtype)
-    assert isinstance(right.dtype, GeometryDtype), msg.format(left.dtype)
+    if check_dtype:
+        msg = "dtype should be a GeometryDtype, got {0}"
+        assert isinstance(left.dtype, GeometryDtype), msg.format(left.dtype)
+        assert isinstance(right.dtype, GeometryDtype), msg.format(left.dtype)
 
     if check_index_type:
         assert isinstance(left.index, type(right.index))
 
-    if check_dtype:
-        assert left.dtype == right.dtype, "dtype: %s != %s" % (left.dtype, right.dtype)
-
     if check_series_type:
         assert isinstance(left, GeoSeries)
         assert isinstance(left, type(right))
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index a09871fb..9645e7e2 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -12,7 +12,7 @@ from shapely.geometry import Point
 
 import geopandas
 from geopandas import GeoDataFrame, GeoSeries, read_file
-from geopandas.array import GeometryArray, GeometryDtype
+from geopandas.array import GeometryArray, GeometryDtype, from_shapely
 
 from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
 from geopandas.tests.util import PACKAGE_DIR, connect, create_postgis, validate_boro_df
@@ -77,6 +77,12 @@ class TestDataFrame:
         assert not isinstance(df["geometry"], GeoSeries)
         assert isinstance(df["location"], GeoSeries)
 
+        df["buff"] = df.buffer(1)
+        assert isinstance(df["buff"], GeoSeries)
+
+        df["array"] = from_shapely([Point(x, y) for x, y in zip(range(5), range(5))])
+        assert isinstance(df["array"], GeoSeries)
+
         data["geometry"] = [Point(x + 1, y - 1) for x, y in zip(range(5), range(5))]
         df = GeoDataFrame(data, crs=self.crs)
         assert isinstance(df.geometry, GeoSeries)
diff --git a/geopandas/tests/test_testing.py b/geopandas/tests/test_testing.py
index 173b0d39..143faacc 100644
--- a/geopandas/tests/test_testing.py
+++ b/geopandas/tests/test_testing.py
@@ -1,8 +1,10 @@
 import numpy as np
 
 from shapely.geometry import Point, Polygon
+from pandas import Series
 
 from geopandas import GeoDataFrame, GeoSeries
+from geopandas.array import from_shapely
 
 from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
 import pytest
@@ -20,12 +22,32 @@ s2 = GeoSeries(
     ]
 )
 
+
+s3 = Series(
+    [
+        Polygon([(0, 2), (0, 0), (2, 0), (2, 2)]),
+        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
+    ]
+)
+
+a = from_shapely(
+    [
+        Polygon([(0, 2), (0, 0), (2, 0), (2, 2)]),
+        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
+    ]
+)
+
+s4 = Series(a)
+
 df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1})
 df2 = GeoDataFrame({"col1": [1, 2], "geometry": s2})
 
 
 def test_geoseries():
     assert_geoseries_equal(s1, s2)
+    assert_geoseries_equal(s1, s3, check_series_type=False, check_dtype=False)
+    assert_geoseries_equal(s3, s2, check_series_type=False, check_dtype=False)
+    assert_geoseries_equal(s1, s4, check_series_type=False)
 
     with pytest.raises(AssertionError):
         assert_geoseries_equal(s1, s2, check_less_precise=True)