Repository: geopandas/geopandas
Pull request #3275 (2024-05-06 22:23:54)

ENH: support writing the geoarrow-based encodings of GeoParquet
~WIP, on top of https://github.com/geopandas/geopandas/pull/3219~

~Will probably break out the reading side separately as a first PR, since that does not rely on https://github.com/geopandas/geopandas/pull/3219~ (reading side: https://github.com/geopandas/geopandas/pull/3278)

This implements the writing side for geoarrow encodings in GeoParquet (https://github.com/geopandas/geopandas/issues/3253)


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7a9c7e5a..9375b628 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -50,18 +50,21 @@ New features and improvements:
   the specified columns (#3101).
 - Added support to ``read_file`` for the ``mask`` keyword for the pyogrio engine (#3062).
 - Added support to ``read_file`` for the ``columns`` keyword for the fiona engine (#3133).
-- Added support to ``read_parquet`` for reading files using the GeoArrow-based native geometry encoding of GeoParquet 1.1 (#3253).
+- Added support to ``to_parquet`` and ``read_parquet`` for writing and reading files
+  using the GeoArrow-based native geometry encoding of GeoParquet 1.1 (#3253, #3275).
 - Add `sort` keyword to `clip` method for GeoSeries and GeoDataFrame to allow optional
   preservation of the original order of observations. (#3233)
 - Added `show_bbox`, `drop_id` and `to_wgs84` arguments to allow further customization of
   `GeoSeries.to_json` (#3226)
-- `explore` now supports `GeoDataFrame`s with additional columns containing datetimes, uuids and 
+- `explore` now supports `GeoDataFrame`s with additional columns containing datetimes, uuids and
   other non JSON serializable objects (#3261).
 - The `GeoSeries.fillna` method now supports the `limit` keyword (#3290).
 - Added support for `bbox` covering encoding in geoparquet. Can filter reading of parquet
 files based on a bounding box, and write out a bounding box column to parquet files (#3282)
 - `align` keyword in binary methods now defaults to `None`, treated as True. Explicit True
   will silence the warning about mismachted indices. (#3212)
+- `GeoSeries.set_crs` can now be used to remove CRS information by passing
+  `crs=None, allow_override=True`. (#3316)
 
 Backwards incompatible API changes:
 
@@ -127,6 +130,8 @@ New deprecations:
   ```
 - The `geopandas.use_pygeos` option has been deprecated and will be removed in GeoPandas
   1.1 (#3283)
+- Manual overriding of an existing CRS of a GeoSeries using the `crs.setter` has been deprecated
+  and will be disabled in future. Use the `set_crs()` method instead (#3085).
 
 Bug fixes:
 
diff --git a/doc/environment.yml b/doc/environment.yml
index eea05d5e..7d37d9e5 100644
--- a/doc/environment.yml
+++ b/doc/environment.yml
@@ -5,9 +5,6 @@ dependencies:
   - python
   - pandas
   - shapely
-  - fiona
-  - pyproj
-  - pyogrio
   - geopy
   - matplotlib
   - mapclassify
@@ -17,30 +14,30 @@ dependencies:
   - ipython
   - pillow
   - mock
-  - cartopy
   - pyepsg
-  - contextily
-  - rasterio
-  - geoplot
   - sphinx-gallery
   - jinja2
   - doc2dash
   - matplotlib-scalebar
   # specify additional dependencies to reduce solving for conda
-  - gdal
-  - libgdal
-  - proj
   - geos
   - nbsphinx
   - jupyter_client
   - ipykernel
   - myst-parser
   - folium
-  - libpysal
   - xyzservices
   - packaging
-  - pointpats
   - pip
   - pip:
       - sphinx-toggleprompt
       - geodatasets
+      # proj/pyproj deps from pip to avoid Proj env warning from conda+readthedocs
+      - pyproj
+      - cartopy
+      - fiona
+      - pyogrio
+      - rasterio
+      - contextily
+      - pointpats
+      - libpysal
diff --git a/doc/source/conf.py b/doc/source/conf.py
index f8ddd866..da5107a5 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -383,8 +383,8 @@ intersphinx_mapping = {
         "https://fiona.readthedocs.io/en/stable/objects.inv",
     ),
     "folium": (
-        "https://python-visualization.github.io/folium/",
-        "https://python-visualization.github.io/folium/objects.inv",
+        "https://python-visualization.github.io/folium/latest",
+        "https://python-visualization.github.io/folium/latest/objects.inv",
     ),
     "geoplot": (
         "https://residentmario.github.io/geoplot/index.html",
diff --git a/doc/source/docs/user_guide/fiona_to_pyogrio.rst b/doc/source/docs/user_guide/fiona_to_pyogrio.rst
index 97869489..4d53c4e2 100644
--- a/doc/source/docs/user_guide/fiona_to_pyogrio.rst
+++ b/doc/source/docs/user_guide/fiona_to_pyogrio.rst
@@ -62,6 +62,7 @@ Pyogrio writes EMPTY and None geometries as such to e.g. GPKG files, Fiona write
 as None.
 
 .. ipython:: python
+    :okwarning:
 
     import shapely
     
diff --git a/environment.yml b/environment.yml
index 44b7251f..5b3480dd 100644
--- a/environment.yml
+++ b/environment.yml
@@ -14,6 +14,5 @@ dependencies:
   - pyepsg
   - contextily
   - rasterio
-  - geoplot
   - folium
   - packaging
diff --git a/geopandas/array.py b/geopandas/array.py
index cb394063..32f338a7 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -851,7 +851,7 @@ class GeometryArray(ExtensionArray):
         warnings.warn(
             "The 'unary_union' attribute is deprecated, "
             "use the 'union_all' method instead.",
-            FutureWarning,
+            DeprecationWarning,
             stacklevel=2,
         )
         return self.union_all()
diff --git a/geopandas/base.py b/geopandas/base.py
index 42688dea..f33cdd0b 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2108,7 +2108,7 @@ GeometryCollection
         warn(
             "The 'unary_union' attribute is deprecated, "
             "use the 'union_all()' method instead.",
-            FutureWarning,
+            DeprecationWarning,
             stacklevel=2,
         )
 
@@ -4539,7 +4539,7 @@ GeometryCollection
         The operation works on a 1-to-1 row-wise manner:
 
         .. image:: ../../../_static/binary_op-01.svg
-        :align: center
+            :align: center
 
         Parameters
         ----------
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index c745278a..19e7f2ee 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -53,7 +53,10 @@ def _ensure_geometry(data, crs=None):
         if data.crs is None and crs is not None:
             # Avoids caching issues/crs sharing issues
             data = data.copy()
-            data.crs = crs
+            if isinstance(data, GeometryArray):
+                data.crs = crs
+            else:
+                data.array.crs = crs
         return data
     else:
         if isinstance(data, Series):
@@ -388,7 +391,10 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         # Check that we are using a listlike of geometries
         level = _ensure_geometry(level, crs=crs)
         # ensure_geometry only sets crs on level if it has crs==None
-        level.crs = crs
+        if isinstance(level, GeoSeries):
+            level.array.crs = crs
+        else:
+            level.crs = crs
         # update _geometry_column_name prior to assignment
         # to avoid default is None warning
         frame._geometry_column_name = geo_column_name
@@ -515,6 +521,14 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
             )
 
         if hasattr(self.geometry.values, "crs"):
+            if self.crs is not None:
+                warnings.warn(
+                    "Overriding the CRS of a GeoDataFrame that already has CRS. "
+                    "This unsafe behavior will be deprecated in future versions. "
+                    "Use GeoDataFrame.set_crs method instead",
+                    stacklevel=2,
+                    category=DeprecationWarning,
+                )
             self.geometry.values.crs = value
         else:
             # column called 'geometry' without geometry
@@ -1228,7 +1242,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         """
         from geopandas.io.geoarrow import ArrowTable, geopandas_to_arrow
 
-        table = geopandas_to_arrow(
+        table, _ = geopandas_to_arrow(
             self,
             index=index,
             geometry_encoding=geometry_encoding,
@@ -1242,8 +1256,9 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         path,
         index=None,
         compression="snappy",
-        schema_version=None,
+        geometry_encoding="WKB",
         write_covering_bbox=False,
+        schema_version=None,
         **kwargs,
     ):
         """Write a GeoDataFrame to the Parquet format.
@@ -1265,6 +1280,10 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             output except `RangeIndex` which is stored as metadata only.
         compression : {'snappy', 'gzip', 'brotli', None}, default 'snappy'
             Name of the compression to use. Use ``None`` for no compression.
+        geometry_encoding : {'WKB', 'geoarrow'}, default 'WKB'
+            The encoding to use for the geometry columns. Defaults to "WKB"
+            for maximum interoperability. Specify "geoarrow" to use one of the
+            native GeoArrow-based single-geometry type encodings.
         schema_version : {'0.1.0', '0.4.0', '1.0.0', None}
             GeoParquet specification version; if not provided will default to
             latest supported version.
@@ -1303,6 +1322,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             self,
             path,
             compression=compression,
+            geometry_encoding=geometry_encoding,
             index=index,
             schema_version=schema_version,
             write_covering_bbox=write_covering_bbox,
@@ -1463,16 +1483,20 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         If there are multiple geometry columns within the GeoDataFrame, only
         the CRS of the active geometry column is set.
 
-        NOTE: The underlying geometries are not transformed to this CRS. To
+        Pass ``None`` to remove CRS from the active geometry column.
+
+        Notes
+        -----
+        The underlying geometries are not transformed to this CRS. To
         transform the geometries to a new CRS, use the ``to_crs`` method.
 
         Parameters
         ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
+        crs : pyproj.CRS | None, optional
             The value can be anything accepted
             by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
             such as an authority string (eg "EPSG:4326") or a WKT string.
-        epsg : int, optional if `crs` is specified
+        epsg : int, optional
             EPSG code specifying the projection.
         inplace : bool, default False
             If True, the CRS of the GeoDataFrame will be changed in place
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 7ae43bac..505304e7 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -236,6 +236,18 @@ class GeoSeries(GeoPandasBase, Series):
     def append(self, *args, **kwargs) -> GeoSeries:
         return self._wrapped_pandas_method("append", *args, **kwargs)
 
+    @GeoPandasBase.crs.setter
+    def crs(self, value):
+        if self.crs is not None:
+            warnings.warn(
+                "Overriding the CRS of a GeoSeries that already has CRS. "
+                "This unsafe behavior will be deprecated in future versions. "
+                "Use GeoSeries.set_crs method instead.",
+                stacklevel=2,
+                category=DeprecationWarning,
+            )
+        self.geometry.values.crs = value
+
     @property
     def geometry(self) -> GeoSeries:
         return self
@@ -642,7 +654,6 @@ class GeoSeries(GeoPandasBase, Series):
         from geopandas import GeoDataFrame
 
         data = GeoDataFrame({"geometry": self}, index=self.index)
-        data.crs = self.crs
         data.to_file(filename, driver, index=index, **kwargs)
 
     #
@@ -987,12 +998,16 @@ class GeoSeries(GeoPandasBase, Series):
         """
         Set the Coordinate Reference System (CRS) of a ``GeoSeries``.
 
-        NOTE: The underlying geometries are not transformed to this CRS. To
+        Pass ``None`` to remove CRS from the ``GeoSeries``.
+
+        Notes
+        -----
+        The underlying geometries are not transformed to this CRS. To
         transform the geometries to a new CRS, use the ``to_crs`` method.
 
         Parameters
         ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
+        crs : pyproj.CRS | None, optional
             The value can be anything accepted
             by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
             such as an authority string (eg "EPSG:4326") or a WKT string.
@@ -1060,8 +1075,6 @@ class GeoSeries(GeoPandasBase, Series):
             crs = CRS.from_user_input(crs)
         elif epsg is not None:
             crs = CRS.from_epsg(epsg)
-        else:
-            raise ValueError("Must pass either crs or epsg.")
 
         if not allow_override and self.crs is not None and not self.crs == crs:
             raise ValueError(
@@ -1074,7 +1087,7 @@ class GeoSeries(GeoPandasBase, Series):
             result = self.copy()
         else:
             result = self
-        result.crs = crs
+        result.array.crs = crs
         return result
 
     def to_crs(
diff --git a/geopandas/io/arrow.py b/geopandas/io/arrow.py
index 756aedc8..9fcebeb7 100644
--- a/geopandas/io/arrow.py
+++ b/geopandas/io/arrow.py
@@ -25,6 +25,7 @@ GEOARROW_ENCODINGS = [
     "multipolygon",
 ]
 SUPPORTED_ENCODINGS = ["WKB"] + GEOARROW_ENCODINGS
+
 # reference: https://github.com/opengeospatial/geoparquet
 
 # Metadata structure:
@@ -110,7 +111,9 @@ def _get_geometry_types(series):
     return sorted([_geometry_type_names[idx] for idx in geometry_types])
 
 
-def _create_metadata(df, schema_version=None, write_covering_bbox=False):
+def _create_metadata(
+    df, schema_version=None, geometry_encoding=None, write_covering_bbox=False
+):
     """Create and encode geo metadata dict.
 
     Parameters
@@ -128,8 +131,13 @@ def _create_metadata(df, schema_version=None, write_covering_bbox=False):
     -------
     dict
     """
-
-    schema_version = schema_version or METADATA_VERSION
+    if schema_version is None:
+        if geometry_encoding and any(
+            encoding != "WKB" for encoding in geometry_encoding.values()
+        ):
+            schema_version = "1.1.0"
+        else:
+            schema_version = METADATA_VERSION
 
     if schema_version not in SUPPORTED_VERSIONS:
         raise ValueError(
@@ -158,7 +166,7 @@ def _create_metadata(df, schema_version=None, write_covering_bbox=False):
                 _remove_id_from_member_of_ensembles(crs)
 
         column_metadata[col] = {
-            "encoding": "WKB",
+            "encoding": geometry_encoding[col],
             "crs": crs,
             geometry_types_name: geometry_types,
         }
@@ -181,7 +189,7 @@ def _create_metadata(df, schema_version=None, write_covering_bbox=False):
     return {
         "primary_column": df._geometry_column_name,
         "columns": column_metadata,
-        "version": schema_version or METADATA_VERSION,
+        "version": schema_version,
         "creator": {"library": "geopandas", "version": geopandas.__version__},
     }
 
@@ -315,7 +323,13 @@ def _validate_geo_metadata(metadata):
                         raise ValueError("Metadata for bbox column is malformed.")
 
 
-def _geopandas_to_arrow(df, index=None, schema_version=None, write_covering_bbox=None):
+def _geopandas_to_arrow(
+    df,
+    index=None,
+    geometry_encoding="WKB",
+    schema_version=None,
+    write_covering_bbox=None,
+):
     """
     Helper function with main, shared logic for to_parquet/to_feather.
     """
@@ -325,13 +339,20 @@ def _geopandas_to_arrow(df, index=None, schema_version=None, write_covering_bbox
 
     _validate_dataframe(df)
 
-    # create geo metadata before altering incoming data frame
-    geo_metadata = _create_metadata(
-        df, schema_version=schema_version, write_covering_bbox=write_covering_bbox
-    )
+    if schema_version is not None:
+        if geometry_encoding != "WKB" and schema_version != "1.1.0":
+            raise ValueError(
+                "'geoarrow' encoding is only supported with schema version >= 1.1.0"
+            )
 
-    table = geopandas_to_arrow(
-        df, geometry_encoding="WKB", index=index, interleaved=True
+    table, geometry_encoding_dict = geopandas_to_arrow(
+        df, geometry_encoding=geometry_encoding, index=index, interleaved=False
+    )
+    geo_metadata = _create_metadata(
+        df,
+        schema_version=schema_version,
+        geometry_encoding=geometry_encoding_dict,
+        write_covering_bbox=write_covering_bbox,
     )
 
     if write_covering_bbox:
@@ -355,6 +376,7 @@ def _to_parquet(
     path,
     index=None,
     compression="snappy",
+    geometry_encoding="WKB",
     schema_version=None,
     write_covering_bbox=False,
     **kwargs,
@@ -383,6 +405,10 @@ def _to_parquet(
         output except `RangeIndex` which is stored as metadata only.
     compression : {'snappy', 'gzip', 'brotli', None}, default 'snappy'
         Name of the compression to use. Use ``None`` for no compression.
+    geometry_encoding : {'WKB', 'geoarrow'}, default 'WKB'
+        The encoding to use for the geometry columns. Defaults to "WKB"
+        for maximum interoperability. Specify "geoarrow" to use one of the
+        native GeoArrow-based single-geometry type encodings.
     schema_version : {'0.1.0', '0.4.0', '1.0.0', None}
         GeoParquet specification version; if not provided will default to
         latest supported version.
@@ -412,6 +438,7 @@ def _to_parquet(
     table = _geopandas_to_arrow(
         df,
         index=index,
+        geometry_encoding=geometry_encoding,
         schema_version=schema_version,
         write_covering_bbox=write_covering_bbox,
     )
diff --git a/geopandas/io/geoarrow.py b/geopandas/io/geoarrow.py
index e67d1dc1..004152ca 100644
--- a/geopandas/io/geoarrow.py
+++ b/geopandas/io/geoarrow.py
@@ -121,6 +121,8 @@ def geopandas_to_arrow(
 
     table = pa.Table.from_pandas(df_attr, preserve_index=index)
 
+    geometry_encoding_dict = {}
+
     if geometry_encoding.lower() == "geoarrow":
         if Version(pa.__version__) < Version("10.0.0"):
             raise ValueError("Converting to 'geoarrow' requires pyarrow >= 10.0.")
@@ -135,6 +137,11 @@ def geopandas_to_arrow(
                 interleaved=interleaved,
             )
             table = table.set_column(i, field, geom_arr)
+            geometry_encoding_dict[col] = (
+                field.metadata[b"ARROW:extension:name"]
+                .decode()
+                .removeprefix("geoarrow.")
+            )
 
     elif geometry_encoding.lower() == "wkb":
         # Encode all geometry columns to WKB
@@ -143,12 +150,13 @@ def geopandas_to_arrow(
                 np.asarray(df[col].array), field_name=col, crs=df[col].crs
             )
             table = table.set_column(i, field, wkb_arr)
+            geometry_encoding_dict[col] = "WKB"
 
     else:
         raise ValueError(
             f"Expected geometry encoding 'WKB' or 'geoarrow' got {geometry_encoding}"
         )
-    return table
+    return table, geometry_encoding_dict
 
 
 def construct_wkb_array(
diff --git a/readthedocs.yml b/readthedocs.yml
index a92d29fb..4f6b06fd 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -8,6 +8,8 @@ build:
       # we need the tags for versioneer to work
       - git fetch origin --depth 150
       - git fetch --tags
+    post_build:
+      - mamba list
 python:
   install:
   - method: pip


Test patch:
diff --git a/geopandas/io/tests/test_arrow.py b/geopandas/io/tests/test_arrow.py
index 86cc3078..94b58122 100644
--- a/geopandas/io/tests/test_arrow.py
+++ b/geopandas/io/tests/test_arrow.py
@@ -66,7 +66,7 @@ def file_format(request):
 
 def test_create_metadata(naturalearth_lowres):
     df = read_file(naturalearth_lowres)
-    metadata = _create_metadata(df)
+    metadata = _create_metadata(df, geometry_encoding={"geometry": "WKB"})
 
     assert isinstance(metadata, dict)
     assert metadata["version"] == METADATA_VERSION
@@ -89,6 +89,10 @@ def test_create_metadata(naturalearth_lowres):
     assert metadata["creator"]["library"] == "geopandas"
     assert metadata["creator"]["version"] == geopandas.__version__
 
+    # specifying non-WKB encoding sets default schema to 1.1.0
+    metadata = _create_metadata(df, geometry_encoding={"geometry": "point"})
+    assert metadata["version"] == "1.1.0"
+
 
 def test_create_metadata_with_z_geometries():
     geometry_types = [
@@ -133,18 +137,18 @@ def test_create_metadata_with_z_geometries():
             ],
         },
     )
-    metadata = _create_metadata(df)
+    metadata = _create_metadata(df, geometry_encoding={"geometry": "WKB"})
     assert sorted(metadata["columns"]["geometry"]["geometry_types"]) == sorted(
         geometry_types
     )
     # only 3D geometries
-    metadata = _create_metadata(df.iloc[1::2])
+    metadata = _create_metadata(df.iloc[1::2], geometry_encoding={"geometry": "WKB"})
     assert all(
         geom_type.endswith(" Z")
         for geom_type in metadata["columns"]["geometry"]["geometry_types"]
     )
 
-    metadata = _create_metadata(df.iloc[5:7])
+    metadata = _create_metadata(df.iloc[5:7], geometry_encoding={"geometry": "WKB"})
     assert metadata["columns"]["geometry"]["geometry_types"] == [
         "MultiPolygon",
         "Polygon Z",
@@ -169,11 +173,17 @@ def test_crs_metadata_datum_ensemble():
     assert pyproj.CRS(crs_json) == crs
 
 
-def test_write_metadata_invalid_spec_version():
+def test_write_metadata_invalid_spec_version(tmp_path):
     gdf = geopandas.GeoDataFrame(geometry=[box(0, 0, 10, 10)], crs="EPSG:4326")
     with pytest.raises(ValueError, match="schema_version must be one of"):
         _create_metadata(gdf, schema_version="invalid")
 
+    with pytest.raises(
+        ValueError,
+        match="'geoarrow' encoding is only supported with schema version >= 1.1.0",
+    ):
+        gdf.to_parquet(tmp_path, schema_version="1.0.0", geometry_encoding="geoarrow")
+
 
 def test_encode_metadata():
     metadata = {"a": "b"}
@@ -340,6 +350,7 @@ def test_to_parquet_does_not_pass_engine_along(mock_to_parquet):
         df,
         "",
         compression="snappy",
+        geometry_encoding="WKB",
         index=None,
         schema_version=None,
         write_covering_bbox=False,
@@ -636,7 +647,7 @@ def test_missing_crs(tmpdir, file_format, naturalearth_lowres):
     reader, writer = file_format
 
     df = read_file(naturalearth_lowres)
-    df.crs = None
+    df.geometry.array.crs = None
 
     filename = os.path.join(str(tmpdir), "test.pq")
     writer(df, filename)
@@ -981,7 +992,6 @@ def test_parquet_read_partitioned_dataset_fsspec(tmpdir, naturalearth_lowres):
     ["point", "linestring", "polygon", "multipoint", "multilinestring", "multipolygon"],
 )
 def test_read_parquet_geoarrow(geometry_type):
-
     result = geopandas.read_parquet(
         DATA_PATH
         / "arrow"
@@ -997,6 +1007,24 @@ def test_read_parquet_geoarrow(geometry_type):
     assert_geodataframe_equal(result, expected, check_crs=True)
 
 
+@pytest.mark.parametrize(
+    "geometry_type",
+    ["point", "linestring", "polygon", "multipoint", "multilinestring", "multipolygon"],
+)
+def test_geoarrow_roundtrip(tmp_path, geometry_type):
+
+    df = geopandas.read_parquet(
+        DATA_PATH
+        / "arrow"
+        / "geoparquet"
+        / f"data-{geometry_type}-encoding_wkb.parquet"
+    )
+
+    df.to_parquet(tmp_path / "test.parquet", geometry_encoding="geoarrow")
+    result = geopandas.read_parquet(tmp_path / "test.parquet")
+    assert_geodataframe_equal(result, df, check_crs=True)
+
+
 def test_to_parquet_bbox_structure_and_metadata(tmpdir, naturalearth_lowres):
     # check metadata being written for covering.
     from pyarrow import parquet
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index d17f28a0..4d79e16a 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -461,8 +461,7 @@ def test_to_file_crs(tmpdir, engine, nybb_filename):
     assert result.crs == "epsg:3857"
 
     # specify CRS for gdf without one
-    df2 = df.copy()
-    df2.crs = None
+    df2 = df.set_crs(None, allow_override=True)
     df2.to_file(tempfilename, crs=2263, engine=engine)
     df = GeoDataFrame.from_file(tempfilename, engine=engine)
     assert df.crs == "epsg:2263"
@@ -555,7 +554,7 @@ def test_empty_crs(tmpdir, driver, ext, engine):
 
     if ext == ".geojson":
         # geojson by default assumes epsg:4326
-        df.crs = "EPSG:4326"
+        df.geometry.array.crs = "EPSG:4326"
 
     assert_geodataframe_equal(result, df)
 
diff --git a/geopandas/io/tests/test_geoarrow.py b/geopandas/io/tests/test_geoarrow.py
index 68881b6f..c9502264 100644
--- a/geopandas/io/tests/test_geoarrow.py
+++ b/geopandas/io/tests/test_geoarrow.py
@@ -36,7 +36,6 @@ def pa_array(array):
 
 
 def assert_table_equal(left, right, check_metadata=True):
-
     geom_type = left["geometry"].type
     # in case of Points (directly the inner fixed_size_list or struct type)
     # -> there are NaNs for empties -> we need to compare them separately
@@ -154,7 +153,7 @@ def test_geoarrow_export(geometry_type, dim, geometry_encoding, interleaved):
     df["geometry"] = GeoSeries.from_wkb(df["geometry"])
     df["row_number"] = df["row_number"].astype("int32")
     df = GeoDataFrame(df)
-    df.geometry.crs = None
+    df.geometry.array.crs = None
 
     # Read the expected data
     if geometry_encoding == "WKB":
@@ -349,7 +348,7 @@ def test_geoarrow_export_with_extension_types(geometry_type, dim):
     df["geometry"] = GeoSeries.from_wkb(df["geometry"])
     df["row_number"] = df["row_number"].astype("int32")
     df = GeoDataFrame(df)
-    df.geometry.crs = None
+    df.geometry.array.crs = None
 
     pytest.importorskip("geoarrow.pyarrow")
 
diff --git a/geopandas/io/tests/test_sql.py b/geopandas/io/tests/test_sql.py
index 10712d61..524e7e49 100644
--- a/geopandas/io/tests/test_sql.py
+++ b/geopandas/io/tests/test_sql.py
@@ -490,7 +490,7 @@ class TestIO:
         table = "nybb"
 
         # Write to db
-        df_nybb.crs = None
+        df_nybb.geometry.array.crs = None
         with pytest.warns(UserWarning, match="Could not parse CRS from the GeoDataF"):
             write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
         # Validate that srid is -1
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 19c93d12..d9df3fa1 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -665,7 +665,9 @@ def test_unary_union():
         shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1)]),
     ]
     G = from_shapely(geoms)
-    with pytest.warns(FutureWarning, match="The 'unary_union' attribute is deprecated"):
+    with pytest.warns(
+        DeprecationWarning, match="The 'unary_union' attribute is deprecated"
+    ):
         u = G.unary_union()
 
     expected = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index f637b376..476fed14 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -213,7 +213,7 @@ class TestGeometryArrayCRS:
         assert s.values.crs == self.osgb
 
         # manually change CRS
-        s.crs = 4326
+        s = s.set_crs(4326, allow_override=True)
         assert s.crs == self.wgs
         assert s.values.crs == self.wgs
 
@@ -262,7 +262,7 @@ class TestGeometryArrayCRS:
         arr = from_shapely(self.geoms)
         s = GeoSeries(arr, crs=27700)
         df = GeoDataFrame(geometry=s)
-        df.crs = 4326
+        df = df.set_crs(crs="epsg:4326", allow_override=True)
         assert df.crs == self.wgs
         assert df.geometry.crs == self.wgs
         assert df.geometry.values.crs == self.wgs
@@ -418,7 +418,7 @@ class TestGeometryArrayCRS:
             FutureWarning, match="You are adding a column named 'geometry'"
         ):
             df["geometry"] = scalar
-        df.crs = 4326
+        df = df.set_crs(4326)
         assert df.crs == self.wgs
         assert df.geometry.crs == self.wgs
         assert df.geometry.values.crs == self.wgs
@@ -741,6 +741,7 @@ class TestSetCRS:
         assert non_naive.crs == "EPSG:3857"
         assert result.crs == "EPSG:3857"
 
-        # raise error when no crs is passed
-        with pytest.raises(ValueError):
-            naive.set_crs(crs=None, epsg=None)
+        # set CRS to None
+        result = non_naive.set_crs(crs=None, allow_override=True)
+        assert result.crs is None
+        assert non_naive.crs == "EPSG:3857"
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index 6cf16c35..de0ae17b 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -70,7 +70,7 @@ def test_dissolve_retains_existing_crs(nybb_polydf):
 
 
 def test_dissolve_retains_nonexisting_crs(nybb_polydf):
-    nybb_polydf.crs = None
+    nybb_polydf.geometry.array.crs = None
     test = nybb_polydf.dissolve("manhattan_bronx")
     assert test.crs is None
 
diff --git a/geopandas/tests/test_explore.py b/geopandas/tests/test_explore.py
index 2a872249..f1b15d65 100644
--- a/geopandas/tests/test_explore.py
+++ b/geopandas/tests/test_explore.py
@@ -378,7 +378,7 @@ class TestExplore:
     def test_no_crs(self):
         """Naive geometry get no tiles"""
         df = self.world.copy()
-        df.crs = None
+        df.geometry.array.crs = None
         m = df.explore()
         assert "openstreetmap" not in m.to_dict()["children"].keys()
 
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 475cf3db..8efb81ab 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -373,6 +373,20 @@ class TestDataFrame:
         with pytest.raises(AttributeError, match=msg_geo_col_missing):
             df.geometry
 
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="Requires pyproj")
+    def test_override_existing_crs_warning(self):
+        with pytest.warns(
+            DeprecationWarning,
+            match="Overriding the CRS of a GeoSeries that already has CRS",
+        ):
+            self.df.geometry.crs = "epsg:2100"
+
+        with pytest.warns(
+            DeprecationWarning,
+            match="Overriding the CRS of a GeoDataFrame that already has CRS",
+        ):
+            self.df.crs = "epsg:4326"
+
     def test_active_geometry_name(self):
         # default single active called "geometry"
         assert self.df.active_geometry_name == "geometry"
@@ -401,8 +415,7 @@ class TestDataFrame:
 
         if compat.HAS_PYPROJ:
             # assert crs is / is not preserved on mixed dataframes
-            df_nocrs = df.copy()
-            df_nocrs.crs = None
+            df_nocrs = df.copy().set_crs(None, allow_override=True)
             res1, res2 = df.align(df_nocrs)
             assert_geodataframe_equal(res1, df)
             assert res1.crs is not None
@@ -429,10 +442,8 @@ class TestDataFrame:
         assert_geodataframe_equal(res2, exp2)
 
         if compat.HAS_PYPROJ:
-            df2_nocrs = df2.copy()
-            df2_nocrs.crs = None
-            exp2_nocrs = exp2.copy()
-            exp2_nocrs.crs = None
+            df2_nocrs = df2.copy().set_crs(None, allow_override=True)
+            exp2_nocrs = exp2.copy().set_crs(None, allow_override=True)
             res1, res2 = df1.align(df2_nocrs)
             assert_geodataframe_equal(res1, exp1)
             assert res1.crs is not None
@@ -466,7 +477,7 @@ class TestDataFrame:
         assert coord == [970217.0223999023, 145643.33221435547]
 
     def test_to_json_no_crs(self):
-        self.df.crs = None
+        self.df.geometry.array.crs = None
         with pytest.raises(ValueError, match="CRS is not set"):
             self.df.to_json(to_wgs84=True)
 
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 38a943f3..dc428a53 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -318,14 +318,12 @@ class TestGeomMethods:
             result = getattr(gdf, op)
         fcmp(result, expected)
 
-    # TODO re-enable for all operations once we use pyproj > 2
-    # def test_crs_warning(self):
-    #     # operations on geometries should warn for different CRS
-    #     no_crs_g3 = self.g3.copy()
-    #     no_crs_g3.crs = None
-    #     with pytest.warns(UserWarning):
-    #         self._test_binary_topological('intersection', self.g3,
-    #                                       self.g3, no_crs_g3)
+    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    def test_crs_warning(self):
+        # operations on geometries should warn for different CRS
+        no_crs_g3 = self.g3.copy().set_crs(None, allow_override=True)
+        with pytest.warns(UserWarning):
+            self._test_binary_topological("intersection", self.g3, self.g3, no_crs_g3)
 
     def test_alignment_warning(self):
         with pytest.warns(
@@ -508,7 +506,7 @@ class TestGeomMethods:
         p2 = Polygon([(2, 0), (3, 0), (3, 1)])
         g = GeoSeries([p1, p2])
         with pytest.warns(
-            FutureWarning, match="The 'unary_union' attribute is deprecated"
+            DeprecationWarning, match="The 'unary_union' attribute is deprecated"
         ):
             result = g.unary_union
         assert result == g.union_all()
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 6f809c73..39590176 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -21,7 +21,6 @@ from shapely.geometry.base import BaseGeometry
 
 import geopandas._compat as compat
 from geopandas import GeoDataFrame, GeoSeries, clip, read_file
-from geopandas._compat import HAS_PYPROJ
 from geopandas.array import GeometryArray, GeometryDtype
 
 import pytest
@@ -39,8 +38,7 @@ class TestSeries:
         self.sq = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
         self.g1 = GeoSeries([self.t1, self.sq])
         self.g2 = GeoSeries([self.sq, self.t1])
-        self.g3 = GeoSeries([self.t1, self.t2])
-        self.g3.crs = "epsg:4326"
+        self.g3 = GeoSeries([self.t1, self.t2], crs="epsg:4326")
         self.g4 = GeoSeries([self.t2, self.t1])
         self.na = GeoSeries([self.t1, self.t2, Polygon()])
         self.na_none = GeoSeries([self.t1, self.t2, None])
@@ -83,19 +81,16 @@ class TestSeries:
         assert a1["B"].equals(a2["B"])
         assert a1["C"] is None
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_align_crs(self):
-        a1 = self.a1
-        a1.crs = "epsg:4326"
-        a2 = self.a2
-        a2.crs = "epsg:31370"
+        a1 = self.a1.set_crs("epsg:4326")
+        a2 = self.a2.set_crs("epsg:31370")
 
         res1, res2 = a1.align(a2)
         assert res1.crs == "epsg:4326"
         assert res2.crs == "epsg:31370"
 
-        a2.crs = None
-        res1, res2 = a1.align(a2)
+        res1, res2 = a1.align(a2.set_crs(None, allow_override=True))
         assert res1.crs == "epsg:4326"
         assert res2.crs is None
 
@@ -265,7 +260,7 @@ class TestSeries:
         assert np.all(self.g3.contains(self.g3.representative_point()))
         assert np.all(self.g4.contains(self.g4.representative_point()))
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_transform(self):
         utm18n = self.landmarks.to_crs(epsg=26918)
         lonlat = utm18n.to_crs(epsg=4326)
@@ -286,14 +281,14 @@ class TestSeries:
             "EPSG:32618"
         )
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_estimate_utm_crs__out_of_bounds(self):
         with pytest.raises(RuntimeError, match="Unable to determine UTM CRS"):
             GeoSeries(
                 [Polygon([(0, 90), (1, 90), (2, 90)])], crs="EPSG:4326"
             ).estimate_utm_crs()
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_estimate_utm_crs__missing_crs(self):
         with pytest.raises(RuntimeError, match="crs must be set"):
             GeoSeries([Polygon([(0, 90), (1, 90), (2, 90)])]).estimate_utm_crs()
@@ -329,7 +324,7 @@ class TestSeries:
         assert self.g1.__geo_interface__["type"] == "FeatureCollection"
         assert len(self.g1.__geo_interface__["features"]) == self.g1.shape[0]
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_proj4strings(self):
         # As string
         reprojected = self.g3.to_crs("+proj=utm +zone=30")
@@ -342,8 +337,7 @@ class TestSeries:
         assert geom_almost_equals(self.g3, reprojected_back)
 
         # Set to equivalent string, convert, compare to original
-        copy = self.g3.copy()
-        copy.crs = "epsg:4326"
+        copy = self.g3.copy().set_crs("epsg:4326", allow_override=True)
         reprojected = copy.to_crs({"proj": "utm", "zone": "30"})
         reprojected_back = reprojected.to_crs(epsg=4326)
         assert geom_almost_equals(self.g3, reprojected_back)
@@ -503,7 +497,7 @@ class TestSeries:
         expected = GeoSeries([Point(0, 2, -1), Point(3, 5, 4)])
         assert_geoseries_equal(expected, GeoSeries.from_xy(x, y, z))
 
-    @pytest.mark.skipif(HAS_PYPROJ, reason="pyproj installed")
+    @pytest.mark.skipif(compat.HAS_PYPROJ, reason="pyproj installed")
     def test_set_crs_pyproj_error(self):
         with pytest.raises(
             ImportError, match="The 'pyproj' package is required for set_crs"
@@ -541,13 +535,22 @@ def test_isna_empty_geoseries():
     assert_series_equal(result, pd.Series([], dtype="bool"))
 
 
-@pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+@pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
 def test_geoseries_crs():
-    gs = GeoSeries()
-    gs.crs = "IGNF:ETRS89UTM28"
+    gs = GeoSeries().set_crs("IGNF:ETRS89UTM28")
     assert gs.crs.to_authority() == ("IGNF", "ETRS89UTM28")
 
 
+@pytest.mark.skipif(not compat.HAS_PYPROJ, reason="Requires pyproj")
+def test_geoseries_override_existing_crs_warning():
+    gs = GeoSeries(crs="epsg:4326")
+    with pytest.warns(
+        DeprecationWarning,
+        match="Overriding the CRS of a GeoSeries that already has CRS",
+    ):
+        gs.crs = "epsg:2100"
+
+
 # -----------------------------------------------------------------------------
 # # Constructor tests
 # -----------------------------------------------------------------------------
@@ -658,7 +661,7 @@ class TestConstructor:
         assert s.name == g.name
         assert s.index is g.index
 
-    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not available")
+    @pytest.mark.skipif(not compat.HAS_PYPROJ, reason="pyproj not available")
     def test_from_series_no_set_crs_on_construction(self):
         # https://github.com/geopandas/geopandas/issues/2492
         # also when passing Series[geometry], ensure we don't change crs of
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index 3defe9cc..d810e22a 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -81,7 +81,7 @@ def test_overlay(dfs_index, how):
         expected = read_file(
             os.path.join(DATA, "polys", "df1_df2-{0}.geojson".format(name))
         )
-        expected.crs = None
+        expected.geometry.array.crs = None
         for col in expected.columns[expected.dtypes == "int32"]:
             expected[col] = expected[col].astype("int64")
         return expected
diff --git a/geopandas/tests/test_testing.py b/geopandas/tests/test_testing.py
index 6ff207cf..cb9d0b68 100644
--- a/geopandas/tests/test_testing.py
+++ b/geopandas/tests/test_testing.py
@@ -47,9 +47,9 @@ df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1})
 df2 = GeoDataFrame({"col1": [1, 2], "geometry": s2})
 
 s4 = s1.copy()
-s4.crs = 4326
+s4.array.crs = 4326
 s5 = s2.copy()
-s5.crs = 27700
+s5.array.crs = 27700
 
 s6 = GeoSeries(
     [
@@ -104,7 +104,7 @@ def test_geodataframe():
 
     assert_geodataframe_equal(df5, df4, check_like=True)
     if HAS_PYPROJ:
-        df5.geom2.crs = 3857
+        df5["geom2"] = df5.geom2.set_crs(3857, allow_override=True)
         with pytest.raises(AssertionError):
             assert_geodataframe_equal(df5, df4, check_like=True)