Repository: geopandas/geopandas
Pull request #2577 (2022-10-02 11:05:15)

ENH: set the default geometry column name to none
This is a proof of concept based upon discussion in https://github.com/geopandas/geopandas/issues/2517#issuecomment-1212852039 to change the default geometry column name to None, rather than assuming it is "geometry".

This will be failing some tests related to the current ability to set the crs of a GeoDataFrame without an active geometry column. This is meant to be deprecated in 0.12 though. Edit: opened #2578 to draft this

martinfleis: I like this change. Thanks for looking into that!
jorisvandenbossche: There are still a bunch of actual failures related to this change (eg https://github.com/geopandas/geopandas/actions/runs/3302728571/jobs/5449873350)

One of them is caused by this snippet:

```
df = GeoDataFrame()
df["geometry"] = scalar
df.crs = 4326
```

where the last line now fails to set the CRS because there is no active geometry column. And that's no longer the case because when creating the empty `df`, the default geometry column name is not "geometry", and so the second line does not set this active column.

Do we want to keep that snippet working? Or are we fine with this breaking?
m-richards: Ah sorry, I think I was aware of this but forgot when I wrote that comment.

>Do we want to keep that snippet working? Or are we fine with this breaking?

I personally think this isn't ideal to keep permitting, but I could also understand it's a breaking change and probably should have associated deprecation warnings. Seems tricky to deprecate though would have to intercept getitem specifically for "geometry" when the active geometry column is none (and perhaps that has some confusing edge cases associated with it?).
jorisvandenbossche: > Seems tricky to deprecate though would have to intercept getitem specifically for "geometry" when the active geometry column is none

That's maybe actually possible? We already override `__setitem__` on GeoDataFrame, so then inside we could do something like (untested pseudocode):

```
if self._geometry_column_name is None and isinstance(key, str) and key == "geometry" and "geometry" not in self.columns:
    warnings.warn(
        "you are adding a "geometry" column to a GeoDataFrame without active geometry column. "
        "Currently, this makes this the active geometry column, but in the future that will no longer "
        "happen automatically. Do ... instead"
    )
    self._geometry_column_name = "geometry"
...
m-richards: I've just had another go of reworking this, to try and deal with issuing a UserWarning in the right cases. To achieve this, the default geometry column name unfortunately is no longer None, but rather a sentinel object, so this is a bit uglier. I was a bit worried about issuing warnings in the wrong situations, rather than only on init - so locally I've been testing with `pytest -W error::FutureWarning` and making sure the new warning isn't issued in any unexpected places.

There's still an issue with multiindexes, which I'm hoping will be resolved by the other PR, and I've also just noticed that there's a handful of extra failures in some of the older CI envs. But I'll worrry about that after I get some input of this is still the right way forward.
jorisvandenbossche: > To achieve this, the default geometry column name unfortunately is no longer None, but rather a sentinel object, so this is a bit uglier.

No problem with using a custom sentinel (that's an approach we use in pandas as well if we want to have a default optional value where None cannot be used because it's a valid option). But just wondering: where there cases where `None` was a valid value for `_geometry_column_name` where we don't want to warn? Do you have an example of that?
jorisvandenbossche: > But just wondering: where there cases where `None` was a valid value for `_geometry_column_name` where we don't want to warn? Do you have an example of that?

Ah, if you drop the active geometry column (but there is still another geometry dtype column), this gets set to None. Is that a case that should be handled differently than when explicitly constructing a geodataframe without geometry column?

Patch:
diff --git a/.github/workflows/versioneer.yml b/.github/workflows/versioneer.yml
index a1814fbc..644056be 100644
--- a/.github/workflows/versioneer.yml
+++ b/.github/workflows/versioneer.yml
@@ -32,7 +32,7 @@ jobs:
           git reset -- geopandas/__init__.py
           git checkout -- geopandas/__init__.py
       - name: Create PR
-        uses: peter-evans/create-pull-request@v4
+        uses: peter-evans/create-pull-request@v5
         with:
           title: "Update Versioneer"
           branch: update-versioneer
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 90e028b3..93aef7b3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,8 @@ New features and improvements:
   using ``engine="pyogrio"`` (#2788).
 - Added a ``to_wgs84`` keyword to ``to_json`` allowing automatic re-projecting to follow
   the 2016 GeoJSON specification (#416).
+- Improve error messages when accessing the `geometry` attribute of GeoDataFrame without an active geometry column 
+  related to the default name `"geometry"` being provided in the constructor (#2577)
 
 Deprecations and compatibility notes:
 
@@ -40,6 +42,9 @@ Bug fixes:
   Pyogrio if provided a URL as input (#2796)
 - Fix `copy()` downcasting GeoDataFrames without an active geometry column to a
   DataFrame (#2775)
+- Fix geometry column name propagation when GeoDataFrame columns are a multiindex (#2088)
+- Fix `iterfeatures()` method of GeoDataFrame to correctly handle non-scalar values
+  when `na='drop'` is specified (#2811)
 
 Notes on (optional) dependencies:
 
diff --git a/doc/source/docs/advanced_guide.rst b/doc/source/docs/advanced_guide.rst
index 549d99a8..f6e54206 100644
--- a/doc/source/docs/advanced_guide.rst
+++ b/doc/source/docs/advanced_guide.rst
@@ -18,3 +18,4 @@ specification in the :doc:`API reference <reference>`.
 
   user_guide/missing_empty
   user_guide/reproject_fiona
+  user_guide/pygeos_to_shapely
diff --git a/doc/source/docs/user_guide/pygeos_to_shapely.rst b/doc/source/docs/user_guide/pygeos_to_shapely.rst
new file mode 100644
index 00000000..94887abb
--- /dev/null
+++ b/doc/source/docs/user_guide/pygeos_to_shapely.rst
@@ -0,0 +1,78 @@
+Migration from PyGEOS geometry backend to Shapely 2.0
+=====================================================
+
+Since the 0.8 version, GeoPandas includes an experimental support of PyGEOS as an
+alternative geometry backend to Shapely. Recently, PyGEOS codebase was merged into the
+Shapely project and released as part of Shapely 2.0. GeoPandas will therefore
+deprecate support of the PyGEOS backend and will go forward with Shapely 2.0 as the
+only geometry engine exposing GEOS functionality.
+
+Given that historically the PyGEOS engine is automatically used if the package is installed (this behaviour will change in GeoPandas 0.14 where Shapely 2.0 will be used by default if installed), some downstream code may depend on
+PyGEOS geometries being available as underlying data of a ``GeometryArray``.
+
+This guide outlines the migration from the PyGEOS-based code to the Shapely-based code.
+
+Migration period
+----------------
+
+The migration is planned for three releases spanning approximately one year, starting
+with 0.13 released in the second quarter of 2023.
+
+GeoPandas 0.13
+^^^^^^^^^^^^^^
+
+- PyGEOS is still used as a default backend over Shapely (1.8 or 2.0) if installed,
+  with a ``FutureWarning`` warning about upcoming changes.
+
+GeoPandas 0.14
+^^^^^^^^^^^^^^
+
+- The default backend will become Shapely 2.0 and the PyGEOS will be used only
+  if Shapely 1.8 is installed instead of 2.0 or newer. The PyGEOS backend will still be
+  supported, but a user will need to opt in using the environment variable
+  ``USE_PYGEOS`` as explained in the
+  `installation instructions <../getting_started/install>`__.
+
+GeoPandas 1.0
+^^^^^^^^^^^^^
+
+- GeoPandas will remove support of both PyGEOS and Shapely<2.
+
+How to prepare your code for transition
+---------------------------------------
+
+If you don't use PyGEOS explicitly, there nothing to be done as GeoPandas internals will
+take care of the transition. If you use PyGEOS directly and access an array of PyGEOS
+geometries using ``GeoSeries.values.data``, you will need to make some changes to avoid
+code breakage.
+
+The recommended way is using Shapely vectorized operations on the ``GeometryArray``
+instead of accessing the NumPy array of geometries and using PyGEOS/Shapely operations
+on the array.
+
+This is a common pattern used with GeoPandas 0.12 (or earlier), that should now be avoided in new code:
+
+.. code-block:: python
+
+    >>> import pygeos
+    >>> geometries = gdf.geometry.values.data
+    >>> mrr = pygeos.minimum_rotated_rectangle(geometries)
+
+The recommended way of refactoring this code would look like this (with Geopandas 0.12 or later):
+
+.. code-block:: python
+
+    >>> import shapely  # shapely 2.0
+    >>> mrr = shapely.minimum_rotated_rectangle(gdf.geometry.array)
+
+This code will work no matter which geometry backend GeoPandas actually uses, because on
+the ``GeometryArray`` level, it always returns Shapely geometry. Although keep in mind, that
+it may involve additional overhead cost of converting PyGEOS geometry to Shapely
+geometry.
+
+Note that while in most cases, a simple replacement of ``pygeos`` with ``shapely``
+together with a change of ``gdf.geometry.values.data`` to ``gdf.geometry.values`` or
+analogous ``gdf.geometry.array``  should work, there are some differences between the
+API of PyGEOS and that of Shapely. Please consult the
+`Migrating from PyGEOS <https://shapely.readthedocs.io/en/stable/migration_pygeos.html>`__
+document for details.
\ No newline at end of file
diff --git a/geopandas/datasets/__init__.py b/geopandas/datasets/__init__.py
index b56d202a..bed996ad 100644
--- a/geopandas/datasets/__init__.py
+++ b/geopandas/datasets/__init__.py
@@ -1,5 +1,7 @@
 import os
 
+from warnings import warn
+
 __all__ = ["available", "get_path"]
 
 _module_path = os.path.dirname(__file__)
@@ -25,9 +27,30 @@ def get_path(dataset):
 naturalearth_lowres/naturalearth_lowres.shp'
 
     """
+    ne_message = "https://www.naturalearthdata.com/downloads/110m-cultural-vectors/."
+    nybb_message = (
+        "the geodatasets package.\n\nfrom geodatasets import get_path\n"
+        "path_to_file = get_path('nybb')\n"
+    )
+    depr_warning = (
+        "The geopandas.dataset module is deprecated and will be removed in GeoPandas "
+        f"1.0. You can get the original '{dataset}' data from "
+        f"{ne_message if 'natural' in dataset else nybb_message}"
+    )
+
     if dataset in _available_dir:
+        warn(
+            depr_warning,
+            FutureWarning,
+            stacklevel=2,
+        )
         return os.path.abspath(os.path.join(_module_path, dataset, dataset + ".shp"))
     elif dataset in _available_zip:
+        warn(
+            depr_warning,
+            FutureWarning,
+            stacklevel=2,
+        )
         fpath = os.path.abspath(os.path.join(_module_path, _available_zip[dataset]))
         return "zip://" + fpath
     else:
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 0169dd49..0045eead 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -20,9 +20,6 @@ from . import _compat as compat
 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
@@ -128,7 +125,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
     _internal_names = DataFrame._internal_names + ["geometry"]
     _internal_names_set = set(_internal_names)
 
-    _geometry_column_name = DEFAULT_GEO_COLUMN_NAME
+    _geometry_column_name = None
 
     def __init__(self, data=None, *args, geometry=None, crs=None, **kwargs):
         with compat.ignore_shapely2_warnings():
@@ -241,6 +238,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
     def _set_geometry(self, col):
         if not pd.api.types.is_list_like(col):
             raise ValueError("Must use a list-like to set the geometry property")
+        self._persist_old_default_geometry_colname()
         self.set_geometry(col, inplace=True)
 
     geometry = property(
@@ -308,12 +306,11 @@ 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
+        if geo_column_name is None:
+            geo_column_name = "geometry"
         if isinstance(col, (Series, list, np.ndarray, GeometryArray)):
             level = col
         elif hasattr(col, "ndim") and col.ndim > 1:
@@ -333,7 +330,6 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
             if drop:
                 to_remove = col
-                geo_column_name = self._geometry_column_name
             else:
                 geo_column_name = col
 
@@ -350,8 +346,11 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
         # Check that we are using a listlike of geometries
         level = _ensure_geometry(level, crs=crs)
-        frame[geo_column_name] = level
+        # update _geometry_column_name prior to assignment
+        # to avoid default is None warning
         frame._geometry_column_name = geo_column_name
+        frame[geo_column_name] = level
+
         if not inplace:
             return frame
 
@@ -445,22 +444,22 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
     @crs.setter
     def crs(self, value):
         """Sets the value of the crs"""
-        if self._geometry_column_name not in self:
+        if self._geometry_column_name is None:
             raise ValueError(
                 "Assigning CRS to a GeoDataFrame without a geometry column is not "
                 "supported. Use GeoDataFrame.set_geometry to set the active "
                 "geometry column.",
             )
+
+        if hasattr(self.geometry.values, "crs"):
+            self.geometry.values.crs = value
         else:
-            if hasattr(self.geometry.values, "crs"):
-                self.geometry.values.crs = value
-            else:
-                # column called 'geometry' without geometry
-                raise ValueError(
-                    "Assigning CRS to a GeoDataFrame without an active geometry "
-                    "column is not supported. Use GeoDataFrame.set_geometry to set "
-                    "the active geometry column.",
-                )
+            # column called 'geometry' without geometry
+            raise ValueError(
+                "Assigning CRS to a GeoDataFrame without an active geometry "
+                "column is not supported. Use GeoDataFrame.set_geometry to set "
+                "the active geometry column.",
+            )
 
     def __setstate__(self, state):
         # overriding DataFrame method for compat with older pickles (CRS handling)
@@ -897,16 +896,22 @@ individually so that features may have different properties
 
         if len(properties_cols) > 0:
             # convert to object to get python scalars.
-            properties = self[properties_cols].astype(object).values
+            properties_cols = self[properties_cols]
+            properties = properties_cols.astype(object).values
+            na_mask = pd.isna(properties_cols).values
+
             if na == "null":
-                properties[pd.isnull(self[properties_cols]).values] = None
+                properties[na_mask] = None
 
             for i, row in enumerate(properties):
                 geom = geometries[i]
 
                 if na == "drop":
+                    na_mask_row = na_mask[i]
                     properties_items = {
-                        k: v for k, v in zip(properties_cols, row) if not pd.isnull(v)
+                        k: v
+                        for k, v, na in zip(properties_cols, row, na_mask_row)
+                        if not na
                     }
                 else:
                     properties_items = {k: v for k, v in zip(properties_cols, row)}
@@ -1444,6 +1449,18 @@ individually so that features may have different properties
         return a GeoDataFrame.
         """
         result = super().__getitem__(key)
+        # Custom logic to avoid waiting for pandas GH51895
+        # result is not geometry dtype for multi-indexes
+        if (
+            pd.api.types.is_scalar(key)
+            and key == ""
+            and isinstance(self.columns, pd.MultiIndex)
+            and isinstance(result, Series)
+            and not is_geometry_type(result)
+        ):
+            loc = self.columns.get_loc(key)
+            # squeeze stops multilevel columns from returning a gdf
+            result = self.iloc[:, loc].squeeze(axis="columns")
         geo_col = self._geometry_column_name
         if isinstance(result, Series) and isinstance(result.dtype, GeometryDtype):
             result.__class__ = GeoSeries
@@ -1456,17 +1473,43 @@ individually so that features may have different properties
                 result.__class__ = DataFrame
         return result
 
+    def _persist_old_default_geometry_colname(self):
+        """Internal util to temporarily persist the default geometry column
+        name of 'geometry' for backwards compatibility."""
+        # self.columns check required to avoid this warning in __init__
+        if self._geometry_column_name is None and "geometry" not in self.columns:
+            msg = (
+                "You are adding a column named 'geometry' to a GeoDataFrame "
+                "constructed without an active geometry column. Currently, "
+                "this automatically sets the active geometry column to 'geometry' "
+                "but in the future that will no longer happen. Instead, either "
+                "provide geometry to the GeoDataFrame constructor "
+                "(GeoDataFrame(... geometry=GeoSeries()) or use "
+                "`set_geometry('geometry')` "
+                "to explicitly set the active geometry column."
+            )
+            warnings.warn(msg, category=FutureWarning, stacklevel=3)
+            self._geometry_column_name = "geometry"
+
     def __setitem__(self, key, value):
         """
         Overwritten to preserve CRS of GeometryArray in cases like
         df['geometry'] = [geom... for geom in df.geometry]
         """
-        if not pd.api.types.is_list_like(key) and key == self._geometry_column_name:
+
+        if (
+            not pd.api.types.is_list_like(key)
+            and key == self._geometry_column_name
+            or key == "geometry"
+            and self._geometry_column_name is None
+        ):
             if pd.api.types.is_scalar(value) or isinstance(value, BaseGeometry):
                 value = [value] * self.shape[0]
             try:
                 crs = getattr(self, "crs", None)
                 value = _ensure_geometry(value, crs=crs)
+                if key == "geometry":
+                    self._persist_old_default_geometry_colname()
             except TypeError:
                 warnings.warn("Geometry column does not contain geometry.")
         super().__setitem__(key, value)
@@ -1515,6 +1558,13 @@ individually so that features may have different properties
         result = super().apply(
             func, axis=axis, raw=raw, result_type=result_type, args=args, **kwargs
         )
+        # pandas <1.4 re-attach last geometry col if lost
+        if (
+            not compat.PANDAS_GE_14
+            and isinstance(result, GeoDataFrame)
+            and result._geometry_column_name is None
+        ):
+            result._geometry_column_name = self._geometry_column_name
         # Reconstruct gdf if it was lost by apply
         if (
             isinstance(result, DataFrame)
diff --git a/pyproject.toml b/pyproject.toml
index 84be25f9..6cbcc162 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -77,5 +77,6 @@ xfail_strict = true
 
 filterwarnings = [
     "ignore:distutils Version classes are deprecated.*:DeprecationWarning:pandas.*",
-    "ignore:distutils Version classes are deprecated.*:DeprecationWarning:numpy.*"
+    "ignore:distutils Version classes are deprecated.*:DeprecationWarning:numpy.*",
+    "ignore:The geopandas.dataset module is deprecated"
 ]


Test patch:
diff --git a/geopandas/io/tests/test_arrow.py b/geopandas/io/tests/test_arrow.py
index 030ffbca..e999cca5 100644
--- a/geopandas/io/tests/test_arrow.py
+++ b/geopandas/io/tests/test_arrow.py
@@ -575,6 +575,15 @@ def test_missing_crs(tmpdir, file_format):
     assert_geodataframe_equal(df, pq_df, check_crs=True)
 
 
+def test_default_geo_col_writes(tmp_path):
+    # edge case geo col name None writes successfully
+    df = GeoDataFrame({"a": [1, 2]})
+    df.to_parquet(tmp_path / "test.pq")
+    # cannot be round tripped as gdf due to invalid geom col
+    pq_df = pd_read_parquet(tmp_path / "test.pq")
+    assert_frame_equal(df, pq_df)
+
+
 @pytest.mark.skipif(
     Version(pyarrow.__version__) >= Version("0.17.0"),
     reason="Feather only supported for pyarrow >= 0.17",
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index 20bd088b..5958cb6b 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -286,6 +286,7 @@ class TestGeometryArrayCRS:
         s = GeoSeries(arr, crs=27700)
         df = GeoDataFrame()
         df = df.set_geometry(s)
+        assert df._geometry_column_name == "geometry"
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
@@ -309,8 +310,16 @@ class TestGeometryArrayCRS:
         assert df.geometry.crs == self.wgs
         assert df.geometry.values.crs == self.wgs
 
-        # geometry column without geometry
+        # geometry column name None on init
         df = GeoDataFrame({"geometry": [0, 1]})
+        with pytest.raises(
+            ValueError,
+            match="Assigning CRS to a GeoDataFrame without a geometry",
+        ):
+            df.crs = 27700
+
+        # geometry column without geometry
+        df = GeoDataFrame({"geometry": [Point(0, 1)]}).assign(geometry=[0])
         with pytest.raises(
             ValueError,
             match="Assigning CRS to a GeoDataFrame without an active geometry",
@@ -337,14 +346,20 @@ class TestGeometryArrayCRS:
         arr = from_shapely(self.geoms)
         s = GeoSeries(arr, crs=27700)
         df = GeoDataFrame()
-        df["geometry"] = s
+        with pytest.warns(
+            FutureWarning, match="You are adding a column named 'geometry'"
+        ):
+            df["geometry"] = s
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
 
         arr = from_shapely(self.geoms, crs=27700)
         df = GeoDataFrame()
-        df["geometry"] = arr
+        with pytest.warns(
+            FutureWarning, match="You are adding a column named 'geometry'"
+        ):
+            df["geometry"] = arr
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
@@ -385,7 +400,10 @@ class TestGeometryArrayCRS:
     )
     def test_scalar(self, scalar):
         df = GeoDataFrame()
-        df["geometry"] = scalar
+        with pytest.warns(
+            FutureWarning, match="You are adding a column named 'geometry'"
+        ):
+            df["geometry"] = scalar
         df.crs = 4326
         assert df.crs == self.wgs
         assert df.geometry.crs == self.wgs
diff --git a/geopandas/tests/test_datasets.py b/geopandas/tests/test_datasets.py
index 0e4e02fe..2b4b14a9 100644
--- a/geopandas/tests/test_datasets.py
+++ b/geopandas/tests/test_datasets.py
@@ -8,4 +8,5 @@ import pytest
     "test_dataset", ["naturalearth_lowres", "naturalearth_cities", "nybb"]
 )
 def test_read_paths(test_dataset):
-    assert isinstance(read_file(get_path(test_dataset)), GeoDataFrame)
+    with pytest.warns(FutureWarning, match="The geopandas.dataset module is"):
+        assert isinstance(read_file(get_path(test_dataset)), GeoDataFrame)
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index b7ffaa48..6d6a6718 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -341,15 +341,16 @@ class TestDataFrame:
 
     def test_get_geometry_invalid(self):
         df = GeoDataFrame()
+        # no column "geometry" ever added
         df["geom"] = self.df.geometry
-        msg_geo_col_missing = "is not present. "
+        msg_geo_col_none = "active geometry column to use has not been set. "
 
-        with pytest.raises(AttributeError, match=msg_geo_col_missing):
+        with pytest.raises(AttributeError, match=msg_geo_col_none):
             df.geometry
-        df2 = self.df.copy()
-        df2["geom2"] = df2.geometry
-        df2 = df2[["BoroCode", "BoroName", "geom2"]]
-        with pytest.raises(AttributeError, match=msg_geo_col_missing):
+        # "geometry" originally present but dropped (but still a gdf)
+        col_subset_drop_geometry = ["BoroCode", "BoroName", "geom2"]
+        df2 = self.df.copy().assign(geom2=self.df.geometry)[col_subset_drop_geometry]
+        with pytest.raises(AttributeError, match="is not present."):
             df2.geometry
 
         msg_other_geo_cols_present = "There are columns with geometry data type"
@@ -884,6 +885,24 @@ class TestDataFrame:
         with pytest.raises(AttributeError):
             list(df.iterfeatures())
 
+    def test_geodataframe_iterfeatures_non_scalars(self):
+        # When some features in geodataframe are non-scalar values
+        df = GeoDataFrame(
+            {"geometry": [Point(1, 2)], "non-scalar": [[1, 2]], "test_col": None}
+        )
+        # null
+        expected = {"non-scalar": [1, 2], "test_col": None}
+        result = list(df.iterfeatures(na="null"))[0].get("properties")
+        assert expected == result
+        # drop
+        expected = {"non-scalar": [1, 2]}
+        result = list(df.iterfeatures(na="drop"))[0].get("properties")
+        assert expected == result
+        # keep
+        expected = {"non-scalar": [1, 2], "test_col": None}
+        result = list(df.iterfeatures(na="keep"))[0].get("properties")
+        assert expected == result
+
     def test_geodataframe_geojson_no_bbox(self):
         geo = self.df._to_geo(na="null", show_bbox=False)
         assert "bbox" not in geo.keys()
@@ -1288,15 +1307,112 @@ class TestConstructor:
         df.columns = pd.MultiIndex.from_product([["geometry"], [0, 1]])
         # don't error in constructor
         gdf = GeoDataFrame(df)
-        # Getting the .geometry column gives GeoDataFrame for both columns
-        # (but with first MultiIndex level removed)
-        # TODO should this give an error instead?
-        result = gdf.geometry
-        assert result.shape == gdf.shape
-        assert result.columns.tolist() == [0, 1]
-        assert_frame_equal(result, gdf["geometry"])
-        result = gdf[["geometry"]]
-        assert_frame_equal(result, gdf if dtype == "geometry" else pd.DataFrame(gdf))
+        with pytest.raises(AttributeError, match=".*geometry .* has not been set.*"):
+            gdf.geometry
+        res_gdf = gdf.set_geometry(("geometry", 0))
+        assert res_gdf.shape == gdf.shape
+        assert isinstance(res_gdf.geometry, GeoSeries)
+
+    def test_default_geo_colname_none(self):
+        match = "You are adding a column named 'geometry' to a GeoDataFrame"
+        gdf = GeoDataFrame({"a": [1, 2]})
+
+        gdf2 = gdf.copy()
+        geo_col = GeoSeries.from_xy([1, 3], [3, 3])
+        with pytest.warns(FutureWarning, match=match):
+            gdf2["geometry"] = geo_col
+        assert gdf2._geometry_column_name == "geometry"
+        gdf4 = gdf.copy()
+        with pytest.warns(FutureWarning, match=match):
+            gdf4.geometry = geo_col
+        assert gdf4._geometry_column_name == "geometry"
+
+        # geo col name should only change if we add geometry
+        gdf5 = gdf.copy()
+        with pytest.warns(
+            UserWarning, match="Geometry column does not contain geometry"
+        ):
+            gdf5["geometry"] = "foo"
+        assert gdf5._geometry_column_name is None
+        gdf3 = gdf.copy().assign(geometry=geo_col)
+        assert gdf3._geometry_column_name == "geometry"
+
+        # Check that adding a GeoSeries to a column called "geometry" to a
+        # gdf without an active geometry column some time after the init does not
+        # warn / set the active geometry column
+        gdf6 = gdf.copy()
+        gdf6["geom2"] = geo_col
+        gdf6["geom3"] = geo_col
+        gdf6 = gdf6.set_geometry("geom2")
+        subset = gdf6[["a", "geom3"]]  # this has a missing active geometry col
+        assert subset._geometry_column_name == "geom2"
+        subset["geometry"] = geo_col
+        # adding column called geometry shouldn't auto-set
+        assert subset._geometry_column_name == "geom2"
+
+    def test_multiindex_geometry_colname_2_level(self):
+        # GH1763 https://github.com/geopandas/geopandas/issues/1763
+        crs = "EPSG:4326"
+        df = pd.DataFrame(
+            [[1, 0], [0, 1]], columns=[["location", "location"], ["x", "y"]]
+        )
+        x_col = df["location", "x"]
+        y_col = df["location", "y"]
+
+        gdf = GeoDataFrame(df, crs=crs, geometry=points_from_xy(x_col, y_col))
+        assert gdf.crs == crs
+        assert gdf.geometry.crs == crs
+        assert gdf.geometry.dtype == "geometry"
+        assert gdf._geometry_column_name == "geometry"
+        assert gdf.geometry.name == "geometry"
+
+    def test_multiindex_geometry_colname_3_level(self):
+        # GH1763 https://github.com/geopandas/geopandas/issues/1763
+        # Note 3-level case uses different code paths in pandas, it is not redundant
+        crs = "EPSG:4326"
+        df = pd.DataFrame(
+            [[1, 0], [0, 1]],
+            columns=[
+                ["foo", "foo"],
+                ["location", "location"],
+                ["x", "y"],
+            ],
+        )
+
+        x_col = df["foo", "location", "x"]
+        y_col = df["foo", "location", "y"]
+
+        gdf = GeoDataFrame(df, crs=crs, geometry=points_from_xy(x_col, y_col))
+        assert gdf.crs == crs
+        assert gdf.geometry.crs == crs
+        assert gdf.geometry.dtype == "geometry"
+        assert gdf._geometry_column_name == "geometry"
+        assert gdf.geometry.name == "geometry"
+
+    def test_multiindex_geometry_colname_3_level_new_col(self):
+        crs = "EPSG:4326"
+        df = pd.DataFrame(
+            [[1, 0], [0, 1]],
+            columns=[
+                ["foo", "foo"],
+                ["location", "location"],
+                ["x", "y"],
+            ],
+        )
+
+        x_col = df["foo", "location", "x"]
+        y_col = df["foo", "location", "y"]
+        df["geometry"] = GeoSeries.from_xy(x_col, y_col)
+        df2 = df.copy()
+        gdf = df.set_geometry("geometry", crs=crs)
+        assert gdf.crs == crs
+        assert gdf._geometry_column_name == "geometry"
+        assert gdf.geometry.name == "geometry"
+        # test again setting with tuple col name
+        gdf = df2.set_geometry(("geometry", "", ""), crs=crs)
+        assert gdf.crs == crs
+        assert gdf._geometry_column_name == ("geometry", "", "")
+        assert gdf.geometry.name == ("geometry", "", "")
 
 
 def test_geodataframe_crs():
diff --git a/geopandas/tests/test_op_output_types.py b/geopandas/tests/test_op_output_types.py
index 365d3627..74a1587e 100644
--- a/geopandas/tests/test_op_output_types.py
+++ b/geopandas/tests/test_op_output_types.py
@@ -1,7 +1,6 @@
 import pandas as pd
 import pyproj
 import pytest
-import geopandas._compat as compat
 
 from shapely.geometry import Point
 import numpy as np
@@ -226,12 +225,11 @@ def test_apply(df):
     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)
-    expected_geo_col_name = geo_name if compat.PANDAS_GE_14 else "geometry"
 
     res = df[["geometry2", "value1"]].apply(identity)
-    assert_obj_no_active_geo_col(res, GeoDataFrame, expected_geo_col_name)
+    assert_obj_no_active_geo_col(res, GeoDataFrame, geo_name)
     assert_obj_no_active_geo_col(
-        df[["geometry2"]].apply(identity), GeoDataFrame, expected_geo_col_name
+        df[["geometry2"]].apply(identity), GeoDataFrame, geo_name
     )
     assert_object(df[["value1"]].apply(identity), pd.DataFrame)
 
@@ -260,7 +258,7 @@ def test_apply(df):
 
 
 def test_apply_axis1_secondary_geo_cols(df):
-    geo_name = df.geometry.name if compat.PANDAS_GE_14 else "geometry"
+    geo_name = df.geometry.name
 
     def identity(x):
         return x
@@ -308,15 +306,13 @@ def test_expanddim_in_unstack():
         index=pd.MultiIndex.from_tuples([("A", "a"), ("A", "b"), ("B", "a")]),
     )
     unstack = s.unstack()
-    expected_geo_col_name = None if compat.PANDAS_GE_12 else "geometry"
-    assert_obj_no_active_geo_col(
-        unstack, GeoDataFrame, geo_colname=expected_geo_col_name
-    )
+    expected_geo_name = None
+    assert_obj_no_active_geo_col(unstack, GeoDataFrame, geo_colname=expected_geo_name)
 
     # https://github.com/geopandas/geopandas/issues/2486
     s.name = "geometry"
     unstack = s.unstack()
-    assert_obj_no_active_geo_col(unstack, GeoDataFrame, expected_geo_col_name)
+    assert_obj_no_active_geo_col(unstack, GeoDataFrame, expected_geo_name)
 
 
 # indexing /  constructor_sliced tests