Repository: geopandas/geopandas

Issue:
BUG: Losing the geometry column should convert to a DataFrame
Branched out from the discussion in #1653 

- [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.

---

#### Code Sample, a copy-pastable example

From @jorisvandenbossche 's comment in #1653 

```python
df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_cities"))

type(df.drop("geometry", axis=1)) # geopandas.geodataframe.GeoDataFrame
type(df.reindex(columns=["name"])) # geopandas.geodataframe.GeoDataFrame
```

#### Problem description

When a GeoDataFrame loses its geometry column (via `.drop()` or `.reindex()`, for instance), the result should be a regular pandas DataFrame. This makes intuitive sense and would allow easier compatibility with the pandas API, as well as preventing issues such as the one found in #1653.

Right now, these methods are simply inherited from the DataFrame class but we can overwrite them in the GeoDataFrame to automate this downcasting when possible.

#### Expected Output

```python
type(df.drop("geometry", axis=1)) # pandas.core.frame.DataFrame
type(df.reindex(columns=["name"])) # pandas.core.frame.DataFrame
```



Discussion:
m-richards: Getting pretty close now. One issue I've noticed is stuff like this
```python
from geopandas import GeoDataFrame
from shapely.geometry import Point

N = 10

args = [{"geometry": Point(x, y), 'baz':'foo'
         } for x, y in zip(range(N), range(N))]
df = GeoDataFrame(args)

a = df.join(df, rsuffix='foo_')
b = df.join(df.drop(columns=['geometry']), rsuffix='foo_')
c = df.drop(columns=['geometry']).join(df, rsuffix='foo_')
d = df.drop(columns=['geometry']).join(df.drop(columns=['geometry']), rsuffix='foo_')

print(type(a))  # GeoDataFrame
print(type(b))  # GeoDataFrame
print(type(c))  # DataFrame (GeoDataFrame on master)
print(type(d))  # DataFrame
```
case c has changed from the old behaviour as dropping the geometry now makes the left object a dataframe instead. I'm not sure if that behaviour is desired or not? I'm also not sure how we could fix it if it's not desired, since the type of join is determined by the left, and hence the finalize will be from the left - and we have no control over that. 

I suspect there are going to be other issues like this that perhaps aren't tested - I only caught this one because we use join in `GeoDataFrame.explode` and we lost `attrs`. 

There's also an issue of the case of the empty GeoDataFrame. Currently on my branch 
`GeoDataFrame().copy()` returns an empty `DataFrame` - as there's no geometry -and so do any other cases that call the fallback constructor. Some of these cases (perhaps all) could be fixed using `__finalize__` but that does introduce a logical inconsistency if GeoDataFrames losing their set geometry become dataframes, but an empty GeoDataFrame does not. Perhaps the geometry is None case becomes special?

My preferred solution would be to alter the checks on empty files and empty input to `gdf.to_file` to include a geometry column e.g.:
```python
def test_to_file_empty(tmpdir):
    input_empty_df = GeoDataFrame() # ->  input_empty_df = GeoDataFrame(columns=['geometry'])
    tempfilename = os.path.join(str(tmpdir), "test.shp")
    with pytest.raises(ValueError, match="Cannot write empty DataFrame to file."):
        input_empty_df.to_file(tempfilename)
```
since a geodataframe with a geometry column containing no rows won't be cast. As far as I can see, situtations where this arises practically - after applying some kind of filter to a GeoDataFrame where there are no matches, that GeoDataFrame should still have columns set. I could be missing something there though.






m-richards: I've spotted what seems to be a rather difficult to solve issue with this approach - and that is inplace operations. 

I was having trouble seeing a way to get GeoDataFrames losing their geometry column to become DataFrames, without breaking `gdf.rename` and `df.copy` amongst other things (which is why I've left those tests as the only failing in the last commit). 

In trying to fix that, I've noticed this:
```python
df = GeoDataFrame({
    "geometry": from_shapely([Point(1, 1), Point(2, 2)]),
    "value": [1, 2]}
)
df2 = df.copy()
df2.drop(columns=['geometry'], inplace=True)
print(type(df2))  # GeoDataFrame

drop = df.drop(columns=['geometry'])
print(type(drop))  # DataFrame
```
This is not something we can fix with `_constructor` because it actually works - in `drop`
https://github.com/pandas-dev/pandas/blob/b687cd4d9e520666a956a60849568a98dd00c672/pandas/core/generic.py#L3934-L3939:
```python
        for axis, labels in axes.items():
            if labels is not None:
                obj = obj._drop_axis(labels, axis, level=level, errors=errors)
        if inplace:
            self._update_inplace(obj)
```
`obj` is a DataFrame on my branch, but `self` is a GeoDataFrame and `self._update_inplace` returns a `GeoDataFrame` and contains the comment that it deliberately doesn't call finalize. So that seems to kill progress in this direction I think, we'd still need to overload methods to handle this. 


jorisvandenbossche: Thanks for the further updates!

> [join] case c has changed from the old behaviour as dropping the geometry now makes the left object a dataframe instead.  I'm not sure if that behaviour is desired or not? I'm also not sure how we could fix it if it's not desired, since the type of join is determined by the left, and hence the finalize will be from the left - and we have no control over that.

That's probably OK, I think. It's also not `join` that changed behaviour, right? It's the `drop` that changed, and if you would have called `join` from a pandas DataFrame (`df.join(gdf)`) on master it would also result in a DataFrame and not GeoDataFrame (no change here).

> There's also an issue of the case of the empty GeoDataFrame. 

Hmm, that's something to think about some more. 
Special casing an empty GeoDataFrame might make sense (so it can be preserved). On the other hand, it might actually also be fine to have those become DataFrame's. In the end, there is no geometry column, so most geospatial methods would raise an error anyway.  

I suppose that being consistent (and updating tests to include a geometry column in tests with empty frames) as you propose makes most sense. Especially this:

> As far as I can see, situtations where this arises practically - after applying some kind of filter to a GeoDataFrame where there are no matches, that GeoDataFrame should still have columns set. 

feels correct to me.

> I've spotted what seems to be a rather difficult to solve issue with this approach - and that is inplace operations.

I think we can ignore this case somewhat. It's a corner case anyway, and the user is essentially doing something inplace that doesn't make much sense for a GeoDataFrame to be inplace. I think it's fine that this gives an unusable GeoDataFrame as result, since the user has an easy solution to it (i.e. not using inplace).

As you mention, it would require to still overwrite all those methods anyway (instead of purely relying on the `_constructor` mechanism). Personally, I don't think that's worth it just for supporting the `inplace` keyword (which might get deprecated in pandas). I would say that if we want to support it, we should bring it up with pandas and give this as an argument that `_update_inplace` should actually call `__finalize__`.

m-richards: So I thought I'd give a quick update on this. I've been making slow progress, but think this is fairly close to ticking all the right boxes now (but I've said that before...). I'm currently testing things on yet another branch (https://github.com/m-richards/geopandas/pull/4) but I'll keep this PR going, and will cherry-pick things or rebase & force into this branch once things are in a clean state. 


Given that tthis PR deals with things discussed in 4 issues
>  #1911, (and also the referenced #1856, #1153 and #544).
I think it's worth outlining exactly what the intended behaviour is - and hopefully this is also captured by the tests.

**GeoDataFrame conversion when geometry columns are lost**
The basic premise is that indexing operations e.g `iloc`, `loc`, `drop`, `reindex` should behave the same as `__getitem__` - in the situations that make sense:
https://github.com/geopandas/geopandas/blob/4e2c8343e58ff78642e8ecedee6f70840bda7a9e/geopandas/geodataframe.py#L1315-L1330
Specifically
- If indexing retains the set geometry column, the type `GeoDataFrame` is preserved. 
- If indexing loses the geometry column, the type is downcast to a `pd.DataFrame`. This downcasting should occur regardless of whether there are any GeoSeries in the result gdf whicha re not the geometry column - the presence of the geometry column is the important part, not the presence of geometry data. (Regardless of the type here, `GeoDataFrame.set_geometry` or the monkey patched `DataFrame.set_geometry` is needed to get a "working" `GeoDataFrame` and the decision to downcast mirrors the `__getitem__` implementation. 
- A concrete example:
```python
df = gpd.GeoDataFrame({
    "geometry": from_shapely([Point(1, 1), Point(2, 2), Point(3, 3)]),
    "value": [1, 2, 3]}
)

df['geometry2'] = df['geometry']
# subset with set geom col
print(type(df[['geometry', 'value']]))  # GeoDataFrame
print(type(df.loc[:, ['geometry', 'value']]))  # GeoDataFrame
# subset with unset geom col
print(type(df[['geometry2', 'value']]))  # pd.DataFrame
print(type(df.loc[:, ['geometry2', 'value']]))  # GeoDataFrame -> should be pd.DataFrame
# subset with no geom col
print(type(df[['value']]))  # pd.DataFrame
print(type(df.loc[:, ['value']]))  # GeoDataFrame -> should be pd.DataFrame
```
This is the core behaviour change, and originally I'd thought that reconciling `squeeze()` and `to_frame` probably belonged in the same category. I now think that should be a separate PR.

This additionally will change the return types of other pandas methods -which make use of `_constructor`, and I had dabbled with trying to ensure they were all consistent here. I now think that's also probably out of scope, if there are edge cases they can


m-richards: The `to_frame` stuff is not quite as neatly removable as I thought. 

For `test_geoseries.test_expanddim` to retain it's behaviour of promoting `GeoSeries.reset_index()` or `to_frame()` to a `GeoDataFrame` under the above changes we need to do something about the fact that if `GeoSeries.name is None`, the column name for the introduce geometry will be `0` from `reset_index` and so this will get immediately downcast by `GeoDataFrame._constructor`. So it probably makes sense to split this off into a separate PR for fixing `GeoSeries._expanddim` as in #2063  (probably why I made that issue but I'd forgotten).
m-richards: So (m-richards#4) is now closer to green, at least on windows...  - again more as a proof of concept than that is the branch that should be merged (well pending comments to the above I guess). 

The thing I've missed (which is what I get for trying to develop on windows) is there's problems with spatial indices, and tests are skiped on windows. `sindex` getting copied in situations where it should be the same. I'm not super familiar with the mechanisms of the spatial indices. From a quick glance they're tied to the the underlying `GeometryArray`, which would suggest that the new `_constructor` is needlessly making copies, but I haven't dug into that fully.
jorisvandenbossche: @m-richards thanks for all this digging! For the practical way to move forward here (given this PR and the additional changes in your new branch): in the end I am fine to follow what you think is best / easiest. But to make a suggestion:

- I would start with a PR (eg this one) with just a basic `_geodataframe_constructor_with_fallback` implementation (checking that it has a geometry column, and if not converting to DataFrame). That will not fix everything, but I think should already help in some cases, right? 
- We can maybe keep all the tests you already added in this PR, but add some skips / comment out some parts that are not yet working (I can start reviewing the tests you added here in this PR)
- And then for the additional changes you have been looking into at https://github.com/m-richards/geopandas/pull/4, you can do follow-up PRs focusing on certain aspects / tests.

(it might not be that easy to keep things "separated", though, so we have to see a bit what is practical)

m-richards: 
> * I would start with a PR (eg this one) with just a basic `_geodataframe_constructor_with_fallback` implementation (checking that it has a geometry column, and if not converting to DataFrame). That will not fix everything, but I think should already help in some cases, right?

If you're happy with this, I think this is definitely the way forward, the no geometry case is simpler, will help in some cases, and doesn't have as many cross dependencies. Guess I was trying to avoid an "incomplete" solution, but it will be too much to manage otherwise.
jorisvandenbossche: Sounds good. So we can basically take the current state of this PR as a starting point for that? (I can review it a some more detail then)
m-richards: Ah, perhaps I should do a quick check first to avoid wasting your time, I need to remember where this branch is at again.
m-richards: @jorisvandenbossche This would be better to look at now, some extra commented out bits, and might be a few changes which over-reach, I can't remember if they're all exactly required for the changes that are in this PR now. 

This is the simpler bit, once we try to feed though where the active geometry column is things get messier
m-richards: Have been looking into version compatibility with pandas, after this push, the only failing tests should be convert_dtypes test I added for pd10 environment.

The concatenation of series is hard coded to return a DataFrame:
https://github.com/pandas-dev/pandas/blob/b687cd4d9e520666a956a60849568a98dd00c672/pandas/core/reshape/concat.py#L470-L477.

Given this, can we just mark this test as and expected failure on pandas 1.0? We could overload concat and essentially backport the method from pandas 1.1, but this seems a lot to fix an edge case for a nice to have but non critical behaviour thats only in one version of pandas (we don't have this problem in 37-minimal on pandas 0.25)
jorisvandenbossche: > Given this, can we just mark this test as and expected failure on pandas 1.0?

Yep, that's fine (I don't think there is a need to do a complicated backport from pandas)


m-richards: Thanks @jorisvandenbossche for having at look at this again and giving it some TLC. That's a much better solution to the line lengths.

I think the extra todos `TODO: should this set the geometry column to geometry2 or return a DataFrame?` are nice to have as that's still something I'm undecided on what the right answer is - but that's something for another PR.


jorisvandenbossche: @m-richards and many thanks to you for the bulk of the PR! 
I think it should be ready now to merge (I would like to include it in the release, as it fixes some actual bugs)
jorisvandenbossche: So some follow-ups:

- Decide what to do with results that loose the "geometry" column, but still have another geometry dtype column (currently this is a bit inconsistent, `__getitem__` gives DataFrame, but most other methods result in an "invalid" GeoDataFrame with `_geometry_column_name` not set)
- Decide what to do with "empty" GeoDataFrames (`GeoDataFrame()`), and that thus have no geometry column name, those will be converted into a DataFrame in an operation (eg `copy()`)
- Update `GeoSeries._constructor_expanddim` so that `to_frame()` works correctly?

Do you remember others? You also still have your other branches that did more extensive changes that are worth rebasing once this is merged?
jorisvandenbossche: @martinfleis no need to go through the detailed code discussion above (and details about interaction with pandas), but maybe a quick sanity check that you are OK with the general behavioural changes: when a method (like drop, loc/iloc, reindex, ..) drops the geometry column, we will now get back a DataFrame instead of a (kind-of-invalid) GeoDataFrame.
m-richards: Have had a look through my notes and the branches, here's what I've come up with:
* Special case of the first point, what should `df[['geometry2']]` return, if behaviour consistent with `__getitem__` is returned (specifically thinking of `.squeeze()` and making this consistent/ useful (although I am yet to find a situation where I've found this to be a problem)
* Implicit in what you've said in referencing the other branches - and in some other PRs, make `_geodataframe_constructor_with_fallback` smarter, to be aware of active geometry column. This will be needed to mimic `__getitem__` - I've already done some testing of this. What is probably more important however is if we can track the active geometry column in this way, it might be possible to clean up some of the special case handling of this in `apply`, `__finalize__`, `astype` and potentially other methods.

Points which are related to the resolution of inconsistency with "`__getitem__` behaviour" (which are probably secondary to what you are asking here, but I still think they're worth noting)
* If the "`__getitem__` behaviour" is opted for, does this have implications for `.rename`, and does that make `df.rename(columns={"geometry":"geom"})` returning a geodataframe with an invalid geometry inconsistent with the rest of the api? I'm not sure if there is a historical reason why `rename` doesn't also do `rename_geometry` under the covers, but that might be the logical direction if this is percieved as inconsistent (I remember this tripping me up when I was new to geopandas).
* A similar but less related tangent, what should happen on direct assignment to a geometry column it becomes invalid, e.g. `df['geometry'] =3`. I guess both of these could also be rephrased as an explicit question about api, "When should a GeoDataFrame return a DataFrame?" [I don't think either of these cases should return dataframes, but I think if return types are conditional, it would be nice if the logic for this was laid out specifically, rather than something one discovers by chance/ in confusion.
jorisvandenbossche: @m-richards I opened a new issue with your last comment: https://github.com/geopandas/geopandas/issues/2133 (will add my list to it as well)
jorisvandenbossche: On second thought: we should maybe keep this for after the 0.10 release since there are still some follow-up issues (that might still change behaviour in related cases)? If we merge this directly after 0.10 we have plenty of time to ensure we covered all cases for the next release
m-richards: > On second thought: we should maybe keep this for after the 0.10 release since there are still some follow-up issues (that might still change behaviour in related cases)? If we merge this directly after 0.10 we have plenty of time to ensure we covered all cases for the next release

That makes sense to me, as much as I would like to see this in a release, it would be good to have a complete solution catching all cases, and those changes should probably happen contained within a single release, rather than spread across them.
m-richards: Hey @jorisvandenbossche, is this mergable now that 0.10 is out? (sorry if this is an unecessary ping)
jorisvandenbossche: Yes, this is certainly mergeable. The main reason I was waiting a bit longer was to make it easier to do bug-fix releases from the main master branch (and once we merge this, we have to work with a separate branch and backport fixes to that, which is no problem to be clear, just a bit more work). 

Patch:
diff --git a/geopandas/array.py b/geopandas/array.py
index fd7b7fc6..bf3106a9 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -422,13 +422,16 @@ class GeometryArray(ExtensionArray):
             return self.__dict__
 
     def __setstate__(self, state):
-        if compat.USE_PYGEOS:
-            geoms = pygeos.from_wkb(state[0])
+        if not isinstance(state, dict):
+            # pickle file saved with pygeos
+            geoms = vectorized.from_wkb(state[0])
             self._crs = state[1]
             self._sindex = None  # pygeos.STRtree could not be pickled yet
             self.data = geoms
             self.base = None
         else:
+            if compat.USE_PYGEOS:
+                state["data"] = vectorized.from_shapely(state["data"])
             if "_crs" not in state:
                 state["_crs"] = None
             self.__dict__.update(state)
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 787a8932..e29b2213 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -13,7 +13,7 @@ from pyproj import CRS
 
 from geopandas.array import GeometryArray, GeometryDtype, from_shapely, to_wkb, to_wkt
 from geopandas.base import GeoPandasBase, is_geometry_type
-from geopandas.geoseries import GeoSeries
+from geopandas.geoseries import GeoSeries, _geoseries_constructor_with_fallback
 import geopandas.io
 from geopandas.explore import _explore
 from . import _compat as compat
@@ -23,6 +23,20 @@ from ._decorator import doc
 DEFAULT_GEO_COLUMN_NAME = "geometry"
 
 
+def _geodataframe_constructor_with_fallback(*args, **kwargs):
+    """
+    A flexible constructor for GeoDataFrame._constructor, which falls back
+    to returning a DataFrame (if a certain operation does not preserve the
+    geometry column)
+    """
+    df = GeoDataFrame(*args, **kwargs)
+    geometry_cols_mask = df.dtypes == "geometry"
+    if len(geometry_cols_mask) == 0 or geometry_cols_mask.sum() == 0:
+        df = pd.DataFrame(df)
+
+    return df
+
+
 def _ensure_geometry(data, crs=None):
     """
     Ensure the data is of geometry dtype or converted to it.
@@ -276,6 +290,9 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
             frame = self
         else:
             frame = self.copy()
+            # if there is no previous self.geometry, self.copy() will downcast
+            if type(frame) == DataFrame:
+                frame = GeoDataFrame(frame)
 
         to_remove = None
         geo_column_name = self._geometry_column_name
@@ -1390,20 +1407,26 @@ individually so that features may have different properties
         result = super().apply(
             func, axis=axis, raw=raw, result_type=result_type, args=args, **kwargs
         )
-        if (
-            isinstance(result, GeoDataFrame)
-            and self._geometry_column_name in result.columns
-            and isinstance(result[self._geometry_column_name].dtype, GeometryDtype)
-        ):
-            # apply calls _constructor which resets geom col name to geometry
-            result._geometry_column_name = self._geometry_column_name
-            if self.crs is not None and result.crs is None:
-                result.set_crs(self.crs, inplace=True)
+        # Reconstruct gdf if it was lost by apply
+        if self._geometry_column_name in result.columns:
+            # axis=1 apply will split GeometryDType to object, try and cast back
+            try:
+                result = result.set_geometry(self._geometry_column_name)
+            except TypeError:
+                pass
+            else:
+                if self.crs is not None and result.crs is None:
+                    result.set_crs(self.crs, inplace=True)
+
         return result
 
     @property
     def _constructor(self):
-        return GeoDataFrame
+        return _geodataframe_constructor_with_fallback
+
+    @property
+    def _constructor_sliced(self):
+        return _geoseries_constructor_with_fallback
 
     def __finalize__(self, other, method=None, **kwargs):
         """propagate metadata from other to self"""
@@ -1671,6 +1694,7 @@ individually so that features may have different properties
         df = (
             df_copy.drop(df_copy._geometry_column_name, axis=1)
             .join(exploded_geom)
+            .set_geometry(self._geometry_column_name)
             .__finalize__(self)
         )
 
@@ -1688,8 +1712,7 @@ individually so that features may have different properties
         if f"__{level_str}" in df.columns:
             df = df.rename(columns={f"__{level_str}": level_str})
 
-        geo_df = df.set_geometry(self._geometry_column_name)
-        return geo_df
+        return df
 
     # overrides the pandas astype method to ensure the correct return type
     def astype(self, dtype, copy=True, errors="raise", **kwargs):
@@ -2209,3 +2232,6 @@ def _dataframe_set_geometry(self, col, drop=False, inplace=False, crs=None):
 
 
 DataFrame.set_geometry = _dataframe_set_geometry
+
+if compat.PANDAS_GE_10 and not compat.PANDAS_GE_11:  # i.e. on pandas 1.0.x
+    _geodataframe_constructor_with_fallback._from_axes = GeoDataFrame._from_axes
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index d60d1fea..1f3bde69 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -627,13 +627,6 @@ class GeoSeries(GeoPandasBase, Series):
                 result.set_crs(self.crs, inplace=True)
         return result
 
-    def __finalize__(self, other, method=None, **kwargs):
-        """propagate metadata from other to self"""
-        # NOTE: backported from pandas master (upcoming v0.13)
-        for name in self._metadata:
-            object.__setattr__(self, name, getattr(other, name, None))
-        return self
-
     def isna(self):
         """
         Detect missing values.
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 2241925e..c693ca14 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -343,6 +343,7 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
             )
             .set_index(index_right)
             .drop(["_key_left", "_key_right"], axis=1)
+            .set_geometry(right_df.geometry.name)
         )
         if isinstance(index_right, list):
             joined.index.names = right_index_name


Test patch:
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index 89430d40..f3712824 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -214,7 +214,7 @@ def test_to_file_int64(tmpdir, df_points):
 
 
 def test_to_file_empty(tmpdir):
-    input_empty_df = GeoDataFrame()
+    input_empty_df = GeoDataFrame(columns=["geometry"])
     tempfilename = os.path.join(str(tmpdir), "test.shp")
     with pytest.raises(ValueError, match="Cannot write empty DataFrame to file."):
         input_empty_df.to_file(tempfilename)
diff --git a/geopandas/io/tests/test_pickle.py b/geopandas/io/tests/test_pickle.py
index 0327714d..bd27adad 100644
--- a/geopandas/io/tests/test_pickle.py
+++ b/geopandas/io/tests/test_pickle.py
@@ -2,6 +2,7 @@
 See generate_legacy_storage_files.py for the creation of the legacy files.
 
 """
+from contextlib import contextmanager
 from distutils.version import LooseVersion
 import glob
 import os
@@ -36,12 +37,14 @@ def legacy_pickle(request):
     return request.param
 
 
-@pytest.fixture
-def with_use_pygeos_false():
+@contextmanager
+def with_use_pygeos(option):
     orig = geopandas.options.use_pygeos
-    geopandas.options.use_pygeos = not orig
-    yield
-    geopandas.options.use_pygeos = orig
+    geopandas.options.use_pygeos = option
+    try:
+        yield
+    finally:
+        geopandas.options.use_pygeos = orig
 
 
 @pytest.mark.skipif(
@@ -70,13 +73,41 @@ def test_round_trip_current(tmpdir, current_pickle_data):
         assert isinstance(result.has_sindex, bool)
 
 
-@pytest.mark.skipif(not compat.HAS_PYGEOS, reason="requires pygeos to test #1745")
-def test_pygeos_switch(tmpdir, with_use_pygeos_false):
-    gdf_crs = geopandas.GeoDataFrame(
+def _create_gdf():
+    return geopandas.GeoDataFrame(
         {"a": [0.1, 0.2, 0.3], "geometry": [Point(1, 1), Point(2, 2), Point(3, 3)]},
         crs="EPSG:4326",
     )
-    path = str(tmpdir / "gdf_crs.pickle")
-    gdf_crs.to_pickle(path)
-    result = pd.read_pickle(path)
-    assert_geodataframe_equal(result, gdf_crs)
+
+
+@pytest.mark.skipif(not compat.HAS_PYGEOS, reason="requires pygeos to test #1745")
+def test_pygeos_switch(tmpdir):
+    # writing and reading with pygeos disabled
+    with with_use_pygeos(False):
+        gdf = _create_gdf()
+        path = str(tmpdir / "gdf_crs1.pickle")
+        gdf.to_pickle(path)
+        result = pd.read_pickle(path)
+        assert_geodataframe_equal(result, gdf)
+
+    # writing without pygeos, reading with pygeos
+    with with_use_pygeos(False):
+        gdf = _create_gdf()
+        path = str(tmpdir / "gdf_crs1.pickle")
+        gdf.to_pickle(path)
+
+    with with_use_pygeos(True):
+        result = pd.read_pickle(path)
+        gdf = _create_gdf()
+        assert_geodataframe_equal(result, gdf)
+
+    # writing with pygeos, reading without pygeos
+    with with_use_pygeos(True):
+        gdf = _create_gdf()
+        path = str(tmpdir / "gdf_crs1.pickle")
+        gdf.to_pickle(path)
+
+    with with_use_pygeos(False):
+        result = pd.read_pickle(path)
+        gdf = _create_gdf()
+        assert_geodataframe_equal(result, gdf)
diff --git a/geopandas/tests/test_op_output_types.py b/geopandas/tests/test_op_output_types.py
new file mode 100644
index 00000000..3be6b357
--- /dev/null
+++ b/geopandas/tests/test_op_output_types.py
@@ -0,0 +1,216 @@
+import pandas as pd
+import pyproj
+import pytest
+
+from shapely.geometry import Point
+
+from geopandas import GeoDataFrame, GeoSeries
+
+
+crs_osgb = pyproj.CRS(27700)
+crs_wgs = pyproj.CRS(4326)
+
+
+N = 10
+
+
+@pytest.fixture(params=["geometry", "point"])
+def df(request):
+    geo_name = request.param
+
+    df = GeoDataFrame(
+        [
+            {
+                "value1": x + y,
+                "value2": x * y,
+                geo_name: Point(x, y),  # rename this col in tests
+            }
+            for x, y in zip(range(N), range(N))
+        ],
+        crs=crs_wgs,
+        geometry=geo_name,
+    )
+    # want geometry2 to be a GeoSeries not Series, test behaviour of non geom col
+    df["geometry2"] = df[geo_name].set_crs(crs_osgb, allow_override=True)
+    return df
+
+
+def _check_metadata_gdf(gdf, geo_name="geometry", crs=crs_wgs):
+    assert gdf._geometry_column_name == geo_name
+    assert gdf.geometry.name == geo_name
+    assert gdf.crs == crs
+
+
+def _check_metadata_gs(gs, name="geometry", crs=crs_wgs):
+    assert gs.name == name
+    assert gs.crs == crs
+
+
+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
+    and crs match. If geo_name is None, then we expect a GeoDataFrame
+    where the geometry column is invalid/ isn't set. This is never desirable,
+    but is a reality of this first stage of implementation.
+    """
+    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:
+            with pytest.raises(AttributeError, match="No geometry data set yet"):
+                result.geometry.name  # be explicit that geometry is invalid here
+    elif expected_type == GeoSeries:
+        _check_metadata_gs(result, name=geo_name, crs=crs)
+
+
+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"]], pd.DataFrame)
+    assert_object(df[["geometry2"]], pd.DataFrame)
+    assert_object(df[["value1"]], pd.DataFrame)
+    # Series
+    assert_object(df[geo_name], GeoSeries, geo_name)
+    assert_object(df["geometry2"], GeoSeries, "geometry2", crs=crs_osgb)
+    assert_object(df["value1"], pd.Series)
+
+
+def test_loc(df):
+    geo_name = df.geometry.name
+    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)
+    # TODO: should this set the geometry column to geometry2 or return a DataFrame?
+    assert_object(df.loc[:, ["geometry2", "value1"]], GeoDataFrame, None)
+    assert_object(df.loc[:, ["geometry2"]], GeoDataFrame, None)
+    # Ideally this would mirror __getitem__ and below would be true
+    # assert_object(df.loc[:, ["geometry2", "value1"]], pd.DataFrame)
+    # assert_object(df.loc[:, ["geometry2"]], pd.DataFrame)
+    assert_object(df.loc[:, ["value1"]], pd.DataFrame)
+    # Series
+    assert_object(df.loc[:, geo_name], GeoSeries, geo_name)
+    assert_object(df.loc[:, "geometry2"], GeoSeries, "geometry2", crs=crs_osgb)
+    assert_object(df.loc[:, "value1"], pd.Series)
+
+
+def test_iloc(df):
+    geo_name = df.geometry.name
+    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)
+    # TODO: should this set the geometry column to geometry2 or return a DataFrame?
+    assert_object(df.iloc[:, [3, 0]], GeoDataFrame, None)
+    assert_object(df.iloc[:, [3]], GeoDataFrame, None)
+    # Ideally this would mirror __getitem__ and below would be true
+    # assert_object(df.iloc[:, [3, 0]], pd.DataFrame)
+    # assert_object(df.iloc[:, [3]], pd.DataFrame)
+    assert_object(df.iloc[:, [0]], pd.DataFrame)
+    # Series
+    assert_object(df.iloc[:, 2], GeoSeries, geo_name)
+    assert_object(df.iloc[:, 3], GeoSeries, "geometry2", crs=crs_osgb)
+    assert_object(df.iloc[:, 0], pd.Series)
+
+
+def test_squeeze(df):
+    geo_name = df.geometry.name
+    assert_object(df[[geo_name]].squeeze(), GeoSeries, geo_name)
+
+    # Not ideal behaviour, but this is consistent with __getitem__
+    assert_object(df[["geometry2"]].squeeze(), pd.Series)
+
+
+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
+
+    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
+
+    res3 = df["value1"].to_frame()
+    assert_object(res3, pd.DataFrame)
+
+
+def test_reindex(df):
+    geo_name = df.geometry.name
+    assert_object(df.reindex(columns=["value1", "value2"]), pd.DataFrame)
+    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)
+    # TODO: should this set the geometry column to geometry2 or return a DataFrame?
+    assert_object(df.reindex(columns=["geometry2", "value1"]), GeoDataFrame, None)
+    assert_object(df.reindex(columns=["geometry2"]), GeoDataFrame, None)
+    # Ideally this would mirror __getitem__ and below would be true
+    # assert_object(df.reindex(columns=["geometry2", "value1"]), pd.DataFrame)
+    # assert_object(df.reindex(columns=["geometry2"]), pd.DataFrame)
+    assert_object(df.reindex(columns=["value1"]), pd.DataFrame)
+
+    # reindexing the rows always preserves the GeoDataFrame
+    assert_object(df.reindex(index=[0, 1, 20]), GeoDataFrame, geo_name)
+
+    # reindexing both rows and columns
+    assert_object(
+        df.reindex(index=[0, 1, 20], columns=[geo_name]), GeoDataFrame, geo_name
+    )
+    assert_object(df.reindex(index=[0, 1, 20], columns=["value1"]), pd.DataFrame)
+
+
+def test_drop(df):
+    geo_name = df.geometry.name
+    assert_object(df.drop(columns=[geo_name, "geometry2"]), pd.DataFrame)
+    assert_object(df.drop(columns=["value1", "value2"]), GeoDataFrame, geo_name)
+    cols = ["value1", "value2", "geometry2"]
+    assert_object(df.drop(columns=cols), GeoDataFrame, geo_name)
+    # TODO: should this set the geometry column to geometry2 or return a DataFrame?
+    assert_object(df.drop(columns=[geo_name, "value2"]), GeoDataFrame, None)
+    assert_object(df.drop(columns=["value1", "value2", geo_name]), GeoDataFrame, None)
+    # Ideally this would mirror __getitem__ and below would be true
+    # assert_object(df.drop(columns=[geo_name, "value2"]), pd.DataFrame)
+    # assert_object(df.drop(columns=["value1", "value2", geo_name]), pd.DataFrame)
+    assert_object(df.drop(columns=["geometry2", "value2", geo_name]), pd.DataFrame)
+
+
+def test_apply(df):
+    geo_name = df.geometry.name
+
+    def identity(x):
+        return x
+
+    # axis = 0
+    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), pd.DataFrame)
+    assert_object(df[["geometry2"]].apply(identity), pd.DataFrame)
+    assert_object(df[["value1"]].apply(identity), pd.DataFrame)
+
+    # axis = 0, Series
+    assert_object(df[geo_name].apply(identity), GeoSeries, geo_name)
+    assert_object(df["geometry2"].apply(identity), GeoSeries, "geometry2", crs=crs_osgb)
+    assert_object(df["value1"].apply(identity), pd.Series)
+
+    # axis = 0, Series, no longer geometry
+    assert_object(df[geo_name].apply(lambda x: str(x)), pd.Series)
+    assert_object(df["geometry2"].apply(lambda x: str(x)), pd.Series)
+
+    # axis = 1
+    assert_object(df[["value1", "value2"]].apply(identity, axis=1), pd.DataFrame)
+    assert_object(
+        df[[geo_name, "geometry2"]].apply(identity, axis=1), GeoDataFrame, geo_name
+    )
+    assert_object(df[[geo_name]].apply(identity, axis=1), GeoDataFrame, geo_name)
+    assert_object(df[["geometry2", "value1"]].apply(identity, axis=1), pd.DataFrame)
+    assert_object(df[["geometry2"]].apply(identity, axis=1), pd.DataFrame)
+    assert_object(df[["value1"]].apply(identity, axis=1), pd.DataFrame)
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 9b553d4c..709ad3da 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -139,8 +139,9 @@ def test_reindex(s, df):
     assert isinstance(res.geometry, GeoSeries)
     assert_frame_equal(res, df[["value1", "geometry"]])
 
-    # TODO df.reindex(columns=['value1', 'value2']) still returns GeoDataFrame,
-    # should it return DataFrame instead ?
+    res = df.reindex(columns=["value1", "value2"])
+    assert type(res) == pd.DataFrame
+    assert_frame_equal(res, df[["value1", "value2"]])
 
 
 def test_take(s, df):
@@ -587,13 +588,12 @@ def test_apply_no_geometry_result(df, crs):
     if crs:
         df = df.set_crs(crs)
     result = df.apply(lambda col: col.astype(str), axis=0)
-    # TODO this should actually not return a GeoDataFrame
-    assert isinstance(result, GeoDataFrame)
+    assert type(result) is pd.DataFrame
     expected = df.astype(str)
     assert_frame_equal(result, expected)
 
     result = df.apply(lambda col: col.astype(str), axis=1)
-    assert isinstance(result, GeoDataFrame)
+    assert type(result) is pd.DataFrame
     assert_frame_equal(result, expected)
 
 
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index cce912f6..d1d88519 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -303,7 +303,7 @@ class TestPointPlotting:
         with pytest.warns(UserWarning):
             ax = s.plot()
         assert len(ax.collections) == 0
-        df = GeoDataFrame([])
+        df = GeoDataFrame([], columns=["geometry"])
         with pytest.warns(UserWarning):
             ax = df.plot()
         assert len(ax.collections) == 0