Repository: geopandas/geopandas

Issue:
ENH: Should round trip I/O of GPKG automatically use schema for non-numeric dtypes?
#### Question about geopandas

I'm surprised that writing to GPKG format, then reloading from disk, does not complete a round trip for all dtypes. This schema information is written to the file, but gpd.read_file() results in `object` dtypes rather than `string` and `datetime64[ns]`.

I noticed that gf.to_parquet('test.parquet') and gf = read_parquet('test.parquet') does retain `string` and `datetime64[ns]`, so that is a current workaround, but it would be convenient if GPKG did this as well for workflows blending geopandas and QGIS, etc (https://github.com/geopandas/geopandas/issues/1437#issuecomment-631636305).

```python
import geopandas as gpd
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'), ignore_fields=["iso_a3", "continent"], rows=2)
world['date'] = gpd.pd.to_datetime('2021-10-12')
world['name'] = world.name.astype('string')
print(world.dtypes)
world.to_file('world.gpkg', driver='GPKG')
```

```
pop_est                int64
name                  string
gdp_md_est           float64
geometry            geometry
date          datetime64[ns]
dtype: object
```

```python
round_trip = gpd.read_file('world.gpkg')
round_trip.dtypes
```
```
pop_est          int64
name            object
gdp_md_est     float64
date            object
geometry      geometry
dtype: object
```

```python
with fiona.open('world.gpkg') as f:
    print(f.schema)
```

```
{'properties': OrderedDict([('pop_est', 'int'), ('name', 'str'), ('gdp_md_est', 'float'), ('date', 'datetime')]), 'geometry': 'Unknown'}
```


<details>
  <summary>gpd.show_versions()</summary>
  
```
SYSTEM INFO
-----------
python     : 3.9.6 | packaged by conda-forge | (default, Jul 11 2021, 03:39:48)  [GCC 9.3.0]
executable : /home/scott/miniconda3/envs/isce2gimp/bin/python
machine    : Linux-5.4.0-80-generic-x86_64-with-glibc2.31

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.9.1
GEOS lib   : /home/scott/miniconda3/envs/isce2gimp/lib/libgeos_c.so
GDAL       : 3.2.2
GDAL data dir: /home/scott/miniconda3/envs/isce2gimp/lib/python3.9/site-packages/fiona/gdal_data
PROJ       : 8.0.1
PROJ data dir: /home/scott/miniconda3/envs/isce2gimp/share/proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.10.1
pandas     : 1.3.3
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      : None
psycopg2   : None
geoalchemy2: None
pyarrow    : 5.0.0
pygeos     : None
```

</details>



Discussion:
m-richards: There was still some issues with microseconds that now seem to be sorted in terms of the tests, the issue is actually on the read side of fiona, rather than write.
```geojson
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "b": "2021-11-21T01:07:43.018" }, "geometry": { "type": "Point", "coordinates": [ 0.0, 0.0 ] } }
]
}
```

```python
In [1]: import fiona
In [2]: with fiona.open('test_datetime.geojson') as src:
    ...:     for f in src:
    ...:         print(f)
    ...:
{'type': 'Feature', 'id': '0', 'properties': OrderedDict([('b', '2021-11-21T01:07:43.018002')]), 'geometry': {'type': 'Point', 'coordinates': (0.0, 0.0)}}
```


m-richards: I think that should be everything resolved now. There is one remaining slightly pathological edge case, as I applied the suggestion of dropping microseconds on read. If you had e.g a GeoJSON written by some other tool than fiona (I guess actually other than GDAL), you would have previously been able to read up to nanosecond precision, but after this change, that would be truncated to milliseconds. 
jorisvandenbossche: > you would have previously been able to read up to nanosecond precision, but after this change, that would be truncated to milliseconds.

That might actually also be "a good thing"? If this GeoJSON file is read through GDAL, I think the nanoseconds will be incorrectly read back anyway given the precision limitation of float32? 

Patch:
diff --git a/doc/source/docs/user_guide/interactive_mapping.ipynb b/doc/source/docs/user_guide/interactive_mapping.ipynb
index 65944edc..58fe74dc 100644
--- a/doc/source/docs/user_guide/interactive_mapping.ipynb
+++ b/doc/source/docs/user_guide/interactive_mapping.ipynb
@@ -2,6 +2,8 @@
  "cells": [
   {
    "cell_type": "markdown",
+   "id": "c554e753",
+   "metadata": {},
    "source": [
     "# Interactive mapping\n",
     "\n",
@@ -10,48 +12,66 @@
     "Creating maps for interactive exploration mirrors the API of [static plots](../reference/api/geopandas.GeoDataFrame.plot.html) in an [explore()](../reference/api/geopandas.GeoDataFrame.explore.html) method of a GeoSeries or GeoDataFrame.\n",
     "\n",
     "Loading some example data:"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "caf2fbd5",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "import geopandas\n",
     "\n",
     "nybb = geopandas.read_file(geopandas.datasets.get_path('nybb'))\n",
     "world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))\n",
     "cities = geopandas.read_file(geopandas.datasets.get_path('naturalearth_cities'))"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "56bf1bcf",
+   "metadata": {},
    "source": [
     "The simplest option is to use `GeoDataFrame.explore()`:"
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "6b484ecc",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "nybb.explore()"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "7a797389",
+   "metadata": {},
    "source": [
     "Interactive plotting offers largely the same customisation as static one plus some features on top of that. Check the code below which plots a customised choropleth map. You can use `\"BoroName\"` column with NY boroughs names as an input of the choropleth, show (only) its name in the tooltip on hover but show all values on click. You can also pass custom background tiles (either a name supported by folium, a name recognized by `xyzservices.providers.query_name()`, XYZ URL or `xyzservices.TileProvider` object), specify colormap (all supported by `matplotlib`) and specify black outline."
-   ],
-   "metadata": {}
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "798bf532",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-info\">\n",
+    "Note\n",
+    "\n",
+    "Note that the GeoDataFrame needs to have a CRS set if you want to use background tiles.\n",
+    "</div>"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "94b4ff24",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "nybb.explore( \n",
     "     column=\"BoroName\", # make choropleth based on \"BoroName\" column\n",
@@ -61,20 +81,22 @@
     "     cmap=\"Set1\", # use \"Set1\" matplotlib colormap\n",
     "     style_kwds=dict(color=\"black\") # use black outline\n",
     "    )"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "5a10291e",
+   "metadata": {},
    "source": [
     "The `explore()` method returns a `folium.Map` object, which can also be passed directly (as you do with `ax` in `plot()`). You can then use folium functionality directly on the resulting map. In the example below, you can plot two GeoDataFrames on the same map and add layer control using folium. You can also add additional tiles allowing you to change the background directly in the map."
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "cba9970b",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "import folium\n",
     "\n",
@@ -100,9 +122,7 @@
     "folium.LayerControl().add_to(m)  # use folium to add layer control\n",
     "\n",
     "m  # show map"
-   ],
-   "outputs": [],
-   "metadata": {}
+   ]
   }
  ],
  "metadata": {
@@ -126,4 +146,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
\ No newline at end of file
+}
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/io/file.py b/geopandas/io/file.py
index 59467818..a4c9e55c 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -236,14 +236,22 @@ def _read_file(filename, bbox=None, mask=None, rows=None, **kwargs):
                 f_filt = features
             # get list of columns
             columns = list(features.schema["properties"])
+            datetime_fields = [
+                k for (k, v) in features.schema["properties"].items() if v == "datetime"
+            ]
             if kwargs.get("ignore_geometry", False):
-                return pd.DataFrame(
+                df = pd.DataFrame(
                     [record["properties"] for record in f_filt], columns=columns
                 )
-
-            return GeoDataFrame.from_features(
-                f_filt, crs=crs, columns=columns + ["geometry"]
-            )
+            else:
+                df = GeoDataFrame.from_features(
+                    f_filt, crs=crs, columns=columns + ["geometry"]
+                )
+            for k in datetime_fields:
+                # fiona only supports up to ms precision, any microseconds are
+                # floating point rounding error
+                df[k] = pd.to_datetime(df[k]).dt.round(freq="ms")
+            return df
 
 
 def read_file(*args, **kwargs):
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..7c13d119 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -1,5 +1,6 @@
 from collections import OrderedDict
 import datetime
+from distutils.version import LooseVersion
 import io
 import os
 import pathlib
@@ -9,6 +10,8 @@ import numpy as np
 import pandas as pd
 
 import fiona
+import pytz
+from pandas.testing import assert_series_equal
 from shapely.geometry import Point, Polygon, box
 
 import geopandas
@@ -21,6 +24,9 @@ from geopandas.tests.util import PACKAGE_DIR, validate_boro_df
 import pytest
 
 
+FIONA_GE_1814 = str(fiona.__version__) >= LooseVersion("1.8.14")  # datetime roundtrip
+
+
 _CRS = "epsg:4326"
 
 
@@ -139,15 +145,48 @@ def test_to_file_bool(tmpdir, driver, ext):
     assert_correct_driver(tempfilename, ext)
 
 
-def test_to_file_datetime(tmpdir):
+TEST_DATE = datetime.datetime(2021, 11, 21, 1, 7, 43, 17500)
+eastern = pytz.timezone("US/Eastern")
+
+datetime_type_tests = (TEST_DATE, eastern.localize(TEST_DATE))
+
+
+@pytest.mark.parametrize(
+    "time", datetime_type_tests, ids=("naive_datetime", "datetime_with_timezone")
+)
+@pytest.mark.parametrize("driver,ext", driver_ext_pairs)
+def test_to_file_datetime(tmpdir, driver, ext, time):
     """Test writing a data file with the datetime column type"""
-    tempfilename = os.path.join(str(tmpdir), "test_datetime.gpkg")
+    if ext in (".shp", ""):
+        pytest.skip(f"Driver corresponding to ext {ext} doesn't support dt fields")
+    if time.tzinfo is not None and FIONA_GE_1814 is False:
+        # https://github.com/Toblerity/Fiona/pull/915
+        pytest.skip("Fiona >= 1.8.14 needed for timezone support")
+
+    tempfilename = os.path.join(str(tmpdir), f"test_datetime{ext}")
     point = Point(0, 0)
-    now = datetime.datetime.now()
-    df = GeoDataFrame({"a": [1, 2], "b": [now, now]}, geometry=[point, point], crs=4326)
-    df.to_file(tempfilename, driver="GPKG")
+
+    df = GeoDataFrame(
+        {"a": [1, 2], "b": [time, time]}, geometry=[point, point], crs=4326
+    )
+    if FIONA_GE_1814:
+        fiona_precision_limit = "ms"
+    else:
+        fiona_precision_limit = "s"
+    df["b"] = df["b"].dt.round(freq=fiona_precision_limit)
+
+    df.to_file(tempfilename, driver=driver)
     df_read = read_file(tempfilename)
-    assert_geoseries_equal(df.geometry, df_read.geometry)
+
+    assert_geodataframe_equal(df.drop(columns=["b"]), df_read.drop(columns=["b"]))
+    if df["b"].dt.tz is not None:
+        # US/Eastern becomes pytz.FixedOffset(-300) when read from file
+        # so compare fairly in terms of UTC
+        assert_series_equal(
+            df["b"].dt.tz_convert(pytz.utc), df_read["b"].dt.tz_convert(pytz.utc)
+        )
+    else:
+        assert_series_equal(df["b"], df_read["b"])
 
 
 @pytest.mark.parametrize("driver,ext", driver_ext_pairs)
@@ -214,7 +253,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