Repository: geopandas/geopandas

Issue:
BUG: duplicate column names causes wrong geojson output
- [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.
- [ ] (optional) I have confirmed this bug exists on the master branch of geopandas.

---

**Note**: Please read [this guide](https://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports) detailing how to provide the necessary information for us to reproduce your bug.

#### Code Sample, a copy-pastable example

```python
import geopandas as gp

gdf = gp.GeoDataFrame(
    {"a": [1], "b": [2], "c": [3]},
    geometry=gp.points_from_xy(x=[1], y=[1])
)
gdf = gdf[["a", "b", "a", "geometry"]]  # Select a column twice

print("\nLooks good!")
print(gdf)

gdf.to_file("./geo.geojson", driver="GeoJSON")

print(f"\nLooks bad!")
with open("./geo.geojson") as f:
    print(f.read())

```

Output:
```text
Looks good!
   a  b  a                 geometry
0  1  2  1  POINT (1.00000 1.00000)

Looks bad!
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "a": 1, "b": 1 }, "geometry": { "type": "Point", "coordinates": [ 1.0, 1.0 ] } }
]
}
```

#### Problem description

The saved geojson file has the wrong attribute data when the source geodataframe has duplicated column names.
Should be ``{ "a": 1, "b": 2 }`` (or perhaps it should fail to save?).
What you get is the incorrect result ``{ "a": 1, "b": 1 }``

#### Output of ``geopandas.show_versions()``

<details>

SYSTEM INFO
-----------
python     : 3.9.2 (default, Mar  3 2021, 11:58:52)  [Clang 10.0.0 ]
executable : /opt/miniconda3/envs/xxx/bin/python
machine    : macOS-10.15.7-x86_64-i386-64bit
GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.9.1
GEOS lib   : /opt/miniconda3/envs/xxx/lib/libgeos_c.dylib
GDAL       : 3.2.2
GDAL data dir: /opt/miniconda3/envs/xxx/share/gdal
PROJ       : 8.0.0
PROJ data dir: /opt/miniconda3/envs/xxx/lib/python3.9/site-packages/pyproj/proj_dir/share/proj
PYTHON DEPENDENCIES
-------------------
geopandas  : 0.9.0
pandas     : 1.2.3
fiona      : 1.8.18
numpy      : 1.20.1
shapely    : 1.7.1
rtree      : 0.9.7
pyproj     : 3.0.1
matplotlib : 3.3.4
mapclassify: None
geopy      : None
psycopg2   : None
geoalchemy2: None
pyarrow    : None
pygeos     : None
</details>


Discussion:

Patch:
diff --git a/doc/source/getting_started/introduction.ipynb b/doc/source/getting_started/introduction.ipynb
index 66ca3b5c..7f13eb04 100644
--- a/doc/source/getting_started/introduction.ipynb
+++ b/doc/source/getting_started/introduction.ipynb
@@ -2,23 +2,25 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "source": [
     "# Introduction to GeoPandas\n",
     "\n",
-    "This quick tutorial provides an introduction to the key concepts of GeoPandas. In a few minutes, we'll describe the basics which allow you to start your projects.\n",
+    "This quick tutorial introduces the key concepts and basic feautures of GeoPandas to help you get started with your projects.\n",
     "\n",
     "## Concepts\n",
     "\n",
-    "GeoPandas, as the name suggests, extends popular data science library [pandas](https://pandas.pydata.org) by adding support for geospatial data. If you are not familiar with `pandas`, we recommend taking a quick look at its [Getting started documentation](https://pandas.pydata.org/docs/getting_started/index.html#getting-started) before proceeding.\n",
+    "GeoPandas, as the name suggests, extends the popular data science library [pandas](https://pandas.pydata.org) by adding support for geospatial data. If you are not familiar with `pandas`, we recommend taking a quick look at its [Getting started documentation](https://pandas.pydata.org/docs/getting_started/index.html#getting-started) before proceeding.\n",
     "\n",
-    "The core data structure in GeoPandas is `geopandas.GeoDataFrame`, a subclass of `pandas.DataFrame` able to store geometry columns and perform spatial operations. Geometries are handled by `geopandas.GeoSeries`, a subclass of `pandas.Series`. Therefore, your `GeoDataFrame` is a combination of `Series` with your data (numerical, boolean, text etc.) and `GeoSeries` with geometries (points, polygons etc.). You can have as many columns with geometries as you wish, there's no limit typical for desktop GIS software.\n",
+    "The core data structure in GeoPandas is the `geopandas.GeoDataFrame`, a subclass of `pandas.DataFrame`, that can store geometry columns and perform spatial operations. The `geopandas.GeoSeries`, a subclass of `pandas.Series`, handles the geometries. Therefore, your `GeoDataFrame` is a combination of `pandas.Series`, with traditional data (numerical, boolean, text etc.), and `geopandas.GeoSeries`, with geometries (points, polygons etc.). You can have as many columns with geometries as you wish; there's no limit typical for desktop GIS software.\n",
     "\n",
     "![geodataframe schema](../_static/dataframe.svg)\n",
     "\n",
-    "Each `GeoSeries` can contain any geometry type (we can even mix them within a single array) and has a `GeoSeries.crs` attribute, which stores information on the projection (CRS stands for Coordinate Reference System). Therefore, each `GeoSeries` in a `GeoDataFrame` can be in a different projection, allowing you to have, for example, multiple versions of the same geometry, just in a different CRS.\n",
+    "Each `GeoSeries` can contain any geometry type (you can even mix them within a single array) and has a `GeoSeries.crs` attribute, which stores information about the projection (CRS stands for Coordinate Reference System). Therefore, each `GeoSeries` in a `GeoDataFrame` can be in a different projection, allowing you to have, for example, multiple versions (different projections) of the same geometry.\n",
     "\n",
-    "One `GeoSeries` within a `GeoDataFrame` is seen as the _active_ geometry, which means that all geometric operations applied to a `GeoDataFrame` use the specified column.\n",
+    "Only one `GeoSeries` in a `GeoDataFrame` is considered the _active_ geometry, which means that all geometric operations applied to a `GeoDataFrame` operate on this _active_ column.\n",
     "\n",
     "\n",
     "<div class=\"alert alert-info\">\n",
@@ -28,15 +30,15 @@
     "</div>\n",
     "\n",
     "\n",
-    "Let's see how this works in practice.\n",
+    "Let's see how some of these concepts work in practice.\n",
     "\n",
     "## Reading and writing files\n",
     "\n",
     "First, we need to read some data.\n",
     "\n",
-    "### Read files\n",
+    "### Reading files\n",
     "\n",
-    "Assuming we have a file containing both data and geometry (e.g. GeoPackage, GeoJSON, Shapefile), we can easily read it using `geopandas.read_file` function, which automatically detects filetype and creates a `GeoDataFrame`. In this example, we'll use the `\"nybb\"` dataset, a map of New York boroughs which is part of GeoPandas installation. Therefore we need to get the path to the actual file. With your file, you specify a path as a string (`\"my_data/my_file.geojson\"`)."
+    "Assuming you have a file containing both data and geometry (e.g. GeoPackage, GeoJSON, Shapefile), you can read it using `geopandas.read_file()`, which automatically detects the filetype and creates a `GeoDataFrame`. This tutorial uses the `\"nybb\"` dataset, a map of New York boroughs, which is part of the GeoPandas installation. Therefore, we use `geopandas.datasets.get_path()` to retrieve the path to the dataset."
    ]
   },
   {
@@ -55,11 +57,13 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "source": [
-    "### Write files\n",
+    "### Writing files\n",
     "\n",
-    "Writing a `GeoDataFrame` back to file is similarly simple, using `GeoDataFrame.to_file`. The default file format is Shapefile, but you can specify your own using `driver` keyword."
+    "To write a `GeoDataFrame` back to file use `GeoDataFrame.to_file()`. The default file format is Shapefile, but you can specify your own with the `driver` keyword."
    ]
   },
   {
@@ -73,7 +77,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "source": [
     "<div class=\"alert alert-info\">\n",
     "User Guide\n",
@@ -83,17 +89,17 @@
     "\n",
     "\n",
     "\n",
-    "## Simple methods\n",
+    "## Simple accessors and methods\n",
     "\n",
     "Now we have our `GeoDataFrame` and can start working with its geometry. \n",
     "\n",
-    "Since we have only one geometry column read from the file, it is automatically seen as the active geometry and methods used on `GeoDataFrame` will be applied to the `\"geometry\"` column.\n",
+    "Since there was only one geometry column in the New York Boroughs dataset, this column automatically becomes the _active_ geometry and spatial methods used on the `GeoDataFrame` will be applied to the `\"geometry\"` column.\n",
     "\n",
     "### Measuring area\n",
     "\n",
-    "To measure the area of each polygon (or MultiPolygon in this specific case), we can use `GeoDataFrame.area` attribute, which returns a `pandas.Series`. Note that `GeoDataFrame.area` is just `GeoSeries.area` applied to an active geometry column.\n",
+    "To measure the area of each polygon (or MultiPolygon in this specific case), access the `GeoDataFrame.area` attribute, which returns a `pandas.Series`. Note that `GeoDataFrame.area` is just `GeoSeries.area` applied to the _active_ geometry column.\n",
     "\n",
-    "But first, we set the names of boroughs as an index, to make the results easier to read."
+    "But first, to make the results easier to read, set the names of the boroughs as the index:"
    ]
   },
   {
@@ -121,7 +127,7 @@
    "source": [
     "### Getting polygon boundary and centroid\n",
     "\n",
-    "To get just the boundary of each polygon (LineString), we can call `GeoDataFrame.boundary`."
+    "To get the boundary of each polygon (LineString), access the `GeoDataFrame.boundary`:"
    ]
   },
   {
@@ -159,7 +165,7 @@
    "source": [
     "### Measuring distance\n",
     "\n",
-    "We can also measure how far is each centroid from the first one."
+    "We can also measure how far each centroid is from the first centroid location."
    ]
   },
   {
@@ -177,9 +183,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "It's still a DataFrame, so we have all the pandas functionality available to use on the geospatial dataset, and to do data manipulations with the attributes and geometry information together.\n",
+    "Note that `geopandas.GeoDataFrame` is a subclass of `pandas.DataFrame`, so we have all the pandas functionality available to use on the geospatial dataset — we can even perform data manipulations with the attributes and geometry information together.\n",
     "\n",
-    "For example, we can calculate average of the distance measured above (by accessing the `'distance'` column, and calling the `mean()` method on it):"
+    "For example, to calculate the average of the distances measured above, access the 'distance' column and call the mean() method on it:"
    ]
   },
   {
@@ -197,7 +203,7 @@
    "source": [
     "## Making maps\n",
     "\n",
-    "GeoPandas can also plot maps, so we can check how our geometries look like in space. The key method here is `GeoDataFrame.plot()`. In the example below, we plot the `\"area\"` we measured earlier using the active geometry column. We also want to show a legend (`legend=True`)."
+    "GeoPandas can also plot maps, so we can check how the geometries appear in space. To plot the active geometry, call `GeoDataFrame.plot()`. To color code by another column, pass in that column as the first argument. In the example below, we plot the active geometry column and color code by the `\"area\"` column. We also want to show a legend (`legend=True`)."
    ]
   },
   {
@@ -275,7 +281,7 @@
     "\n",
     "### Convex hull\n",
     "\n",
-    "If we are interested in the convex hull of our polygons, we can call `GeoDataFrame.convex_hull`."
+    "If we are interested in the convex hull of our polygons, we can access `GeoDataFrame.convex_hull`."
    ]
   },
   {
@@ -418,7 +424,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "outputs": [],
    "source": [
     "gdf = gdf.set_geometry(\"buffered_centroid\")\n",
@@ -432,7 +440,7 @@
    "source": [
     "## Projections\n",
     "\n",
-    "Each `GeoSeries` has the Coordinate Reference System (CRS) accessible as `GeoSeries.crs`. CRS tells GeoPandas where the coordinates of geometries are located on the Earth. In some cases, CRS is geographic, which means that coordinates are in latitude and longitude. In those cases, its CRS is WGS84, with the authority code `EPSG:4326`. Let's see the projection of our NY boroughs `GeoDataFrame`."
+    "Each `GeoSeries` has its Coordinate Reference System (CRS) accessible at `GeoSeries.crs`. The CRS tells GeoPandas where the coordinates of the geometries are located on the earth's surface. In some cases, the CRS is geographic, which means that the coordinates are in latitude and longitude. In those cases, its CRS is WGS84, with the authority code `EPSG:4326`. Let's see the projection of our NY boroughs `GeoDataFrame`."
    ]
   },
   {
@@ -475,12 +483,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Notice the difference in coordinates along the axes of the plot. Where we had 120 000 - 280 000 (feet) before, we have 40.5 - 40.9 (degrees) now. In this case, `boroughs_4326` has a `\"geometry\"` column in WGS84 but all the other (with centroids etc.) remains in the original CRS.\n",
+    "Notice the difference in coordinates along the axes of the plot. Where we had 120 000 - 280 000 (feet) before, we now have 40.5 - 40.9 (degrees). In this case, `boroughs_4326` has a `\"geometry\"` column in WGS84 but all the other (with centroids etc.) remain in the original CRS.\n",
     "\n",
     "<div class=\"alert alert-warning\">\n",
     "Warning\n",
     "    \n",
-    "For operations that rely on distance or area, you always need to use projected CRS (in meters, feet, kilometers etc.) not a geographic one. GeoPandas operations are planar, and degrees reflect the position on a sphere. Therefore the results may not be correct. For example, the result of `gdf.area.sum()` (projected CRS) is 8 429 911 572 ft<sup>2</sup> but the result of `boroughs_4326.area.sum()` (geographic CRS) is 0.083.\n",
+    "For operations that rely on distance or area, you always need to use a projected CRS (in meters, feet, kilometers etc.) not a geographic one (in degrees). GeoPandas operations are planar, whereas degrees reflect the position on a sphere. Therefore, spatial operations using degrees may not yield correct results. For example, the result of `gdf.area.sum()` (projected CRS) is 8 429 911 572 ft<sup>2</sup> but the result of `boroughs_4326.area.sum()` (geographic CRS) is 0.083.\n",
     "</div>\n",
     "\n",
     "<div class=\"alert alert-info\">\n",
@@ -491,17 +499,10 @@
     "\n",
     "## What next?\n",
     "\n",
-    "With GeoPandas we can do much more that this, from [aggregations](../docs/user_guide/aggregation_with_dissolve.rst), to [spatial joins](../docs/user_guide/mergingdata.rst), [geocoding](../docs/user_guide/geocoding.rst) and [much more](../gallery/index.rst).\n",
+    "With GeoPandas we can do much more than what has been introduced so far, from [aggregations](../docs/user_guide/aggregation_with_dissolve.rst), to [spatial joins](../docs/user_guide/mergingdata.rst), to [geocoding](../docs/user_guide/geocoding.rst), and [much more](../gallery/index.rst).\n",
     "\n",
-    "Head to the [User Guide](../docs/user_guide.rst) for to learn more about different functionality of GeoPandas, to the [Examples](../gallery/index.rst) to see how it can be used or the the [API reference](../docs/reference.rst) for the details."
+    "Head over to the [User Guide](../docs/user_guide.rst) to learn more about the different features of GeoPandas, the [Examples](../gallery/index.rst) to see how they can be used, or to the [API reference](../docs/reference.rst) for the details."
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -515,7 +516,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.6"
+   "version": "3.9.1"
   }
  },
  "nbformat": 4,
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 3184c94b..2ec09f0f 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -7,6 +7,7 @@ Uses PyGEOS if available/set, otherwise loops through Shapely geometries.
 import warnings
 
 import numpy as np
+import pandas as pd
 
 import shapely.geometry
 import shapely.geos
@@ -44,9 +45,9 @@ else:
     type_mapping, geometry_type_ids, geometry_type_values = None, None, None
 
 
-def _isna(value):
+def isna(value):
     """
-    Check if scalar value is NA-like (None or np.nan).
+    Check if scalar value is NA-like (None, np.nan or pd.NA).
 
     Custom version that only works for scalars (returning True or False),
     as `pd.isna` also works for array-like input returning a boolean array.
@@ -55,6 +56,8 @@ def _isna(value):
         return True
     elif isinstance(value, float) and np.isnan(value):
         return True
+    elif compat.PANDAS_GE_10 and value is pd.NA:
+        return True
     else:
         return False
 
@@ -127,7 +130,7 @@ def from_shapely(data):
                 out.append(_shapely_to_pygeos(geom))
             else:
                 out.append(geom)
-        elif _isna(geom):
+        elif isna(geom):
             out.append(None)
         else:
             raise TypeError("Input must be valid geometry objects: {0}".format(geom))
@@ -165,7 +168,7 @@ def from_wkb(data):
     out = []
 
     for geom in data:
-        if geom is not None and len(geom):
+        if not isna(geom) and len(geom):
             geom = shapely.wkb.loads(geom)
         else:
             geom = None
@@ -200,7 +203,7 @@ def from_wkt(data):
     out = []
 
     for geom in data:
-        if geom is not None and len(geom):
+        if not isna(geom) and len(geom):
             if isinstance(geom, bytes):
                 geom = geom.decode("utf-8")
             geom = shapely.wkt.loads(geom)
@@ -941,7 +944,7 @@ def transform(data, func):
         result = np.empty(n, dtype=object)
         for i in range(n):
             geom = data[i]
-            if _isna(geom):
+            if isna(geom):
                 result[i] = geom
             else:
                 result[i] = transform(func, geom)
diff --git a/geopandas/array.py b/geopandas/array.py
index e35676c2..d2b1ea2a 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -56,23 +56,6 @@ class GeometryDtype(ExtensionDtype):
 register_extension_dtype(GeometryDtype)
 
 
-def _isna(value):
-    """
-    Check if scalar value is NA-like (None, np.nan or pd.NA).
-
-    Custom version that only works for scalars (returning True or False),
-    as `pd.isna` also works for array-like input returning a boolean array.
-    """
-    if value is None:
-        return True
-    elif isinstance(value, float) and np.isnan(value):
-        return True
-    elif compat.PANDAS_GE_10 and value is pd.NA:
-        return True
-    else:
-        return False
-
-
 def _check_crs(left, right, allow_none=False):
     """
     Check if the projection of both arrays is the same.
@@ -398,8 +381,8 @@ class GeometryArray(ExtensionArray):
             if isinstance(key, numbers.Integral):
                 raise ValueError("cannot set a single element with an array")
             self.data[key] = value.data
-        elif isinstance(value, BaseGeometry) or _isna(value):
-            if _isna(value):
+        elif isinstance(value, BaseGeometry) or vectorized.isna(value):
+            if vectorized.isna(value):
                 # internally only use None as missing value indicator
                 # but accept others
                 value = None
@@ -1005,7 +988,7 @@ class GeometryArray(ExtensionArray):
 
         if mask.any():
             # fill with value
-            if _isna(value):
+            if vectorized.isna(value):
                 value = None
             elif not isinstance(value, BaseGeometry):
                 raise NotImplementedError(
@@ -1326,7 +1309,7 @@ class GeometryArray(ExtensionArray):
         """
         Return for `item in self`.
         """
-        if _isna(item):
+        if vectorized.isna(item):
             if (
                 item is self.dtype.na_value
                 or isinstance(item, self.dtype.type)
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 0b2848ea..505dc406 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -810,6 +810,9 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
         ids = np.array(self.index, copy=False)
         geometries = np.array(self[self._geometry_column_name], copy=False)
 
+        if not self.columns.is_unique:
+            raise ValueError("GeoDataFrame cannot contain duplicated column names.")
+
         properties_cols = self.columns.difference([self._geometry_column_name])
 
         if len(properties_cols) > 0:


Test patch:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index ef3cd143..3770ca27 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -139,15 +139,16 @@ def test_from_wkb():
     assert all(v.equals(t) for v, t in zip(res, points_no_missing))
 
     # missing values
-    # TODO(pygeos) does not support empty strings
-    if compat.USE_PYGEOS:
-        L_wkb.extend([None])
-    else:
-        L_wkb.extend([b"", None])
-    res = from_wkb(L_wkb)
-    assert res[-1] is None
+    # TODO(pygeos) does not support empty strings, np.nan, or pd.NA
+    missing_values = [None]
     if not compat.USE_PYGEOS:
-        assert res[-2] is None
+        missing_values.extend([b"", np.nan])
+
+        if compat.PANDAS_GE_10:
+            missing_values.append(pd.NA)
+
+    res = from_wkb(missing_values)
+    np.testing.assert_array_equal(res, np.full(len(missing_values), None))
 
     # single MultiPolygon
     multi_poly = shapely.geometry.MultiPolygon(
@@ -202,15 +203,16 @@ def test_from_wkt(string_type):
     assert all(v.almost_equals(t) for v, t in zip(res, points_no_missing))
 
     # missing values
-    # TODO(pygeos) does not support empty strings
-    if compat.USE_PYGEOS:
-        L_wkt.extend([None])
-    else:
-        L_wkt.extend([f(""), None])
-    res = from_wkt(L_wkt)
-    assert res[-1] is None
+    # TODO(pygeos) does not support empty strings, np.nan, or pd.NA
+    missing_values = [None]
     if not compat.USE_PYGEOS:
-        assert res[-2] is None
+        missing_values.extend([f(""), np.nan])
+
+        if compat.PANDAS_GE_10:
+            missing_values.append(pd.NA)
+
+    res = from_wkb(missing_values)
+    np.testing.assert_array_equal(res, np.full(len(missing_values), None))
 
     # single MultiPolygon
     multi_poly = shapely.geometry.MultiPolygon(
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 62cb23e9..36781142 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -457,6 +457,15 @@ class TestDataFrame:
         for f in data["features"]:
             assert "id" not in f.keys()
 
+    def test_to_json_with_duplicate_columns(self):
+        df = GeoDataFrame(
+            data=[[1, 2, 3]], columns=["a", "b", "a"], geometry=[Point(1, 1)]
+        )
+        with pytest.raises(
+            ValueError, match="GeoDataFrame cannot contain duplicated column names."
+        ):
+            df.to_json()
+
     def test_copy(self):
         df2 = self.df.copy()
         assert type(df2) is GeoDataFrame
@@ -485,6 +494,16 @@ class TestDataFrame:
         df = GeoDataFrame.from_file(tempfilename)
         assert df.crs == "epsg:2263"
 
+    def test_to_file_with_duplicate_columns(self):
+        df = GeoDataFrame(
+            data=[[1, 2, 3]], columns=["a", "b", "a"], geometry=[Point(1, 1)]
+        )
+        with pytest.raises(
+            ValueError, match="GeoDataFrame cannot contain duplicated column names."
+        ):
+            tempfilename = os.path.join(self.tempdir, "crs.shp")
+            df.to_file(tempfilename)
+
     def test_bool_index(self):
         # Find boros with 'B' in their name
         df = self.df[self.df["BoroName"].str.contains("B")]
@@ -666,6 +685,14 @@ class TestDataFrame:
         result = list(df_only_numerical_cols.iterfeatures(na="keep"))[0]
         assert type(result["properties"]["Shape_Leng"]) is float
 
+        with pytest.raises(
+            ValueError, match="GeoDataFrame cannot contain duplicated column names."
+        ):
+            df_with_duplicate_columns = df[
+                ["Shape_Leng", "Shape_Leng", "Shape_Area", "geometry"]
+            ]
+            list(df_with_duplicate_columns.iterfeatures())
+
         # geometry not set
         df = GeoDataFrame({"values": [0, 1], "geom": [Point(0, 1), Point(1, 0)]})
         with pytest.raises(AttributeError):