Repository: geopandas/geopandas
Pull request #3231 (2024-03-26 11:59:17)

ENH: add attribute requirement with spatial join
Addresses #3049 

Proposed approach to do a spatial join when both geometry and an attribute column are equal. Thought I'd get some eyes on it before I progress further as I'm new to contributing to this repo.

Things left to do:
- more tests
- changelog
- don't add suffix to the `sharedAttribute`
- consider if sharedAttribute should allow more than one attribute column, i.e. accept list and force all attributes in listed columns to be equal
- more details in docstring

martinfleis: Thanks! We should probably get #2353 in before touching sjoin code elsewhere. 

One minor nit - please use snake_case in the code, not camelCase. Thanks!
nicholas-ys-tan: Thanks Martin,

I noted #2353 has been merged. I've since re-based and resolved conflicts.

I've made further updates so it can take lists and tuples too - similar to the ``merge()`` function in pandas. Renamed ``shared_attribute`` (corrected to snake case) to ``on_attribute`` - to be a bit more similar to kwarg `on` used in pandas. Happy to change to any argument you think is preferable. Have also updated to the "_left", "_right" suffix does not get appended to the ``on_attribute`` columns. 

nicholas-ys-tan: Question on desired behaviour: 
If we are doing a `dissolve` on a `Polygon` and a `Point` with an aggfunc to be `area_averaged_mean`, would we:

- consider this invalid and produce an error message or,
- exclude the `Point` data and just do an `area_averaged_mean` of the `Polygon` data
nicholas-ys-tan: > Can you also update the docstring of GeoDataFrame.sjoin in geodataframe.py?
> 
> Otherwise this looks good to my eyes. Thanks!

I've added to the docstring, I also noted the `distance` kwarg wasn't added to the `GeoDataFrame.sjoin` docstring, I've added that here (though wasn't sure if it really should be a separate commit)
martinfleis: @nicholas-ys-tan can you merge main here to resolve conflicts?
m-richards: > For `sjoin` this looks good! I was wondering what the best approach would be, though, but I can imagine that this will depend on the characteristics of your data. But in general you could either first evaluate the spatial predicate and then the attribute match (as you did here), or either first perform the attribute merge (as a way to reduce the set of geometries for which to evaluate the spatial predicate?) and then the spatial predicate. Before I looked at the code in this PR, I was assuming this PR was for the latter, but of course for `sjoin` this should give identical results. So for a first version this is probably just fine. But I do think it might be interesting to explore performance characteristics of both on some example datasets in the future.

I will confess I was not engaging with this PR and the original issue precisely because of this, I would think the best order of operations depends quite a bit of the data in question, and for that reason was wondering if this actually belongs in geopandas. But I can still see value in providing a convenience utility for this for the general case, even if it's not necessarily the most performant in all circumstances.
nicholas-ys-tan: > For `sjoin` this looks good! I was wondering what the best approach would be, though, but I can imagine that this will depend on the characteristics of your data. But in general you could either first evaluate the spatial predicate and then the attribute match (as you did here), or either first perform the attribute merge (as a way to reduce the set of geometries for which to evaluate the spatial predicate?) and then the spatial predicate. Before I looked at the code in this PR, I was assuming this PR was for the latter, but of course for `sjoin` this should give identical results. So for a first version this is probably just fine. But I do think it might be interesting to explore performance characteristics of both on some example datasets in the future.
> 
> However, that also made me wonder if this approach is correct for the `sjoin_nearest`? Because in that case it _does_ matter in which order you do those operations, I think? I would assume that I get "the closest geometry among those with a matching attribute", but in practice it will give "the overall closest geometry _if_ that closest geometry has a matching attribute". That seems an important difference in expected behaviour, which I am not sure has been discussed?
> 
> I would maybe suggest to focus this first PR on the `sjoin` function (where the behaviour is clearer), and leave out `sjoin_nearest` for later.

Thank you @jorisvandenbossche , that's great food for thought re performance and I will do some investigation into that.

Also thank you for pointing out the implications on `sjoin_nearest`, the sequence and its impact on the output was not something that had crossed my mind. I've since removed all `on_attribute` joins from `sjoin_nearest`.

I will open up a new issue for discussion on the ordering of operations how `on_attribute` behaves in `sjoin_nearest`.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9375b628..e64df83e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,14 +1,14 @@
 # Changelog
 
-## Version 1.0.0
+## Version 1.0.0-rc1 (June 17, 2024)
 
 Notes on dependencies:
 
 - GeoPandas 1.0 drops support for shapely<2 and PyGEOS. The only geometry engine that is
   currently supported is shapely >= 2. As a consequence, spatial indexing based on the
-  rtree package has also been removed. (#3035)
+  rtree package has also been removed (#3035).
 - The I/O engine now defaults to Pyogrio which is now installed with GeoPandas instead
-  of Fiona. (#3223)
+  of Fiona (#3223).
 
 New methods:
 
@@ -36,6 +36,9 @@ New methods:
 - Added `polygonize` method exposing both `polygonize` and `polygonize_full` from
   shapely to GeoSeries/GeoDataframe (#2963).
 - Added `is_valid_reason` method from shapely to GeoSeries/GeoDataframe (#3176).
+- Added `to_arrow` method and `from_arrow` class method to
+  GeoSeries/GeoDataFrame to export and import to/from Arrow data with GeoArrow
+  extension types (#3219, #3301).
 
 New features and improvements:
 
@@ -53,18 +56,30 @@ New features and improvements:
 - 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)
+  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)
+  `GeoSeries.to_json` (#3226).
 - `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 ``on_attribute`` option argument to the ``sjoin()``
+  method, allowing to restrict joins to the observations with 
+  matching attributes. (#3231)
 - 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)
+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)
+  will silence the warning about mismatched indices (#3212).
 - `GeoSeries.set_crs` can now be used to remove CRS information by passing
-  `crs=None, allow_override=True`. (#3316)
+  `crs=None, allow_override=True` (#3316).
+- Added ``autolim`` keyword argument to ``GeoSeries.plot()`` and ``GeoDataFrame.plot()`` (#2817).
+- Added `metadata` parameter to `GeoDataFrame.to_file` (#2850)
+- Updated documentation to clarify that passing a named (Geo)Series as the `geometry`
+  argument to the GeoDataFrame constructor will not use the name but will always  
+  produce a GeoDataFrame with an active geometry column named "geometry" (#3337).
+- `read_postgis` will query the spatial_ref_sys table to determine the CRS authority
+  instead of its current behaviour of assuming EPSG. In the event the spiatal_ref_sys
+  table is not present, or the SRID is not present, `read_postgis` will fallback
+  on assuming EPSG CRS authority. (#3329)
 
 Backwards incompatible API changes:
 
@@ -80,7 +95,12 @@ Backwards incompatible API changes:
   the previous active geometry column name. This means that if the new and old names are
   different, then both columns will be preserved in the GeoDataFrame. To replicate the previous
   behaviour, you can instead call `gdf.set_geometry(ser.rename(gdf.active_geometry_name))` (#3237).
-- `delaunay_triangles` now considers all geometries together when creating the Delaunay trianguation
+  Note that this behaviour change does not affect the `GeoDataframe` constructor, passing a named 
+  GeoSeries `ser` to `GeoDataFrame(df, geometry=ser)` will always produce a GeoDataFrame with a
+  geometry column named "geometry" to preserve backwards compatibility. If you would like to
+  instead propagate the name of `ser` when constructing a GeoDataFrame, you can instead call
+  `df.set_geometry(ser)` or `GeoDataFrame(df, geometry=ser).rename_geometry(ser.name)` (#3337).
+- `delaunay_triangles` now considers all geometries together when creating the Delaunay triangulation
   instead of performing the operation element-wise. If you want to generate Delaunay
   triangles for each geometry separately, use ``shapely.delaunay_triangles`` instead. (#3273)
 - Reading a data source that does not have a geometry field using ``read_file``
@@ -109,6 +129,8 @@ Enforced deprecations:
     method directly instead.
   - Removed deprecated GeoSeries/GeoDataFrame `.plot` parameters `axes` and `colormap`, instead use
     `ax` and `cmap` respectively.
+  - Removed compatibility for specifying the `version` keyword in `to_parquet` and `to_feather`.
+    This keyword will now be passed through to pyarrow and use `schema_version` to specify the GeoParquet specification version (#3334).
 
 New deprecations:
 
@@ -130,7 +152,7 @@ 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
+- Manual overriding of an existing CRS of a GeoSeries or GeoDataFrame by setting the `crs` property has been deprecated
   and will be disabled in future. Use the `set_crs()` method instead (#3085).
 
 Bug fixes:
@@ -146,6 +168,8 @@ Bug fixes:
   pyogrio 0.8.1.
 - Fix `to_parquet` to write correct metadata in case of 3D geometries (#2824).
 - Fixes for compatibility with psycopg (#3167).
+- Fix to allow appending dataframes with no CRS to PostGIS tables with no CRS (#3328)
+- Fix plotting of all-empty GeoSeries using `explore` (#3316).
 
 ## Version 0.14.4 (April 26, 2024)
 
diff --git a/doc/source/about/roadmap.md b/doc/source/about/roadmap.md
index 298a926a..0921f074 100644
--- a/doc/source/about/roadmap.md
+++ b/doc/source/about/roadmap.md
@@ -6,62 +6,10 @@ specified date, and some may not happen at all if the implementation proves to b
 against the will of the community or face technical issues preventing their inclusion in
 the code base.
 
-The current roadmap is divided into two milestones. The first milestone aims at a
-release of the first major version of GeoPandas, while the second milestone is a
-longer-term vision covering enhancements that should happen in subsequent releases.
+The current roadmap reflects longer-term vision covering enhancements that should happen
+in upcoming releases.
 
-## Roadmap for GeoPandas 1.0
-
-### Fully vectorized geometry engine
-
-GeoPandas uses `shapely` as its geometry engine, based on scalar geometries, requiring a
-loop-based implementation of most GeoPandas methods. That comes at a significant
-performance cost, which is being resolved in shapely 2.0, a new major release resulting
-from a complete rewrite of the internals using the vectorized implementation prototyped
-in the `PyGEOS` project. At this moment, GeoPandas supports `shapely<2.0`,
-`shapely>=2.0`, and `PyGEOS` as possible geometry engines, which causes friction in the
-development process and uneven performance on the user side based on what geometry
-engine the user happens to be using.
-
-GeoPandas 1.0 will require `shapely>=2.0` and deprecate both older shapely and `PyGEOS`
-engines. This change should simplify the code base allowing more manageable maintenance
-and a lower barrier to entry for new contributors.
-
-### Feature parity with shapely
-
-Even though GeoPandas uses shapely as the geometry engine, not all its functions are
-exposed at a GeoPandas level. This has resulted in a less convenient API and a need to
-switch between `GeoSeries` objects and lists or arrays of geometries, potentially
-risking the data loss or corruption as the CRS is not included in such operations. In
-the first phase, all element-wise operations (e.g. `segmentize`, or
-`minimum_bounding_circle`) should be exposed as `GeoSeries` methods. The feature parity
-should be reached in the second phase, covering all relevant functions.
-
-### Clarity of the API
-
-The first version of the GeoPandas API is nearly ten years old. The PyData ecosystem has
-significantly changed in the meantime, and some of the early decisions may no longer be
-future-proof. Ahead of GeoPandas 1.0, the API will be revised to ensure that all the
-necessary deprecations occur before the major release to provide the stability of the
-API for the coming years.
-
-### Pruned dependencies
-
-GeoPandas offers functionality for every step of a typical geospatial workflow, from
-reading of the GIS file formats to geometry operations and handling of Coordinate
-Reference Systems (CRS) and transformation of geometries between them. However, GIS I/O
-depends on a relatively heavy C++ library `GDAL` and CRS management on another C++
-library `PROJ`, even though not every application based on GeoPandas is necessarily
-geospatial. GeoPandas 1.0 should eliminate the hard dependency on both `GDAL` and `PROJ`
-and offer the basic capability of a GeoDataFrame with a minimal set of dependencies
-limited to `pandas` and `shapely`.
-
-## Beyond GeoPandas 1.0
-
-Additional work is planned for a longer time frame, stretching beyond GeoPandas 1.0
-without a specific target release.
-
-### S2 geometry engine
+## S2 geometry engine
 
 The geometry engine used in GeoPandas is `shapely`, which serves as a Python API for
 `GEOS`. It means that all geometry operations in GeoPandas are planar, using (possibly)
@@ -74,16 +22,7 @@ The GeoPandas community is currently working together with the R-spatial communi
 has already exposed `S2` in an R counterpart of GeoPandas `sf` on Python bindings for
 `S2`, that should be used as a secondary geometry engine in GeoPandas.
 
-### Lighter-weight geospatial I/O
-
-In order to support lighter-weight installations of GeoPandas that do not depend on
-heavier and difficult to install libraries such as GDAL, additional I/O libraries should
-be developed and integrated into GeoPandas as optional dependencies.  These should be
-simpler to install and not require binary dependencies, which would lower the barrier to
-entry for GeoPandas users that need basic I/O support for a limited number of GIS
-formats such as ESRI Shapefiles or GeoPackages.
-
-### Prepared geometries
+## Prepared geometries
 
 GeoPandas is using spatial indexing for the operations that may benefit from it. Further
 performance gains can be achieved using prepared geometries. Preparation creates a
@@ -93,7 +32,7 @@ become less computationally expensive in `shapely` 2.0, GeoPandas should expose
 preparation to the user but, more importantly, use smart automatic geometry preparation
 under the hood.
 
-### Static plotting improvements
+## Static plotting improvements
 
 GeoPandas currently covers a broad range of geospatial tasks, from data exploration to
 advanced analysis. However, one moment may tempt the user to use different software -
diff --git a/doc/source/docs/reference/geodataframe.rst b/doc/source/docs/reference/geodataframe.rst
index 8a7c33d2..63843f95 100644
--- a/doc/source/docs/reference/geodataframe.rst
+++ b/doc/source/docs/reference/geodataframe.rst
@@ -3,8 +3,8 @@ GeoDataFrame
 ============
 .. currentmodule:: geopandas
 
-A ``GeoDataFrame`` is a tabular data structure that contains a column
-which contains a ``GeoSeries`` storing geometry.
+A ``GeoDataFrame`` is a tabular data structure that contains at least
+one ``GeoSeries`` column storing geometry.
 
 Constructor
 -----------
diff --git a/doc/source/docs/reference/tools.rst b/doc/source/docs/reference/tools.rst
index 618fe99b..80f539c0 100644
--- a/doc/source/docs/reference/tools.rst
+++ b/doc/source/docs/reference/tools.rst
@@ -13,6 +13,4 @@ Tools
    tools.geocode
    tools.reverse_geocode
    tools.collect
-   points_from_xy
-   datasets.available
-   datasets.get_path
+   points_from_xy
\ No newline at end of file
diff --git a/doc/source/docs/user_guide.rst b/doc/source/docs/user_guide.rst
index f68c8f8c..f593e751 100644
--- a/doc/source/docs/user_guide.rst
+++ b/doc/source/docs/user_guide.rst
@@ -21,4 +21,5 @@ Advanced topics can be found in the :doc:`Advanced Guide <advanced_guide>` and f
   Aggregation with dissolve <user_guide/aggregation_with_dissolve>
   Merging data <user_guide/mergingdata>
   Geocoding <user_guide/geocoding>
-  Sampling points <user_guide/sampling>
\ No newline at end of file
+  Sampling points <user_guide/sampling>
+  How to... <user_guide/how_to>
\ No newline at end of file
diff --git a/doc/source/docs/user_guide/data_structures.rst b/doc/source/docs/user_guide/data_structures.rst
index d8435cf8..c4c99ec9 100644
--- a/doc/source/docs/user_guide/data_structures.rst
+++ b/doc/source/docs/user_guide/data_structures.rst
@@ -25,7 +25,8 @@ of only one shape (like a single polygon) or multiple shapes that are
 meant to be thought of as one observation (like the many polygons that
 make up the State of Hawaii or a country like Indonesia).
 
-GeoPandas has three basic classes of geometric objects (which are actually *shapely* objects):
+GeoPandas has three basic classes of geometric objects (which are actually
+`Shapely <https://shapely.readthedocs.io/en/stable/manual.html>`__ objects):
 
 * Points / Multi-Points
 * Lines / Multi-Lines
@@ -82,9 +83,14 @@ GeoDataFrame
 
 A :class:`GeoDataFrame` is a tabular data structure that contains a :class:`GeoSeries`.
 
-The most important property of a :class:`GeoDataFrame` is that it always has one :class:`GeoSeries` column that holds a special status. This :class:`GeoSeries` is referred to as the :class:`GeoDataFrame`'s "geometry". When a spatial method is applied to a :class:`GeoDataFrame` (or a spatial attribute like ``area`` is called), this commands will always act on the "geometry" column.
+The most important property of a :class:`GeoDataFrame` is that it always has one :class:`GeoSeries` column that
+holds a special status - the "active geometry column". When a spatial method is applied to a
+:class:`GeoDataFrame` (or a spatial attribute like ``area`` is called), these operations will always act on the
+active geometry column.
 
-The "geometry" column -- no matter its name -- can be accessed through the :attr:`~GeoDataFrame.geometry` attribute (``gdf.geometry``), and the name of the ``geometry`` column can be found by typing ``gdf.geometry.name``.
+The active geometry column -- no matter the name of the corresponding :class:`GeoSeries` --
+can be accessed through the :attr:`~GeoDataFrame.geometry` attribute (``gdf.geometry``),
+and the name of the ``geometry`` column can be found by typing ``gdf.geometry.name`` or ``gdf.active_geometry_name``.
 
 A :class:`GeoDataFrame` may also contain other columns with geometrical (shapely) objects, but only one column can be the active geometry at a time. To change which column is the active geometry column, use the :meth:`GeoDataFrame.set_geometry` method.
 
@@ -136,9 +142,15 @@ Now, you create centroids and make it the geometry:
 Attributes and methods
 ~~~~~~~~~~~~~~~~~~~~~~
 
-Any of the attributes calls or methods described for a :class:`GeoSeries` will work on a :class:`GeoDataFrame` -- effectively, they are just applied to the "geometry" :class:`GeoSeries`.
+Any of the attributes calls or methods described for a :class:`GeoSeries` will work on a :class:`GeoDataFrame` -- they are just applied to the active geometry column :class:`GeoSeries`.
+
+However, :class:`GeoDataFrames <GeoDataFrame>` also have a number few extra methods for:
+
+* :doc:`Reading and writing files <io>`
+* :ref:`Spatial joins <mergingdata.spatial-joins>`
+* :doc:`Spatial aggregations <aggregation_with_dissolve>`
+* :doc:`Geocoding <geocoding>`
 
-However, :class:`GeoDataFrames <GeoDataFrame>` also have a few extra methods for input and output which are described on the :doc:`Reading and writing files <io>` page and for geocoding with are described in :doc:`Geocoding <geocoding>`.
 
 
 .. ipython:: python
@@ -150,8 +162,7 @@ However, :class:`GeoDataFrames <GeoDataFrame>` also have a few extra methods for
 Display options
 ---------------
 
-GeoPandas has an ``options`` attribute with currently a single configuration
-option to control:
+GeoPandas has an ``options`` attribute with global configuration attributes:
 
 .. ipython:: python
 
diff --git a/doc/source/docs/user_guide/geometric_manipulations.rst b/doc/source/docs/user_guide/geometric_manipulations.rst
index 38447e31..b9a121eb 100644
--- a/doc/source/docs/user_guide/geometric_manipulations.rst
+++ b/doc/source/docs/user_guide/geometric_manipulations.rst
@@ -106,7 +106,7 @@ Affine transformations
 
 
 Examples of geometric manipulations
-------------------------------------
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. sourcecode:: python
 
diff --git a/doc/source/docs/user_guide/how_to.rst b/doc/source/docs/user_guide/how_to.rst
new file mode 100644
index 00000000..713eacad
--- /dev/null
+++ b/doc/source/docs/user_guide/how_to.rst
@@ -0,0 +1,30 @@
+.. _how_to:
+
+How to...
+=========
+
+Drop duplicate geometry in all situations
+-----------------------------------------
+
+Using the standard Pandas :meth:`~pandas.DataFrame.drop_duplicates` function on a geometry column can lead to some duplicate
+geometries not being dropped, in certain circumstances. When used on a geometry columnm, the Pandas function compares the
+WKB of each geometry object. This is sensitive to the orders of various components of the geometry - for example, a line
+with co-ordinates in the order left-to-right should be equal to a line with the same co-ordinates in the order right-to-left,
+but the WKB representations will be different. The same applies for the order of rings of polygons and parts in multipart
+geometries.
+
+To deal with this problem, use the :meth:`~geopandas.GeoSeries.normalize` method first to order the co-ordinates in a canonincal form,
+and then use the standard :meth:`~pandas.DataFrame.drop_duplicates` method::
+
+    gdf["geometry"] = gdf.normalize()
+    gdf.drop_duplicates()
+
+The effect of the :meth:`~geopandas.GeoSeries.normalize` method can be seen in the following example::
+
+    >>> geopandas.GeoSeries([
+    ...     shapely.LineString([(0, 0), (1, 0), (2, 0)]),
+    ...     shapely.LineString([(2, 0), (1, 0), (0, 0)]),
+    ... ]).normalize().to_wkt()
+    0    LINESTRING (0 0, 1 0, 2 0)
+    1    LINESTRING (0 0, 1 0, 2 0)
+    dtype: object
diff --git a/doc/source/docs/user_guide/io.rst b/doc/source/docs/user_guide/io.rst
index 2e659a39..35f60f87 100644
--- a/doc/source/docs/user_guide/io.rst
+++ b/doc/source/docs/user_guide/io.rst
@@ -73,7 +73,7 @@ supported by that project::
 You can also read path objects::
 
     import pathlib
-    path_object = pathlib.path(filename)
+    path_object = pathlib.Path(filename)
     df = geopandas.read_file(path_object)
 
 Reading subsets of the data
@@ -173,21 +173,23 @@ Load in a subset of data with a `SQL WHERE clause <https://gdal.org/user/ogr_sql
         where="subborough='Coney Island'",
     )
 
-Supported drivers
-~~~~~~~~~~~~~~~~~
+Supported drivers / file formats
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 When using pyogrio, all drivers supported by the GDAL installation are enabled,
 and you can check those with::
 
     import pyogrio; pyogrio.list_drivers()
 
-Fiona only exposes the default drivers. To display those, type::
+where the values indicate whether reading, writing or both are supported for
+a given driver.
+Fiona only exposes a default subset of drivers. To display those, type::
 
     import fiona; fiona.supported_drivers
 
 There is a `list of available drivers <https://github.com/Toblerity/Fiona/blob/master/fiona/drvsupport.py>`_
-which are unexposed but supported (depending on the GDAL-build). You can activate
-these on runtime by updating the `supported_drivers` dictionary like::
+which are unexposed by default but may be supported (depending on the GDAL-build). You can activate
+these at runtime by updating the `supported_drivers` dictionary like::
 
     fiona.supported_drivers["NAS"] = "raw"
 
@@ -271,10 +273,10 @@ the spatial information.
 
 .. note::
 
-    This is tracking version 1.0.0 of the GeoParquet specification at:
+    The GeoParquet specification is developed at:
     https://github.com/opengeospatial/geoparquet.
 
-    Previous versions are still supported as well. By default, the latest
-    version is used when writing files (older versions can be specified using
-    the ``schema_version`` keyword), and GeoPandas supports reading files
-    of any version.
+    By default, the latest
+    version is used when writing files, but older versions can be specified using
+    the ``schema_version`` keyword. GeoPandas supports reading files
+    encoded using any GeoParquet version.
diff --git a/doc/source/docs/user_guide/mergingdata.rst b/doc/source/docs/user_guide/mergingdata.rst
index b2be6314..b09ea432 100644
--- a/doc/source/docs/user_guide/mergingdata.rst
+++ b/doc/source/docs/user_guide/mergingdata.rst
@@ -77,6 +77,7 @@ that initially has only area ID for each geometry by merging it with a :class:`~
    chicago_shapes = chicago_shapes.merge(chicago_names, on='NID')
    chicago_shapes.head()
 
+.. _mergingdata.spatial-joins:
 
 Spatial joins
 ----------------
diff --git a/doc/source/docs/user_guide/projections.rst b/doc/source/docs/user_guide/projections.rst
index b589c10a..fbba9c15 100644
--- a/doc/source/docs/user_guide/projections.rst
+++ b/doc/source/docs/user_guide/projections.rst
@@ -252,7 +252,7 @@ One possible way to find out the EPSG code is using pyproj for this:
 value if the match is not perfect)
 
 Further, on websites such as `Spatial Reference <https://spatialreference.org/>`__
-and `epsg.io <https://epsg.io/>`__ the descriptions of many CRS can be found
+and `epsg.org <https://epsg.org/>`__ the descriptions of many CRS can be found
 including their EPSG codes and proj4 string definitions.
 
 **Other formats**
@@ -363,7 +363,7 @@ can be a reason that the CRS object is not equal to the expected EPSG code.
 
 Consider the following example of a Canadian projected CRS EPSG:2953. When
 constructing the CRS object from the WKT string as provided on
-`EPSG:2953 <https://epsg.io/2953>`__:
+`EPSG:2953 <https://spatialreference.org/ref/epsg/2953/>`__:
 
 .. code-block:: python
 
diff --git a/doc/source/gallery/plotting_basemap_background.ipynb b/doc/source/gallery/plotting_basemap_background.ipynb
index 6b0ddc12..77da1a89 100644
--- a/doc/source/gallery/plotting_basemap_background.ipynb
+++ b/doc/source/gallery/plotting_basemap_background.ipynb
@@ -57,7 +57,7 @@
     "coordinate reference systems (CRS) of the tiles and the data match.\n",
     "Web map tiles are typically provided in\n",
     "[Web Mercator](https://en.wikipedia.org/wiki/Web_Mercator>)\n",
-    "([EPSG 3857](https://epsg.io/3857)), so let us first check what\n",
+    "([EPSG 3857](https://spatialreference.org/ref/epsg/3857/)), so let us first check what\n",
     "CRS our NYC boroughs are in:"
    ]
   },
diff --git a/doc/source/gallery/polygon_plotting_with_folium.ipynb b/doc/source/gallery/polygon_plotting_with_folium.ipynb
index 1da13782..b6f7b68a 100644
--- a/doc/source/gallery/polygon_plotting_with_folium.ipynb
+++ b/doc/source/gallery/polygon_plotting_with_folium.ipynb
@@ -81,7 +81,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The [epsg:2263](https://epsg.io/2263) crs is a projected coordinate reference system with linear units (ft in this case).\n",
+    "The [EPSG:2263](https://spatialreference.org/ref/epsg/2263/) crs is a projected coordinate reference system with linear units (ft in this case).\n",
     "As folium (i.e. leaflet.js) by default accepts values of latitude and longitude (angular units) as input, we need to project the geometry to a geographic coordinate system first."
    ]
   },
diff --git a/doc/source/getting_started.md b/doc/source/getting_started.md
index e5578371..a9df9b05 100644
--- a/doc/source/getting_started.md
+++ b/doc/source/getting_started.md
@@ -29,7 +29,7 @@ conda install -c conda-forge geopandas
 
 ### Detailed instructions
 
-Do you prefer ``pip install`` or installation from source? Or specific version? See
+Do you prefer ``pip install`` or installation from source? Or want to install a specific version? See
 {doc}`detailed instructions <getting_started/install>`.
 
 ### What now?
@@ -37,7 +37,7 @@ Do you prefer ``pip install`` or installation from source? Or specific version?
 - If you don't have GeoPandas yet, check {doc}`Installation <getting_started/install>`.
 - If you have never used GeoPandas and want to get familiar with it and its core
   functionality quickly, see {doc}`Getting Started Tutorial <getting_started/introduction>`.
-- Detailed illustration how to work with different parts of GeoPandas, how to make maps,
+- Detailed illustration of how to work with different parts of GeoPandas, make maps,
   manage projections, spatially merge data or geocode are part of our
   {doc}`User Guide <docs/user_guide>`.
 - And if you are interested in the complete
diff --git a/doc/source/getting_started/install.rst b/doc/source/getting_started/install.rst
index b6b1d112..6246f998 100644
--- a/doc/source/getting_started/install.rst
+++ b/doc/source/getting_started/install.rst
@@ -16,7 +16,7 @@ To install GeoPandas and all its dependencies, we recommend to use the `conda`_
 package manager. This can be obtained by installing the
 `Anaconda Distribution`_ (a free Python distribution for data science), or
 through `miniconda`_ (minimal distribution only containing Python and the
-`conda`_ package manager). See also the `installation docs
+`conda`_ package manager). See also the `Conda installation docs
 <https://conda.io/docs/user-guide/install/download.html>`__ for more information
 on how to install Anaconda or miniconda locally.
 
@@ -38,7 +38,8 @@ which packages can be installed, in addition to the "*defaults*" channel
 provided by Anaconda.
 Depending on what other packages you are working with, the *defaults* channel
 or *conda-forge* channel may be better for your needs (e.g. some packages are
-available on *conda-forge* and not on *defaults*).
+available on *conda-forge* and not on *defaults*). Generally, *conda-forge* is
+more up to date with the latest versions of geospatial python packages.
 
 GeoPandas and all its dependencies are available on the *conda-forge*
 channel, and can be installed as::
@@ -163,6 +164,7 @@ For plotting, these additional packages may be used:
 
 - `matplotlib`_ (>= 3.5.0)
 - `mapclassify`_ (>= 2.4.0)
+- `folium`_ (for interactive plotting)
 
 
 .. _PyPI: https://pypi.python.org/pypi/geopandas
@@ -210,3 +212,5 @@ For plotting, these additional packages may be used:
 .. _packaging: https://packaging.pypa.io/en/latest/
 
 .. _pointpats: https://pysal.org/pointpats/
+
+.. _folium: https://python-visualization.github.io/folium/stable/
diff --git a/doc/source/getting_started/introduction.ipynb b/doc/source/getting_started/introduction.ipynb
index a23a0946..b50052fc 100644
--- a/doc/source/getting_started/introduction.ipynb
+++ b/doc/source/getting_started/introduction.ipynb
@@ -15,13 +15,14 @@
     "\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 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",
+    "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; unlike in some typical desktop GIS software.\n",
     "\n",
     "![geodataframe schema](../_static/dataframe.svg)\n",
     "\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",
-    "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",
+    "Only one `GeoSeries` in a `GeoDataFrame` is considered the _active geometry column_, which means that all geometric operations applied to a `GeoDataFrame` operate on this column.\n",
+    "The active geometry column is accessed via the `GeoDataFrame.geometry` attribute.\n",
     "\n",
     "\n",
     "<div class=\"alert alert-info\">\n",
@@ -65,7 +66,7 @@
    "source": [
     "### Writing files\n",
     "\n",
-    "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."
+    "To write a `GeoDataFrame` back to file use `GeoDataFrame.to_file()`. By default, the file format is inferred by the file extension (e.g. `.shp` for Shapefile, `.geojson` for GeoJSON), but you can specify this explicitly with the `driver` keyword."
    ]
   },
   {
diff --git a/geopandas/explore.py b/geopandas/explore.py
index 899be57d..bc14fb53 100644
--- a/geopandas/explore.py
+++ b/geopandas/explore.py
@@ -1,3 +1,4 @@
+import warnings
 from packaging.version import Version
 from statistics import mean
 
@@ -341,7 +342,7 @@ def _explore(
         # Get bounds to specify location and map extent
         bounds = gdf.total_bounds
         location = kwargs.pop("location", None)
-        if location is None:
+        if location is None and not np.isnan(bounds).all():
             x = mean([bounds[0], bounds[2]])
             y = mean([bounds[1], bounds[3]])
             location = (y, x)
@@ -394,6 +395,15 @@ def _explore(
         if fit:
             m.fit_bounds([[bounds[1], bounds[0]], [bounds[3], bounds[2]]])
 
+    if gdf.is_empty.all():
+        warnings.warn(
+            "The GeoSeries you are attempting to plot is "
+            "composed of empty geometries. Nothing has been displayed.",
+            UserWarning,
+            stacklevel=3,
+        )
+        return m
+
     for map_kwd in _MAP_KWARGS:
         kwargs.pop(map_kwd, None)
 
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 19e7f2ee..339315d6 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -77,8 +77,8 @@ crs_mismatch_error = (
 
 class GeoDataFrame(GeoPandasBase, DataFrame):
     """
-    A GeoDataFrame object is a pandas.DataFrame that has a column
-    with geometry. In addition to the standard DataFrame constructor arguments,
+    A GeoDataFrame object is a pandas.DataFrame that has one or more columns
+    containing geometry. In addition to the standard DataFrame constructor arguments,
     GeoDataFrame also accepts the following keyword arguments:
 
     Parameters
@@ -87,9 +87,16 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         Coordinate Reference System of the geometry objects. 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.
-    geometry : str or array (optional)
-        If str, column to use as geometry. If array, will be set as 'geometry'
-        column on GeoDataFrame.
+    geometry : str or array-like (optional)
+        Value to use as the active geometry column.
+        If str, treated as column name to use. If array-like, it will be
+        added as new column named 'geometry' on the GeoDataFrame and set as the
+        active geometry column.
+
+        Note that if ``geometry`` is a (Geo)Series with a
+        name, the name will not be used, a column named "geometry" will still be
+        added. To preserve the name, you can use :meth:`~GeoDataFrame.rename_geometry`
+        to update the geometry column name.
 
     Examples
     --------
@@ -194,6 +201,11 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
             ):
                 raise ValueError(crs_mismatch_error)
 
+            if hasattr(geometry, "name") and geometry.name not in ("geometry", None):
+                # __init__ always creates geometry col named "geometry"
+                # rename as `set_geometry` respects the given series name
+                geometry = geometry.rename("geometry")
+
             self.set_geometry(geometry, inplace=True, crs=crs)
 
         if geometry is None and crs:
@@ -816,6 +828,42 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
         return df
 
+    @classmethod
+    def from_arrow(cls, table, geometry=None):
+        """
+        Construct a GeoDataFrame from a Arrow table object based on GeoArrow
+        extension types.
+
+        See https://geoarrow.org/ for details on the GeoArrow specification.
+
+        This functions accepts any tabular Arrow object implementing
+        the `Arrow PyCapsule Protocol`_ (i.e. having an ``__arrow_c_array__``
+        or ``__arrow_c_stream__`` method).
+
+        .. _Arrow PyCapsule Protocol: https://arrow.apache.org/docs/format/CDataInterface/PyCapsuleInterface.html
+
+        .. versionadded:: 1.0
+
+        Parameters
+        ----------
+        table : pyarrow.Table or Arrow-compatible table
+            Any tabular object implementing the Arrow PyCapsule Protocol
+            (i.e. has an ``__arrow_c_array__`` or ``__arrow_c_stream__``
+            method). This table should have at least one column with a
+            geoarrow geometry type.
+        geometry : str, default None
+            The name of the geometry column to set as the active geometry
+            column. If None, the first geometry column found will be used.
+
+        Returns
+        -------
+        GeoDataFrame
+
+        """
+        from geopandas.io._geoarrow import arrow_to_geopandas
+
+        return arrow_to_geopandas(table, geometry=geometry)
+
     def to_json(
         self, na="null", show_bbox=False, drop_id=False, to_wgs84=False, **kwargs
     ):
@@ -1222,7 +1270,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
 
         >>> arrow_table = gdf.to_arrow()
         >>> arrow_table
-        <geopandas.io.geoarrow.ArrowTable object at ...>
+        <geopandas.io._geoarrow.ArrowTable object at ...>
 
         The returned data object needs to be consumed by a library implementing
         the Arrow PyCapsule Protocol. For example, wrapping the data as a
@@ -1240,7 +1288,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
 01010000000000000000000040000000000000F03F]]
 
         """
-        from geopandas.io.geoarrow import ArrowTable, geopandas_to_arrow
+        from geopandas.io._geoarrow import ArrowTable, geopandas_to_arrow
 
         table, _ = geopandas_to_arrow(
             self,
@@ -1263,7 +1311,8 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
     ):
         """Write a GeoDataFrame to the Parquet format.
 
-        Any geometry columns present are serialized to WKB format in the file.
+        By default, all geometry columns present are serialized to WKB format
+        in the file.
 
         Requires 'pyarrow'.
 
@@ -1284,14 +1333,21 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             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.
+            Note: the "geoarrow" option is part of the newer GeoParquet 1.1
+            specification, should be considered as experimental, and may not
+            be supported by all readers.
         write_covering_bbox : bool, default False
             Writes the bounding box column for each row entry with column
             name 'bbox'. Writing a bbox column can be computationally
             expensive, but allows you to specify a `bbox` in :
             func:`read_parquet` for filtered reading.
+            Note: this bbox column is part of the newer GeoParquet 1.1
+            specification and should be considered as experimental. While
+            writing the column is backwards compatible, using it for filtering
+            may not be supported by all readers.
+        schema_version : {'0.1.0', '0.4.0', '1.0.0', '1.1.0', None}
+            GeoParquet specification version; if not provided, will default to
+            latest supported stable version (1.0.0).
         kwargs
             Additional keyword arguments passed to :func:`pyarrow.parquet.write_table`.
 
@@ -1430,6 +1486,9 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             The underlying library that is used to write the file. Currently, the
             supported options are "pyogrio" and "fiona". Defaults to "pyogrio" if
             installed, otherwise tries "fiona".
+        metadata : dict[str, str], default None
+            Optional metadata to be stored in the file. Keys and values must be
+            strings. Supported only for "GPKG" driver.
         **kwargs :
             Keyword args to be passed to the engine, and can be used to write
             to multi-layer data, store data within archives (zip files), etc.
@@ -1813,7 +1872,12 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             return _geodataframe_constructor_with_fallback(
                 pd.DataFrame._from_mgr(mgr, axes)
             )
-        return GeoDataFrame._from_mgr(mgr, axes)
+        gdf = GeoDataFrame._from_mgr(mgr, axes)
+        # _from_mgr doesn't preserve metadata (expect __finalize__ to be called)
+        # still need to mimic __init__ behaviour with geometry=None
+        if (gdf.columns == "geometry").sum() == 1:  # only if "geometry" is single col
+            gdf._geometry_column_name = "geometry"
+        return gdf
 
     @property
     def _constructor_sliced(self):
@@ -2265,6 +2329,16 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
             Suffix to apply to overlapping column names (left GeoDataFrame).
         rsuffix : string, default 'right'
             Suffix to apply to overlapping column names (right GeoDataFrame).
+        distance : number or array_like, optional
+            Distance(s) around each input geometry within which to query the tree
+            for the 'dwithin' predicate. If array_like, must be
+            one-dimesional with length equal to length of left GeoDataFrame.
+            Required if ``predicate='dwithin'``.
+        on_attribute : string, list or tuple
+            Column name(s) to join on as an additional join restriction on top
+            of the spatial predicate. These must be found in both DataFrames.
+            If set, observations are joined only if the predicate applies
+            and values in specified columns match.
 
         Examples
         --------
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 505304e7..ea30f61d 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -556,6 +556,41 @@ class GeoSeries(GeoPandasBase, Series):
             **kwargs,
         )
 
+    @classmethod
+    def from_arrow(cls, arr, **kwargs) -> GeoSeries:
+        """
+        Construct a GeoSeries from a Arrow array object with a GeoArrow
+        extension type.
+
+        See https://geoarrow.org/ for details on the GeoArrow specification.
+
+        This functions accepts any Arrow array object implementing
+        the `Arrow PyCapsule Protocol`_ (i.e. having an ``__arrow_c_array__``
+        method).
+
+        .. _Arrow PyCapsule Protocol: https://arrow.apache.org/docs/format/CDataInterface/PyCapsuleInterface.html
+
+        .. versionadded:: 1.0
+
+        Parameters
+        ----------
+        arr : pyarrow.Array, Arrow array
+            Any array object implementing the Arrow PyCapsule Protocol
+            (i.e. has an ``__arrow_c_array__`` or ``__arrow_c_stream__``
+            method). The type of the array should be one of the
+            geoarrow geometry types.
+        **kwargs
+            Other parameters passed to the GeoSeries constructor.
+
+        Returns
+        -------
+        GeoSeries
+
+        """
+        from geopandas.io._geoarrow import arrow_to_geometry_array
+
+        return cls(arrow_to_geometry_array(arr), **kwargs)
+
     @property
     def __geo_interface__(self) -> Dict:
         """Returns a ``GeoSeries`` as a python feature collection.
@@ -1376,7 +1411,7 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
 
         >>> arrow_array = gser.to_arrow()
         >>> arrow_array
-        <geopandas.io.geoarrow.GeoArrowArray object at ...>
+        <geopandas.io._geoarrow.GeoArrowArray object at ...>
 
         The returned array object needs to be consumed by a library implementing
         the Arrow PyCapsule Protocol. For example, wrapping the data as a
@@ -1394,7 +1429,7 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
         """
         import pyarrow as pa
 
-        from geopandas.io.geoarrow import (
+        from geopandas.io._geoarrow import (
             GeoArrowArray,
             construct_geometry_array,
             construct_wkb_array,
diff --git a/geopandas/io/geoarrow.py b/geopandas/io/_geoarrow.py
similarity index 81%
rename from geopandas/io/geoarrow.py
rename to geopandas/io/_geoarrow.py
index 004152ca..32ccf519 100644
--- a/geopandas/io/geoarrow.py
+++ b/geopandas/io/_geoarrow.py
@@ -10,7 +10,21 @@ from numpy.typing import NDArray
 import shapely
 from shapely import GeometryType
 
+from geopandas import GeoDataFrame
 from geopandas._compat import SHAPELY_GE_204
+from geopandas.array import from_shapely, from_wkb
+
+GEOARROW_ENCODINGS = [
+    "point",
+    "linestring",
+    "polygon",
+    "multipoint",
+    "multilinestring",
+    "multipolygon",
+]
+
+
+## GeoPandas -> GeoArrow
 
 
 class ArrowTable:
@@ -407,6 +421,119 @@ def construct_geometry_array(
         raise ValueError(f"Unsupported type for geoarrow: {geom_type}")
 
 
+## GeoArrow -> GeoPandas
+
+
+def _get_arrow_geometry_field(field):
+    if (meta := field.metadata) is not None:
+        if (ext_name := meta.get(b"ARROW:extension:name", None)) is not None:
+            if ext_name.startswith(b"geoarrow."):
+                if (
+                    ext_meta := meta.get(b"ARROW:extension:metadata", None)
+                ) is not None:
+                    ext_meta = json.loads(ext_meta.decode())
+                return ext_name.decode(), ext_meta
+
+    if isinstance(field.type, pa.ExtensionType):
+        ext_name = field.type.extension_name
+        if ext_name.startswith("geoarrow."):
+            ext_meta_ser = field.type.__arrow_ext_serialize__()
+            if ext_meta_ser:
+                ext_meta = json.loads(ext_meta_ser.decode())
+            else:
+                ext_meta = None
+            return ext_name, ext_meta
+
+    return None
+
+
+def arrow_to_geopandas(table, geometry=None):
+    """
+    Convert Arrow table object to a GeoDataFrame based on GeoArrow extension types.
+
+    Parameters
+    ----------
+    table : pyarrow.Table
+        The Arrow table to convert.
+    geometry : str, default None
+        The name of the geometry column to set as the active geometry
+        column. If None, the first geometry column found will be used.
+
+    Returns
+    -------
+    GeoDataFrame
+
+    """
+    if not isinstance(table, pa.Table):
+        table = pa.table(table)
+
+    geom_fields = []
+
+    for i, field in enumerate(table.schema):
+        geom = _get_arrow_geometry_field(field)
+        if geom is not None:
+            geom_fields.append((i, field.name, *geom))
+
+    if len(geom_fields) == 0:
+        raise ValueError("No geometry column found in the Arrow table.")
+
+    table_attr = table.drop([f[1] for f in geom_fields])
+    df = table_attr.to_pandas()
+
+    for i, col, ext_name, ext_meta in geom_fields:
+        crs = None
+        if ext_meta is not None and "crs" in ext_meta:
+            crs = ext_meta["crs"]
+
+        if ext_name == "geoarrow.wkb":
+            geom_arr = from_wkb(np.array(table[col]), crs=crs)
+        elif ext_name.split(".")[1] in GEOARROW_ENCODINGS:
+
+            geom_arr = from_shapely(
+                construct_shapely_array(table[col].combine_chunks(), ext_name), crs=crs
+            )
+        else:
+            raise TypeError(f"Unknown GeoArrow extension type: {ext_name}")
+
+        df.insert(i, col, geom_arr)
+
+    return GeoDataFrame(df, geometry=geometry or geom_fields[0][1])
+
+
+def arrow_to_geometry_array(arr):
+    """
+    Convert Arrow array object (representing single GeoArrow array) to a
+    geopandas GeometryArray.
+
+    Specifically for GeoSeries.from_arrow.
+    """
+    if Version(pa.__version__) < Version("14.0.0"):
+        raise ValueError("Importing from Arrow requires pyarrow >= 14.0.")
+
+    schema_capsule, array_capsule = arr.__arrow_c_array__()
+    field = pa.Field._import_from_c_capsule(schema_capsule)
+    pa_arr = pa.Array._import_from_c_capsule(field.__arrow_c_schema__(), array_capsule)
+
+    geom_info = _get_arrow_geometry_field(field)
+    if geom_info is None:
+        raise ValueError("No GeoArrow geometry field found.")
+    ext_name, ext_meta = geom_info
+
+    crs = None
+    if ext_meta is not None and "crs" in ext_meta:
+        crs = ext_meta["crs"]
+
+    if ext_name == "geoarrow.wkb":
+        geom_arr = from_wkb(np.array(pa_arr), crs=crs)
+    elif ext_name.split(".")[1] in GEOARROW_ENCODINGS:
+
+        geom_arr = from_shapely(construct_shapely_array(pa_arr, ext_name), crs=crs)
+    else:
+        raise ValueError(f"Unknown GeoArrow extension type: {ext_name}")
+
+    return geom_arr
+
+
 def _get_inner_coords(arr):
     if pa.types.is_struct(arr.type):
         if arr.type.num_fields == 2:
@@ -433,6 +560,9 @@ def construct_shapely_array(arr: pa.Array, extension_name: str):
     with GeoArrow extension type.
 
     """
+    if isinstance(arr, pa.ExtensionArray):
+        arr = arr.storage
+
     if extension_name == "geoarrow.point":
         coords = _get_inner_coords(arr)
         result = shapely.from_ragged_array(GeometryType.POINT, coords, None)
diff --git a/geopandas/io/arrow.py b/geopandas/io/arrow.py
index c0b2f3f7..53cf77ed 100644
--- a/geopandas/io/arrow.py
+++ b/geopandas/io/arrow.py
@@ -335,7 +335,7 @@ def _geopandas_to_arrow(
     """
     from pyarrow import StructArray
 
-    from geopandas.io.geoarrow import geopandas_to_arrow
+    from geopandas.io._geoarrow import geopandas_to_arrow
 
     _validate_dataframe(df)
 
@@ -428,17 +428,6 @@ def _to_parquet(
         "pyarrow.parquet", extra="pyarrow is required for Parquet support."
     )
 
-    if kwargs and "version" in kwargs and kwargs["version"] is not None:
-        if schema_version is None and kwargs["version"] in SUPPORTED_VERSIONS:
-            warnings.warn(
-                "the `version` parameter has been replaced with `schema_version`. "
-                "`version` will instead be passed directly to the underlying "
-                "parquet writer unless `version` is 0.1.0 or 0.4.0.",
-                FutureWarning,
-                stacklevel=2,
-            )
-            schema_version = kwargs.pop("version")
-
     path = _expand_user(path)
     table = _geopandas_to_arrow(
         df,
@@ -491,17 +480,6 @@ def _to_feather(df, path, index=None, compression=None, schema_version=None, **k
     if Version(pyarrow.__version__) < Version("0.17.0"):
         raise ImportError("pyarrow >= 0.17 required for Feather support")
 
-    if kwargs and "version" in kwargs and kwargs["version"] is not None:
-        if schema_version is None and kwargs["version"] in SUPPORTED_VERSIONS:
-            warnings.warn(
-                "the `version` parameter has been replaced with `schema_version`. "
-                "`version` will instead be passed directly to the underlying "
-                "feather writer unless `version` is 0.1.0 or 0.4.0.",
-                FutureWarning,
-                stacklevel=2,
-            )
-            schema_version = kwargs.pop("version")
-
     path = _expand_user(path)
     table = _geopandas_to_arrow(df, index=index, schema_version=schema_version)
     feather.write_feather(table, path, compression=compression, **kwargs)
@@ -511,15 +489,17 @@ def _arrow_to_geopandas(table, geo_metadata=None):
     """
     Helper function with main, shared logic for read_parquet/read_feather.
     """
-    df = table.to_pandas()
-
-    geo_metadata = geo_metadata or _decode_metadata(
-        table.schema.metadata.get(b"geo", b"")
-    )
+    if geo_metadata is None:
+        # Note: this path of not passing metadata is also used by dask-geopandas
+        geo_metadata = _validate_and_decode_metadata(table.schema.metadata)
 
     # Find all geometry columns that were read from the file.  May
     # be a subset if 'columns' parameter is used.
-    geometry_columns = df.columns.intersection(geo_metadata["columns"])
+    geometry_columns = [
+        col for col in geo_metadata["columns"] if col in table.column_names
+    ]
+    result_column_names = list(table.slice(0, 0).to_pandas().columns)
+    geometry_columns.sort(key=result_column_names.index)
 
     if not len(geometry_columns):
         raise ValueError(
@@ -543,6 +523,9 @@ def _arrow_to_geopandas(table, geo_metadata=None):
                 stacklevel=3,
             )
 
+    table_attr = table.drop(geometry_columns)
+    df = table_attr.to_pandas()
+
     # Convert the WKB columns that are present back to geometry.
     for col in geometry_columns:
         col_metadata = geo_metadata["columns"][col]
@@ -556,17 +539,19 @@ def _arrow_to_geopandas(table, geo_metadata=None):
             crs = "OGC:CRS84"
 
         if col_metadata["encoding"] == "WKB":
-            df[col] = from_wkb(df[col].values, crs=crs)
+            geom_arr = from_wkb(np.array(table[col]), crs=crs)
         else:
-            from geopandas.io.geoarrow import construct_shapely_array
+            from geopandas.io._geoarrow import construct_shapely_array
 
-            df[col] = from_shapely(
+            geom_arr = from_shapely(
                 construct_shapely_array(
                     table[col].combine_chunks(), "geoarrow." + col_metadata["encoding"]
                 ),
                 crs=crs,
             )
 
+        df.insert(result_column_names.index(col), col, geom_arr)
+
     return GeoDataFrame(df, geometry=geometry)
 
 
@@ -860,7 +845,6 @@ def _read_feather(path, columns=None, **kwargs):
     path = _expand_user(path)
 
     table = feather.read_table(path, columns=columns, **kwargs)
-    _validate_and_decode_metadata(table.schema.metadata)
     return _arrow_to_geopandas(table)
 
 
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index 4d39ffdb..dae1f46b 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import os
 import urllib.request
 import warnings
@@ -26,7 +28,6 @@ _VALID_URLS.discard("")
 # file:// URIs are supported by fiona/pyogrio -> don't already open + read the file here
 _VALID_URLS.discard("file")
 
-
 fiona = None
 fiona_env = None
 fiona_import_error = None
@@ -103,6 +104,20 @@ def _check_pyogrio(func):
         )
 
 
+def _check_metadata_supported(metadata: str | None, engine: str, driver: str) -> None:
+    if metadata is None:
+        return
+    if driver != "GPKG":
+        raise NotImplementedError(
+            "The 'metadata' keyword is only supported for the GPKG driver."
+        )
+
+    if engine == "fiona" and not FIONA_GE_19:
+        raise NotImplementedError(
+            "The 'metadata' keyword is only supported for Fiona >= 1.9."
+        )
+
+
 def _check_engine(engine, func):
     # if not specified through keyword or option, then default to "pyogrio" if
     # installed, otherwise try fiona
@@ -565,6 +580,7 @@ def _to_file(
     mode="w",
     crs=None,
     engine=None,
+    metadata=None,
     **kwargs,
 ):
     """
@@ -620,6 +636,10 @@ def _to_file(
         supported options are "pyogrio" and "fiona". Defaults to "pyogrio" if
         installed, otherwise tries "fiona". Engine can also be set globally
         with the ``geopandas.options.io_engine`` option.
+    metadata : dict[str, str], default None
+        Optional metadata to be stored in the file. Keys and values must be
+        strings. Only supported for the "GPKG" driver
+        (requires Fiona >= 1.9 or pyogrio >= 0.6).
     **kwargs :
         Keyword args to be passed to the engine, and can be used to write
         to multi-layer data, store data within archives (zip files), etc.
@@ -663,19 +683,20 @@ def _to_file(
             "to a supported format like a well-known text (WKT) using "
             "`GeoSeries.to_wkt()`.",
         )
+    _check_metadata_supported(metadata, engine, driver)
 
     if mode not in ("w", "a"):
         raise ValueError(f"'mode' should be one of 'w' or 'a', got '{mode}' instead")
 
     if engine == "pyogrio":
-        _to_file_pyogrio(df, filename, driver, schema, crs, mode, **kwargs)
+        _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs)
     elif engine == "fiona":
-        _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs)
+        _to_file_fiona(df, filename, driver, schema, crs, mode, metadata, **kwargs)
     else:
         raise ValueError(f"unknown engine '{engine}'")
 
 
-def _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs):
+def _to_file_fiona(df, filename, driver, schema, crs, mode, metadata, **kwargs):
     if not HAS_PYPROJ and crs:
         raise ImportError(
             "The 'pyproj' package is required to write a file with a CRS, but it is not"
@@ -707,10 +728,12 @@ def _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs):
         with fiona.open(
             filename, mode=mode, driver=driver, crs_wkt=crs_wkt, schema=schema, **kwargs
         ) as colxn:
+            if metadata is not None:
+                colxn.update_tags(metadata)
             colxn.writerecords(df.iterfeatures())
 
 
-def _to_file_pyogrio(df, filename, driver, schema, crs, mode, **kwargs):
+def _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs):
     import pyogrio
 
     if schema is not None:
@@ -728,7 +751,7 @@ def _to_file_pyogrio(df, filename, driver, schema, crs, mode, **kwargs):
     if not df.columns.is_unique:
         raise ValueError("GeoDataFrame cannot contain duplicated column names.")
 
-    pyogrio.write_dataframe(df, filename, driver=driver, **kwargs)
+    pyogrio.write_dataframe(df, filename, driver=driver, metadata=metadata, **kwargs)
 
 
 def infer_schema(df):
diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py
index aa3c9f8e..0f99b09e 100644
--- a/geopandas/io/sql.py
+++ b/geopandas/io/sql.py
@@ -1,5 +1,6 @@
 import warnings
 from contextlib import contextmanager
+from functools import lru_cache
 
 import pandas as pd
 
@@ -41,7 +42,7 @@ def _get_conn(conn_or_engine):
         raise ValueError(f"Unknown Connectable: {conn_or_engine}")
 
 
-def _df_to_geodf(df, geom_col="geom", crs=None):
+def _df_to_geodf(df, geom_col="geom", crs=None, con=None):
     """
     Transforms a pandas DataFrame into a GeoDataFrame.
     The column 'geom_col' must be a geometry column in WKB representation.
@@ -58,6 +59,8 @@ def _df_to_geodf(df, geom_col="geom", crs=None):
         such as an authority string (eg "EPSG:4326") or a WKT string.
         If not set, tries to determine CRS from the SRID associated with the
         first geometry in the database, and assigns that to all geometries.
+    con : sqlalchemy.engine.Connection or sqlalchemy.engine.Engine
+        Active connection to the database to query.
     Returns
     -------
     GeoDataFrame
@@ -92,7 +95,28 @@ def _df_to_geodf(df, geom_col="geom", crs=None):
             srid = shapely.get_srid(geoms.iat[0])
             # if no defined SRID in geodatabase, returns SRID of 0
             if srid != 0:
-                crs = "epsg:{}".format(srid)
+                try:
+                    spatial_ref_sys_df = _get_spatial_ref_sys_df(con, srid)
+                except pd.errors.DatabaseError:
+                    warning_msg = (
+                        f"Could not find the spatial reference system table "
+                        f"(spatial_ref_sys) in PostGIS."
+                        f"Trying epsg:{srid} as a fallback."
+                    )
+                    warnings.warn(warning_msg, UserWarning, stacklevel=3)
+                    crs = "epsg:{}".format(srid)
+                else:
+                    if not spatial_ref_sys_df.empty:
+                        auth_name = spatial_ref_sys_df["auth_name"].item()
+                        crs = f"{auth_name}:{srid}"
+                    else:
+                        warning_msg = (
+                            f"Could not find srid {srid} in the "
+                            f"spatial_ref_sys table. "
+                            f"Trying epsg:{srid} as a fallback."
+                        )
+                        warnings.warn(warning_msg, UserWarning, stacklevel=3)
+                        crs = "epsg:{}".format(srid)
 
     return GeoDataFrame(df, crs=crs, geometry=geom_col)
 
@@ -167,7 +191,7 @@ def _read_postgis(
             params=params,
             chunksize=chunksize,
         )
-        return _df_to_geodf(df, geom_col=geom_col, crs=crs)
+        return _df_to_geodf(df, geom_col=geom_col, crs=crs, con=con)
 
     else:
         # read data in chunks and return a generator
@@ -180,7 +204,9 @@ def _read_postgis(
             params=params,
             chunksize=chunksize,
         )
-        return (_df_to_geodf(df, geom_col=geom_col, crs=crs) for df in df_generator)
+        return (
+            _df_to_geodf(df, geom_col=geom_col, crs=crs, con=con) for df in df_generator
+        )
 
 
 def _get_geometry_type(gdf):
@@ -231,7 +257,7 @@ def _get_geometry_type(gdf):
 
 def _get_srid_from_crs(gdf):
     """
-    Get EPSG code from CRS if available. If not, return -1.
+    Get EPSG code from CRS if available. If not, return 0.
     """
 
     # Use geoalchemy2 default for srid
@@ -257,7 +283,7 @@ def _get_srid_from_crs(gdf):
             warnings.warn(warning_msg, UserWarning, stacklevel=2)
 
     if srid is None:
-        srid = -1
+        srid = 0
         warnings.warn(warning_msg, UserWarning, stacklevel=2)
 
     return srid
@@ -437,3 +463,11 @@ def _write_postgis(
             dtype=dtype,
             method=_psql_insert_copy,
         )
+
+
+@lru_cache
+def _get_spatial_ref_sys_df(con, srid):
+    spatial_ref_sys_sql = (
+        f"SELECT srid, auth_name FROM spatial_ref_sys WHERE srid = {srid}"
+    )
+    return pd.read_sql(spatial_ref_sys_sql, con)
diff --git a/geopandas/plotting.py b/geopandas/plotting.py
index 90495833..5c2f416e 100644
--- a/geopandas/plotting.py
+++ b/geopandas/plotting.py
@@ -112,7 +112,15 @@ def _PolygonPatch(polygon, **kwargs):
 
 
 def _plot_polygon_collection(
-    ax, geoms, values=None, color=None, cmap=None, vmin=None, vmax=None, **kwargs
+    ax,
+    geoms,
+    values=None,
+    color=None,
+    cmap=None,
+    vmin=None,
+    vmax=None,
+    autolim=True,
+    **kwargs,
 ):
     """
     Plots a collection of Polygon and MultiPolygon geometries to `ax`
@@ -133,6 +141,8 @@ def _plot_polygon_collection(
         Color to fill the polygons. Cannot be used together with `values`.
     color : single color or sequence of `N` colors
         Sets both `edgecolor` and `facecolor`
+    autolim : bool (default True)
+        Update axes data limits to contain the new geometries.
     **kwargs
         Additional keyword arguments passed to the collection
 
@@ -167,13 +177,21 @@ def _plot_polygon_collection(
         if "norm" not in kwargs:
             collection.set_clim(vmin, vmax)
 
-    ax.add_collection(collection, autolim=True)
+    ax.add_collection(collection, autolim=autolim)
     ax.autoscale_view()
     return collection
 
 
 def _plot_linestring_collection(
-    ax, geoms, values=None, color=None, cmap=None, vmin=None, vmax=None, **kwargs
+    ax,
+    geoms,
+    values=None,
+    color=None,
+    cmap=None,
+    vmin=None,
+    vmax=None,
+    autolim=True,
+    **kwargs,
 ):
     """
     Plots a collection of LineString and MultiLineString geometries to `ax`
@@ -189,6 +207,8 @@ def _plot_linestring_collection(
         have 1:1 correspondence with the geometries (not their components).
     color : single color or sequence of `N` colors
         Cannot be used together with `values`.
+    autolim : bool (default True)
+        Update axes data limits to contain the new geometries.
 
     Returns
     -------
@@ -222,7 +242,7 @@ def _plot_linestring_collection(
         if "norm" not in kwargs:
             collection.set_clim(vmin, vmax)
 
-    ax.add_collection(collection, autolim=True)
+    ax.add_collection(collection, autolim=autolim)
     ax.autoscale_view()
     return collection
 
@@ -291,7 +311,14 @@ def _plot_point_collection(
 
 
 def plot_series(
-    s, cmap=None, color=None, ax=None, figsize=None, aspect="auto", **style_kwds
+    s,
+    cmap=None,
+    color=None,
+    ax=None,
+    figsize=None,
+    aspect="auto",
+    autolim=True,
+    **style_kwds,
 ):
     """
     Plot a GeoSeries.
@@ -327,6 +354,8 @@ def plot_series(
         square appears square in the middle of the plot. This implies an
         Equirectangular projection. If None, the aspect of `ax` won't be changed. It can
         also be set manually (float) as the ratio of y-unit to x-unit.
+    autolim : bool (default True)
+        Update axes data limits to contain the new geometries.
     **style_kwds : dict
         Color options to be passed on to the actual plot function, such
         as ``edgecolor``, ``facecolor``, ``linewidth``, ``markersize``,
@@ -421,7 +450,13 @@ def plot_series(
 
         values_ = values[poly_idx] if cmap else None
         _plot_polygon_collection(
-            ax, polys, values_, facecolor=facecolor, cmap=cmap, **style_kwds
+            ax,
+            polys,
+            values_,
+            facecolor=facecolor,
+            cmap=cmap,
+            autolim=autolim,
+            **style_kwds,
         )
 
     # plot all LineStrings and MultiLineString components in same collection
@@ -431,7 +466,7 @@ def plot_series(
         color_ = expl_color[line_idx] if color_given else color
 
         _plot_linestring_collection(
-            ax, lines, values_, color=color_, cmap=cmap, **style_kwds
+            ax, lines, values_, color=color_, cmap=cmap, autolim=autolim, **style_kwds
         )
 
     # plot all Points in the same collection
@@ -468,6 +503,7 @@ def plot_dataframe(
     classification_kwds=None,
     missing_kwds=None,
     aspect="auto",
+    autolim=True,
     **style_kwds,
 ):
     """
@@ -571,7 +607,8 @@ def plot_dataframe(
         square appears square in the middle of the plot. This implies an
         Equirectangular projection. If None, the aspect of `ax` won't be changed. It can
         also be set manually (float) as the ratio of y-unit to x-unit.
-
+    autolim : bool (default True)
+        Update axes data limits to contain the new geometries.
     **style_kwds : dict
         Style options to be passed on to the actual plot function, such
         as ``edgecolor``, ``facecolor``, ``linewidth``, ``markersize``,
@@ -658,6 +695,7 @@ def plot_dataframe(
             figsize=figsize,
             markersize=markersize,
             aspect=aspect,
+            autolim=autolim,
             **style_kwds,
         )
 
@@ -797,7 +835,14 @@ def plot_dataframe(
     subset = values[poly_idx & np.invert(nan_idx)]
     if not polys.empty:
         _plot_polygon_collection(
-            ax, polys, subset, vmin=mn, vmax=mx, cmap=cmap, **style_kwds
+            ax,
+            polys,
+            subset,
+            vmin=mn,
+            vmax=mx,
+            cmap=cmap,
+            autolim=autolim,
+            **style_kwds,
         )
 
     # plot all LineStrings and MultiLineString components in same collection
@@ -805,7 +850,14 @@ def plot_dataframe(
     subset = values[line_idx & np.invert(nan_idx)]
     if not lines.empty:
         _plot_linestring_collection(
-            ax, lines, subset, vmin=mn, vmax=mx, cmap=cmap, **style_kwds
+            ax,
+            lines,
+            subset,
+            vmin=mn,
+            vmax=mx,
+            cmap=cmap,
+            autolim=autolim,
+            **style_kwds,
         )
 
     # plot all Points in the same collection
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 7a3c70c3..06d7ef74 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -18,6 +18,7 @@ def sjoin(
     lsuffix="left",
     rsuffix="right",
     distance=None,
+    on_attribute=None,
     **kwargs,
 ):
     """Spatial join of two GeoDataFrames.
@@ -50,6 +51,11 @@ def sjoin(
         for the 'dwithin' predicate. If array_like, must be
         one-dimesional with length equal to length of left GeoDataFrame.
         Required if ``predicate='dwithin'``.
+    on_attribute : string, list or tuple
+        Column name(s) to join on as an additional join restriction on top
+        of the spatial predicate. These must be found in both DataFrames.
+        If set, observations are joined only if the predicate applies
+        and values in specified columns match.
 
     Examples
     --------
@@ -103,18 +109,38 @@ def sjoin(
         first = next(iter(kwargs.keys()))
         raise TypeError(f"sjoin() got an unexpected keyword argument '{first}'")
 
-    _basic_checks(left_df, right_df, how, lsuffix, rsuffix)
+    on_attribute = _maybe_make_list(on_attribute)
+
+    _basic_checks(left_df, right_df, how, lsuffix, rsuffix, on_attribute=on_attribute),
 
-    indices = _geom_predicate_query(left_df, right_df, predicate, distance)
+    indices = _geom_predicate_query(
+        left_df, right_df, predicate, distance, on_attribute=on_attribute
+    )
 
     joined, _ = _frame_join(
-        left_df, right_df, indices, None, how, lsuffix, rsuffix, predicate
+        left_df,
+        right_df,
+        indices,
+        None,
+        how,
+        lsuffix,
+        rsuffix,
+        predicate,
+        on_attribute=on_attribute,
     )
 
     return joined
 
 
-def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
+def _maybe_make_list(obj):
+    if isinstance(obj, tuple):
+        return list(obj)
+    if obj is not None and not isinstance(obj, list):
+        return [obj]
+    return obj
+
+
+def _basic_checks(left_df, right_df, how, lsuffix, rsuffix, on_attribute=None):
     """Checks the validity of join input parameters.
 
     `how` must be one of the valid options.
@@ -131,6 +157,8 @@ def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
         left index suffix
     rsuffix : str
         right index suffix
+    on_attribute : list, default None
+        list of column names to merge on along with geometry
     """
     if not isinstance(left_df, GeoDataFrame):
         raise ValueError(
@@ -151,8 +179,28 @@ def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
     if not _check_crs(left_df, right_df):
         _crs_mismatch_warn(left_df, right_df, stacklevel=4)
 
+    if on_attribute:
+        for attr in on_attribute:
+            if (attr not in left_df) and (attr not in right_df):
+                raise ValueError(
+                    f"Expected column {attr} is missing from both of the dataframes."
+                )
+            if attr not in left_df:
+                raise ValueError(
+                    f"Expected column {attr} is missing from the left dataframe."
+                )
+            if attr not in right_df:
+                raise ValueError(
+                    f"Expected column {attr} is missing from the right dataframe."
+                )
+            if attr in (left_df.geometry.name, right_df.geometry.name):
+                raise ValueError(
+                    "Active geometry column cannot be used as an input "
+                    "for on_attribute parameter."
+                )
+
 
-def _geom_predicate_query(left_df, right_df, predicate, distance):
+def _geom_predicate_query(left_df, right_df, predicate, distance, on_attribute=None):
     """Compute geometric comparisons and get matching indices.
 
     Parameters
@@ -161,6 +209,9 @@ def _geom_predicate_query(left_df, right_df, predicate, distance):
     right_df : GeoDataFrame
     predicate : string
         Binary predicate to query.
+    on_attribute: list, default None
+        list of column names to merge on along with geometry
+
 
     Returns
     -------
@@ -200,6 +251,12 @@ def _geom_predicate_query(left_df, right_df, predicate, distance):
         l_idx = l_idx[indexer]
         r_idx = r_idx[indexer]
 
+    if on_attribute:
+        for attr in on_attribute:
+            (l_idx, r_idx), _ = _filter_shared_attribute(
+                left_df, right_df, l_idx, r_idx, attr
+            )
+
     return l_idx, r_idx
 
 
@@ -356,7 +413,15 @@ def _adjust_indexers(indices, distances, original_length, how, predicate):
 
 
 def _frame_join(
-    left_df, right_df, indices, distances, how, lsuffix, rsuffix, predicate
+    left_df,
+    right_df,
+    indices,
+    distances,
+    how,
+    lsuffix,
+    rsuffix,
+    predicate,
+    on_attribute=None,
 ):
     """Join the GeoDataFrames at the DataFrame level.
 
@@ -376,12 +441,18 @@ def _frame_join(
         Suffix to apply to overlapping column names (left GeoDataFrame).
     rsuffix : string
         Suffix to apply to overlapping column names (right GeoDataFrame).
+    on_attribute: list, default None
+        list of column names to merge on along with geometry
+
 
     Returns
     -------
     GeoDataFrame
         Joined GeoDataFrame.
     """
+    if on_attribute:  # avoid renaming or duplicating shared column
+        right_df = right_df.drop(on_attribute, axis=1)
+
     if how in ("inner", "left"):
         right_df = right_df.drop(right_df.geometry.name, axis=1)
     else:  # how == 'right':
@@ -407,7 +478,6 @@ def _frame_join(
     )
     left_df.columns = left_column_names
     right_df.columns = right_column_names
-
     left_index = left_df.columns[:left_nlevels]
     right_index = right_df.columns[:right_nlevels]
 
@@ -444,6 +514,7 @@ def _nearest_query(
     how: str,
     return_distance: bool,
     exclusive: bool,
+    on_attribute: Optional[list] = None,
 ):
     # use the opposite of the join direction for the index
     use_left_as_sindex = how == "right"
@@ -481,9 +552,32 @@ def _nearest_query(
             distances = np.array([], dtype=np.float64)
         else:
             distances = None
+
+    if on_attribute:
+        for attr in on_attribute:
+            (l_idx, r_idx), shared_attribute_rows = _filter_shared_attribute(
+                left_df, right_df, l_idx, r_idx, attr
+            )
+            distances = distances[shared_attribute_rows]
+
     return (l_idx, r_idx), distances
 
 
+def _filter_shared_attribute(left_df, right_df, l_idx, r_idx, attribute):
+    """
+    Returns the indices for the left and right dataframe that share the same entry
+    in the attribute column. Also returns a Boolean `shared_attribute_rows` for rows
+    with the same entry.
+    """
+    shared_attribute_rows = (
+        left_df[attribute].iloc[l_idx].values == right_df[attribute].iloc[r_idx].values
+    )
+
+    l_idx = l_idx[shared_attribute_rows]
+    r_idx = r_idx[shared_attribute_rows]
+    return (l_idx, r_idx), shared_attribute_rows
+
+
 def sjoin_nearest(
     left_df: GeoDataFrame,
     right_df: GeoDataFrame,
@@ -607,6 +701,7 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
     Every operation in GeoPandas is planar, i.e. the potential third
     dimension is not taken into account.
     """
+
     _basic_checks(left_df, right_df, how, lsuffix, rsuffix)
 
     left_df.geometry.values.check_geographic_crs(stacklevel=1)
@@ -615,10 +710,22 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
     return_distance = distance_col is not None
 
     indices, distances = _nearest_query(
-        left_df, right_df, max_distance, how, return_distance, exclusive
+        left_df,
+        right_df,
+        max_distance,
+        how,
+        return_distance,
+        exclusive,
     )
     joined, distances = _frame_join(
-        left_df, right_df, indices, distances, how, lsuffix, rsuffix, None
+        left_df,
+        right_df,
+        indices,
+        distances,
+        how,
+        lsuffix,
+        rsuffix,
+        None,
     )
 
     if return_distance:


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 83303400..4e419a47 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -101,6 +101,7 @@ jobs:
         run: |
           conda install postgis -c conda-forge
           sh ci/scripts/setup_postgres.sh
+          set -o pipefail
           pytest -v -r fEs --color=yes --cov=geopandas --cov-append --cov-report term-missing --cov-report xml geopandas/io/tests/test_sql.py | tee /dev/stderr | if grep SKIPPED >/dev/null;then echo "TESTS SKIPPED, FAILING" && exit 1;fi
 
       - name: Test docstrings
diff --git a/geopandas/io/tests/test_arrow.py b/geopandas/io/tests/test_arrow.py
index 836c7e23..9078191d 100644
--- a/geopandas/io/tests/test_arrow.py
+++ b/geopandas/io/tests/test_arrow.py
@@ -43,6 +43,8 @@ DATA_PATH = pathlib.Path(os.path.dirname(__file__)) / "data"
 pyarrow = pytest.importorskip("pyarrow")
 
 import pyarrow.compute as pc
+import pyarrow.parquet as pq
+from pyarrow import feather
 
 
 @pytest.fixture(
@@ -433,6 +435,37 @@ def test_index(tmpdir, file_format, naturalearth_lowres):
     assert_geodataframe_equal(df.reset_index(drop=True), pq_df)
 
 
+def test_column_order(tmpdir, file_format, naturalearth_lowres):
+    """The order of columns should be preserved in the output."""
+    reader, writer = file_format
+
+    df = read_file(naturalearth_lowres)
+    df = df.set_index("iso_a3")
+    df["geom2"] = df.geometry.representative_point()
+    table = _geopandas_to_arrow(df)
+    custom_column_order = [
+        "iso_a3",
+        "geom2",
+        "pop_est",
+        "continent",
+        "name",
+        "geometry",
+        "gdp_md_est",
+    ]
+    table = table.select(custom_column_order)
+
+    if reader is read_parquet:
+        filename = os.path.join(str(tmpdir), "test_column_order.pq")
+        pq.write_table(table, filename)
+    else:
+        filename = os.path.join(str(tmpdir), "test_column_order.feather")
+        feather.write_feather(table, filename)
+
+    result = reader(filename)
+    assert list(result.columns) == custom_column_order[1:]
+    assert_geodataframe_equal(result, df[custom_column_order[1:]])
+
+
 @pytest.mark.parametrize("compression", ["snappy", "gzip", "brotli", None])
 def test_parquet_compression(compression, tmpdir, naturalearth_lowres):
     """Using compression options should not raise errors, and should
@@ -854,46 +887,6 @@ def test_write_spec_version(tmpdir, format, schema_version):
         assert metadata["columns"]["geometry"]["geometry_types"] == ["Polygon"]
 
 
-@pytest.mark.parametrize(
-    "format,version", product(["feather", "parquet"], [None] + SUPPORTED_VERSIONS)
-)
-def test_write_deprecated_version_parameter(tmpdir, format, version):
-    if format == "feather":
-        from pyarrow.feather import read_table
-
-        version = version or 2
-
-    else:
-        from pyarrow.parquet import read_table
-
-        version = version or "2.6"
-
-    filename = os.path.join(str(tmpdir), f"test.{format}")
-    gdf = geopandas.GeoDataFrame(geometry=[box(0, 0, 10, 10)], crs="EPSG:4326")
-    write = getattr(gdf, f"to_{format}")
-
-    if version in SUPPORTED_VERSIONS:
-        with pytest.warns(
-            FutureWarning,
-            match="the `version` parameter has been replaced with `schema_version`",
-        ):
-            write(filename, version=version)
-
-    else:
-        # no warning raised if not one of the captured versions
-        write(filename, version=version)
-
-    table = read_table(filename)
-    metadata = json.loads(table.schema.metadata[b"geo"])
-
-    if version in SUPPORTED_VERSIONS:
-        # version is captured as a parameter
-        assert metadata["version"] == version
-    else:
-        # version is passed to underlying writer
-        assert metadata["version"] == METADATA_VERSION
-
-
 @pytest.mark.parametrize("version", ["0.1.0", "0.4.0", "1.0.0-beta.1"])
 def test_read_versioned_file(version):
     """
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index 71036746..8072cf6c 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -26,7 +26,6 @@ from pandas.testing import assert_frame_equal, assert_series_equal
 
 try:
     import pyogrio
-
 except ImportError:
     pyogrio = False
 
@@ -1298,6 +1297,54 @@ def test_write_read_file(test_file, engine):
     os.remove(os.path.expanduser(test_file))
 
 
+@pytest.mark.skipif(fiona is False, reason="Fiona not available")
+@pytest.mark.skipif(FIONA_GE_19, reason="Fiona >= 1.9 supports metadata")
+def test_to_file_metadata_unsupported_fiona_version(tmp_path, df_points):
+    metadata = {"title": "test"}
+    tmp_file = tmp_path / "test.gpkg"
+    match = "'metadata' keyword is only supported for Fiona >= 1.9"
+    with pytest.raises(NotImplementedError, match=match):
+        df_points.to_file(tmp_file, driver="GPKG", engine="fiona", metadata=metadata)
+
+
+@pytest.mark.skipif(not FIONA_GE_19, reason="only Fiona >= 1.9 supports metadata")
+def test_to_file_metadata_supported_fiona_version(tmp_path, df_points):
+    metadata = {"title": "test"}
+    tmp_file = tmp_path / "test.gpkg"
+
+    df_points.to_file(tmp_file, driver="GPKG", engine="fiona", metadata=metadata)
+
+    # Check that metadata is written to the file
+    with fiona.open(tmp_file) as src:
+        tags = src.tags()
+        assert tags == metadata
+
+
+@pytest.mark.skipif(pyogrio is False, reason="Pyogrio not available")
+def test_to_file_metadata_pyogrio(tmp_path, df_points):
+    metadata = {"title": "test"}
+    tmp_file = tmp_path / "test.gpkg"
+
+    df_points.to_file(tmp_file, driver="GPKG", engine="pyogrio", metadata=metadata)
+
+    # Check that metadata is written to the file
+    info = pyogrio.read_info(tmp_file)
+    layer_metadata = info["layer_metadata"]
+    assert layer_metadata == metadata
+
+
+@pytest.mark.parametrize(
+    "driver, ext", [("ESRI Shapefile", ".shp"), ("GeoJSON", ".geojson")]
+)
+def test_to_file_metadata_unsupported_driver(driver, ext, tmpdir, df_points, engine):
+    metadata = {"title": "Test"}
+    tempfilename = os.path.join(str(tmpdir), "test" + ext)
+    with pytest.raises(
+        NotImplementedError, match="'metadata' keyword is only supported for"
+    ):
+        df_points.to_file(tempfilename, driver=driver, metadata=metadata)
+
+
 def test_multiple_geom_cols_error(tmpdir, df_nybb):
     df_nybb["geom2"] = df_nybb.geometry
     with pytest.raises(ValueError, match="GeoDataFrame contains multiple geometry"):
diff --git a/geopandas/io/tests/test_geoarrow.py b/geopandas/io/tests/test_geoarrow.py
index c9502264..0087e6b0 100644
--- a/geopandas/io/tests/test_geoarrow.py
+++ b/geopandas/io/tests/test_geoarrow.py
@@ -12,6 +12,7 @@ from shapely import MultiPoint, Point, box
 from geopandas import GeoDataFrame, GeoSeries
 
 import pytest
+from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
 
 pytest.importorskip("pyarrow")
 import pyarrow as pa
@@ -199,6 +200,10 @@ def test_geoarrow_export(geometry_type, dim, geometry_encoding, interleaved):
     assert result_arr.equals(expected["geometry"].chunk(0))
 
 
+@pytest.mark.skipif(
+    Version(shapely.__version__) < Version("2.0.2"),
+    reason="from_ragged_array failing with read-only array input",
+)
 @pytest.mark.parametrize("encoding", ["WKB", "geoarrow"])
 def test_geoarrow_multiple_geometry_crs(encoding):
     pytest.importorskip("pyproj")
@@ -216,6 +221,11 @@ def test_geoarrow_multiple_geometry_crs(encoding):
     )
     assert json.loads(meta2["crs"])["id"]["code"] == 3857
 
+    roundtripped = GeoDataFrame.from_arrow(result)
+    assert_geodataframe_equal(gdf, roundtripped)
+    assert gdf.geometry.crs == "epsg:4326"
+    assert gdf.geom2.crs == "epsg:3857"
+
 
 @pytest.mark.parametrize("encoding", ["WKB", "geoarrow"])
 def test_geoarrow_series_name_crs(encoding):
@@ -361,3 +371,167 @@ def test_geoarrow_export_with_extension_types(geometry_type, dim):
 
         result3 = pa_table(df.to_arrow(geometry_encoding="geoarrow", interleaved=False))
         assert isinstance(result3["geometry"].type, pa.ExtensionType)
+
+
+@pytest.mark.skipif(
+    Version(shapely.__version__) < Version("2.0.2"),
+    reason="from_ragged_array failing with read-only array input",
+)
+@pytest.mark.parametrize("dim", ["xy", "xyz"])
+@pytest.mark.parametrize(
+    "geometry_type",
+    [
+        "point",
+        "linestring",
+        "polygon",
+        "multipoint",
+        "multilinestring",
+        "multipolygon",
+    ],
+)
+def test_geoarrow_import(geometry_type, dim):
+    base_path = DATA_PATH / "geoarrow"
+    suffix = geometry_type + ("_z" if dim == "xyz" else "")
+
+    # Read the example data
+    df = feather.read_feather(base_path / f"example-{suffix}-wkb.arrow")
+    df["geometry"] = GeoSeries.from_wkb(df["geometry"])
+    df = GeoDataFrame(df)
+    df.geometry.crs = None
+
+    table1 = feather.read_table(base_path / f"example-{suffix}-wkb.arrow")
+    result1 = GeoDataFrame.from_arrow(table1)
+    assert_geodataframe_equal(result1, df)
+
+    table2 = feather.read_table(base_path / f"example-{suffix}-interleaved.arrow")
+    result2 = GeoDataFrame.from_arrow(table2)
+    assert_geodataframe_equal(result2, df)
+
+    table3 = feather.read_table(base_path / f"example-{suffix}.arrow")
+    result3 = GeoDataFrame.from_arrow(table3)
+    assert_geodataframe_equal(result3, df)
+
+
+@pytest.mark.skipif(
+    Version(shapely.__version__) < Version("2.0.2"),
+    reason="from_ragged_array failing with read-only array input",
+)
+@pytest.mark.parametrize("encoding", ["WKB", "geoarrow"])
+def test_geoarrow_import_geometry_column(encoding):
+    pytest.importorskip("pyproj")
+    # ensure each geometry column has its own crs
+    gdf = GeoDataFrame(geometry=[box(0, 0, 10, 10)])
+    gdf["centroid"] = gdf.geometry.centroid
+
+    result = GeoDataFrame.from_arrow(pa_table(gdf.to_arrow(geometry_encoding=encoding)))
+    assert_geodataframe_equal(result, gdf)
+    assert result.active_geometry_name == "geometry"
+
+    result = GeoDataFrame.from_arrow(
+        pa_table(gdf[["centroid"]].to_arrow(geometry_encoding=encoding))
+    )
+    assert result.active_geometry_name == "centroid"
+
+    result = GeoDataFrame.from_arrow(
+        pa_table(gdf.to_arrow(geometry_encoding=encoding)), geometry="centroid"
+    )
+    assert result.active_geometry_name == "centroid"
+    assert_geodataframe_equal(result, gdf.set_geometry("centroid"))
+
+
+def test_geoarrow_import_missing_geometry():
+    pytest.importorskip("pyarrow", minversion="14.0.0")
+
+    table = pa.table({"a": [0, 1, 2], "b": [0.1, 0.2, 0.3]})
+    with pytest.raises(ValueError, match="No geometry column found"):
+        GeoDataFrame.from_arrow(table)
+
+    with pytest.raises(ValueError, match="No GeoArrow geometry field found"):
+        GeoSeries.from_arrow(table["a"].chunk(0))
+
+
+def test_geoarrow_import_capsule_interface():
+    # ensure we can import non-pyarrow object
+    pytest.importorskip("pyarrow", minversion="14.0.0")
+    gdf = GeoDataFrame({"col": [1]}, geometry=[box(0, 0, 10, 10)])
+
+    result = GeoDataFrame.from_arrow(gdf.to_arrow())
+    assert_geodataframe_equal(result, gdf)
+
+
+@pytest.mark.parametrize("dim", ["xy", "xyz"])
+@pytest.mark.parametrize(
+    "geometry_type",
+    ["point", "linestring", "polygon", "multipoint", "multilinestring", "multipolygon"],
+)
+def test_geoarrow_import_from_extension_types(geometry_type, dim):
+    # ensure the exported data can be imported by geoarrow-pyarrow and are
+    # recognized as extension types
+    pytest.importorskip("pyproj")
+    base_path = DATA_PATH / "geoarrow"
+    suffix = geometry_type + ("_z" if dim == "xyz" else "")
+
+    # Read the example data
+    df = feather.read_feather(base_path / f"example-{suffix}-wkb.arrow")
+    df["geometry"] = GeoSeries.from_wkb(df["geometry"])
+    df = GeoDataFrame(df, crs="EPSG:3857")
+
+    pytest.importorskip("geoarrow.pyarrow")
+
+    with with_geoarrow_extension_types():
+        result1 = GeoDataFrame.from_arrow(
+            pa_table(df.to_arrow(geometry_encoding="WKB"))
+        )
+        assert_geodataframe_equal(result1, df)
+
+        result2 = GeoDataFrame.from_arrow(
+            pa_table(df.to_arrow(geometry_encoding="geoarrow"))
+        )
+        assert_geodataframe_equal(result2, df)
+
+        result3 = GeoDataFrame.from_arrow(
+            pa_table(df.to_arrow(geometry_encoding="geoarrow", interleaved=False))
+        )
+        assert_geodataframe_equal(result3, df)
+
+
+def test_geoarrow_import_geoseries():
+    pytest.importorskip("pyproj")
+    gp = pytest.importorskip("geoarrow.pyarrow")
+    ser = GeoSeries.from_wkt(["POINT (1 1)", "POINT (2 2)"], crs="EPSG:3857")
+
+    with with_geoarrow_extension_types():
+        arr = gp.array(ser.to_arrow(geometry_encoding="WKB"))
+        result = GeoSeries.from_arrow(arr)
+        assert_geoseries_equal(result, ser)
+
+        arr = gp.array(ser.to_arrow(geometry_encoding="geoarrow"))
+        result = GeoSeries.from_arrow(arr)
+        assert_geoseries_equal(result, ser)
+
+        # the name is lost when going through a pyarrow.Array
+        ser.name = "name"
+        arr = gp.array(ser.to_arrow())
+        result = GeoSeries.from_arrow(arr)
+        assert result.name is None
+        # we can specify the name as one of the kwargs
+        result = GeoSeries.from_arrow(arr, name="test")
+        assert_geoseries_equal(result, ser)
+
+
+def test_geoarrow_import_unknown_geoarrow_type():
+    gdf = GeoDataFrame({"col": [1]}, geometry=[box(0, 0, 10, 10)])
+    table = pa_table(gdf.to_arrow())
+    schema = table.schema
+    new_field = schema.field("geometry").with_metadata(
+        {
+            b"ARROW:extension:name": b"geoarrow.unknown",
+            b"ARROW:extension:metadata": b"{}",
+        }
+    )
+
+    new_schema = pa.schema([schema.field(0), new_field])
+    new_table = table.cast(new_schema)
+
+    with pytest.raises(TypeError, match="Unknown GeoArrow extension type"):
+        GeoDataFrame.from_arrow(new_table)
diff --git a/geopandas/io/tests/test_sql.py b/geopandas/io/tests/test_sql.py
index 524e7e49..d394098f 100644
--- a/geopandas/io/tests/test_sql.py
+++ b/geopandas/io/tests/test_sql.py
@@ -14,11 +14,17 @@ import pandas as pd
 import geopandas
 import geopandas._compat as compat
 from geopandas import GeoDataFrame, read_file, read_postgis
+from geopandas._compat import HAS_PYPROJ
 from geopandas.io.sql import _get_conn as get_conn
 from geopandas.io.sql import _write_postgis as write_postgis
 
 import pytest
-from geopandas.tests.util import create_postgis, create_spatialite, validate_boro_df
+from geopandas.tests.util import (
+    create_postgis,
+    create_spatialite,
+    mock,
+    validate_boro_df,
+)
 
 try:
     from sqlalchemy import text
@@ -773,6 +779,22 @@ class TestIO:
         with pytest.raises(ValueError, match="CRS of the target table"):
             write_postgis(df_nybb2, con=engine, name=table, if_exists="append")
 
+    @pytest.mark.parametrize("engine_postgis", POSTGIS_DRIVERS, indirect=True)
+    def test_append_without_crs(self, engine_postgis, df_nybb):
+        # This test was included in #3328 when the default value for no
+        # CRS was changed from an SRID of -1 to 0. This resolves issues
+        # of appending dataframes to postgis that have no CRS as postgis
+        # no CRS value is 0.
+        engine = engine_postgis
+        df_nybb = df_nybb.set_crs(None, allow_override=True)
+        table = "nybb"
+
+        write_postgis(df_nybb, con=engine, name=table, if_exists="replace")
+        # append another dataframe with no crs
+
+        df_nybb2 = df_nybb
+        write_postgis(df_nybb2, con=engine, name=table, if_exists="append")
+
     @pytest.mark.parametrize("engine_postgis", POSTGIS_DRIVERS, indirect=True)
     @pytest.mark.xfail(
         compat.PANDAS_GE_20 and not compat.PANDAS_GE_202,
@@ -788,3 +810,69 @@ class TestIO:
 
         with pytest.raises(ValueError):
             read_postgis(sql, engine, geom_col="geom")
+
+    @pytest.mark.parametrize("connection_postgis", POSTGIS_DRIVERS, indirect=True)
+    def test_read_non_epsg_crs(self, connection_postgis, df_nybb):
+        con = connection_postgis
+        df_nybb = df_nybb.to_crs(crs="esri:54052")
+        create_postgis(con, df_nybb, srid=54052)
+
+        sql = "SELECT * FROM nybb;"
+        df = read_postgis(sql, con)
+        validate_boro_df(df)
+        assert df.crs == "ESRI:54052"
+
+    @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not installed")
+    @mock.patch("shapely.get_srid")
+    @pytest.mark.parametrize("connection_postgis", POSTGIS_DRIVERS, indirect=True)
+    def test_read_srid_not_in_table(self, mock_get_srid, connection_postgis, df_nybb):
+        # mock a non-existent srid for edge case if shapely has an srid
+        # not present in postgis table.
+        pyproj = pytest.importorskip("pyproj")
+
+        mock_get_srid.return_value = 99999
+
+        con = connection_postgis
+        df_nybb = df_nybb.to_crs(crs="epsg:4326")
+        create_postgis(con, df_nybb)
+
+        sql = "SELECT * FROM nybb;"
+        with pytest.raises(pyproj.exceptions.CRSError, match="crs not found"):
+            with pytest.warns(UserWarning, match="Could not find srid 99999"):
+                read_postgis(sql, con)
+
+    @mock.patch("geopandas.io.sql._get_spatial_ref_sys_df")
+    @pytest.mark.parametrize("connection_postgis", POSTGIS_DRIVERS, indirect=True)
+    def test_read_no_spatial_ref_sys_table_in_postgis(
+        self, mock_get_spatial_ref_sys_df, connection_postgis, df_nybb
+    ):
+        # mock for a non-existent spatial_ref_sys database
+
+        mock_get_spatial_ref_sys_df.side_effect = pd.errors.DatabaseError
+
+        con = connection_postgis
+        df_nybb = df_nybb.to_crs(crs="epsg:4326")
+        create_postgis(con, df_nybb, srid=4326)
+
+        sql = "SELECT * FROM nybb;"
+        with pytest.warns(
+            UserWarning, match="Could not find the spatial reference system table"
+        ):
+            df = read_postgis(sql, con)
+
+        assert df.crs == "EPSG:4326"
+
+    @pytest.mark.parametrize("connection_postgis", POSTGIS_DRIVERS, indirect=True)
+    def test_read_non_epsg_crs_chunksize(self, connection_postgis, df_nybb):
+        """Test chunksize argument with non epsg crs"""
+        chunksize = 2
+        con = connection_postgis
+        df_nybb = df_nybb.to_crs(crs="esri:54052")
+
+        create_postgis(con, df_nybb, srid=54052)
+
+        sql = "SELECT * FROM nybb;"
+        df = pd.concat(read_postgis(sql, con, chunksize=chunksize))
+
+        validate_boro_df(df)
+        assert df.crs == "ESRI:54052"
diff --git a/geopandas/tests/test_explore.py b/geopandas/tests/test_explore.py
index f1b15d65..d2f161cc 100644
--- a/geopandas/tests/test_explore.py
+++ b/geopandas/tests/test_explore.py
@@ -1004,3 +1004,37 @@ class TestExplore:
         df.loc[0, df.geometry.name] = shapely.Point()
         m = df.explore()
         self._fetch_map_string(m)
+
+    def test_all_empty(self):
+        with_crs = gpd.GeoDataFrame(
+            geometry=[shapely.Point(), shapely.Point()], crs=4326
+        )
+        with pytest.warns(
+            UserWarning,
+            match="The GeoSeries you are attempting to plot is composed of empty",
+        ):
+            m = with_crs.explore()
+        out_str = self._fetch_map_string(m)
+        if HAS_PYPROJ:
+            assert "center:[0.0,0.0],crs:L.CRS.EPSG3857" in out_str
+
+        no_crs = gpd.GeoDataFrame(geometry=[shapely.Point(), shapely.Point()])
+        with pytest.warns(
+            UserWarning,
+            match="The GeoSeries you are attempting to plot is composed of empty",
+        ):
+            m = no_crs.explore()
+        out_str = self._fetch_map_string(m)
+        assert "center:[0.0,0.0],crs:L.CRS.Simple" in out_str
+
+    def test_add_all_empty_named_index(self):
+        gdf1 = gpd.GeoDataFrame(geometry=[shapely.Point(0, 0), shapely.Point(1, 1)])
+        gdf2 = gpd.GeoDataFrame(geometry=[shapely.Point(), shapely.Point()])
+        m = gdf1.rename_axis(index="index_name").explore()
+        with pytest.warns(
+            UserWarning,
+            match="The GeoSeries you are attempting to plot is composed of empty",
+        ):
+            m = gdf2.rename_axis(index="index_name").explore(m=m, color="red")
+        out_str = self._fetch_map_string(m)
+        assert "center:[0.5,0.5],crs:L.CRS.Simple" in out_str
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 8efb81ab..c9d74b1f 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -1323,18 +1323,17 @@ class TestConstructor:
         check_geodataframe(gdf)
         assert list(gdf.columns) == ["geometry", "a"]
 
-    def test_preserve_series_name(self):
+    def test_do_not_preserve_series_name_in_constructor(self):
+        # GH3337
+        # GeoDataFrame(... geometry=...) should always create geom col "geometry"
         geoms = [Point(1, 1), Point(2, 2), Point(3, 3)]
         gs = GeoSeries(geoms)
         gdf = GeoDataFrame({"a": [1, 2, 3]}, geometry=gs)
-
         check_geodataframe(gdf, geometry_column="geometry")
-
-        geoms = [Point(1, 1), Point(2, 2), Point(3, 3)]
+        # still get "geometry", even with custom geoseries name
         gs = GeoSeries(geoms, name="my_geom")
         gdf = GeoDataFrame({"a": [1, 2, 3]}, geometry=gs)
-
-        check_geodataframe(gdf, geometry_column="my_geom")
+        check_geodataframe(gdf, geometry_column="geometry")
 
     def test_overwrite_geometry(self):
         # GH602
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 4a9601ad..7230c3ea 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -689,14 +689,16 @@ def test_groupby_metadata(crs, geometry_name):
     df.groupby("value2").apply(func, **kwargs)
     # selecting the non-group columns -> no need to pass the keyword
     if (
-        compat.PANDAS_GE_21 if geometry_name == "geometry" else compat.PANDAS_GE_20
-    ) and not compat.PANDAS_GE_22:
+        compat.PANDAS_GE_22
+        or (compat.PANDAS_GE_20 and geometry_name == "geometry")
+        or not compat.PANDAS_GE_20
+    ):
+        df.groupby("value2")[[geometry_name, "value1"]].apply(func)
+    else:
         # https://github.com/geopandas/geopandas/pull/2966#issuecomment-1878816712
-        # with pandas 2.0 and 2.1 this is failing
+        # with pandas 2.0 and 2.1 with geom col != geometry this is failing
         with pytest.raises(AttributeError):
             df.groupby("value2")[[geometry_name, "value1"]].apply(func)
-    else:
-        df.groupby("value2")[[geometry_name, "value1"]].apply(func)
 
     # actual test with functionality
     res = df.groupby("value2").apply(
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index d47b8987..9df655eb 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -494,6 +494,28 @@ class TestLineStringPlotting:
         )
         self.df3 = GeoDataFrame({"geometry": self.linearrings, "values": values})
 
+    def test_autolim_false(self):
+        """Test linestring plot preserving axes limits."""
+        ax = self.lines[: self.N // 2].plot()
+        ylim = ax.get_ylim()
+        self.lines.plot(ax=ax, autolim=False)
+        assert ax.get_ylim() == ylim
+        ax = self.df[: self.N // 2].plot()
+        ylim = ax.get_ylim()
+        self.df.plot(ax=ax, autolim=False)
+        assert ax.get_ylim() == ylim
+
+    def test_autolim_true(self):
+        """Test linestring plot autoscaling axes limits."""
+        ax = self.lines[: self.N // 2].plot()
+        ylim = ax.get_ylim()
+        self.lines.plot(ax=ax, autolim=True)
+        assert ax.get_ylim() != ylim
+        ax = self.df[: self.N // 2].plot()
+        ylim = ax.get_ylim()
+        self.df.plot(ax=ax, autolim=True)
+        assert ax.get_ylim() != ylim
+
     def test_single_color(self):
         ax = self.lines.plot(color="green")
         _check_colors(self.N, ax.collections[0].get_colors(), ["green"] * self.N)
@@ -637,6 +659,28 @@ class TestPolygonPlotting:
         df_nan = GeoDataFrame({"geometry": t3, "values": [np.nan]})
         self.df3 = pd.concat([self.df, df_nan])
 
+    def test_autolim_false(self):
+        """Test polygon plot preserving axes limits."""
+        ax = self.polys[:1].plot()
+        xlim = ax.get_xlim()
+        self.polys.plot(ax=ax, autolim=False)
+        assert ax.get_xlim() == xlim
+        ax = self.df[:1].plot()
+        xlim = ax.get_xlim()
+        self.df.plot(ax=ax, autolim=False)
+        assert ax.get_xlim() == xlim
+
+    def test_autolim_true(self):
+        """Test polygon plot autoscaling axes limits."""
+        ax = self.polys[:1].plot()
+        xlim = ax.get_xlim()
+        self.polys.plot(ax=ax, autolim=True)
+        assert ax.get_xlim() != xlim
+        ax = self.df[:1].plot()
+        xlim = ax.get_xlim()
+        self.df.plot(ax=ax, autolim=True)
+        assert ax.get_xlim() != xlim
+
     def test_single_color(self):
         ax = self.polys.plot(color="green")
         _check_colors(2, ax.collections[0].get_facecolors(), ["green"] * 2)
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index 0453355b..0e44f87b 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -96,6 +96,52 @@ def dfs(request):
     return [request.param, df1, df2, expected]
 
 
+@pytest.fixture()
+def dfs_shared_attribute():
+    geo_left = [
+        Point(0, 0),
+        Point(1, 1),
+        Point(2, 2),
+        Point(3, 3),
+        Point(4, 4),
+        Point(5, 5),
+        Point(6, 6),
+        Point(7, 7),
+    ]
+    geo_right = [
+        Point(0, 0),
+        Point(1, 1),
+        Point(2, 2),
+        Point(3, 3),
+        Point(4, 4),
+        Point(5, 5),
+        Point(6, 6),
+        Point(7, 7),
+    ]
+    attr_tracker = ["A", "B", "C", "D", "E", "F", "G", "H"]
+
+    left_gdf = geopandas.GeoDataFrame(
+        {
+            "geometry": geo_left,
+            "attr_tracker": attr_tracker,
+            "duplicate_column": [0, 1, 2, 3, 4, 5, 6, 7],
+            "attr1": [True, True, True, True, True, True, True, True],
+            "attr2": [True, True, True, True, True, True, True, True],
+        }
+    )
+
+    right_gdf = geopandas.GeoDataFrame(
+        {
+            "geometry": geo_right,
+            "duplicate_column": [0, 1, 2, 3, 4, 5, 6, 7],
+            "attr1": [True, True, False, False, True, True, False, False],
+            "attr2": [True, True, False, False, False, False, False, False],
+        }
+    )
+
+    return left_gdf, right_gdf
+
+
 class TestSpatialJoin:
     @pytest.mark.parametrize(
         "how, lsuffix, rsuffix, expected_cols",
@@ -147,7 +193,6 @@ class TestSpatialJoin:
         index, df1, df2, expected = dfs
 
         res = sjoin(df1, df2, how="inner", predicate=predicate)
-
         exp = expected[predicate].dropna().copy()
         exp = exp.drop("geometry_y", axis=1).rename(columns={"geometry_x": "geometry"})
         exp[["df1", "df2"]] = exp[["df1", "df2"]].astype("int64")
@@ -168,7 +213,6 @@ class TestSpatialJoin:
             exp.index.names = df1.index.names
         if index == "named-multi-index":
             exp = exp.set_index(["df1_ix1", "df1_ix2"])
-
         assert_frame_equal(res, exp)
 
     @pytest.mark.parametrize(
@@ -424,6 +468,104 @@ class TestSpatialJoin:
         joined = sjoin(pts, polys, predicate=predicate, how="left")
         assert_index_equal(joined.index, pts.index)
 
+    def test_sjoin_shared_attribute(self, naturalearth_lowres, naturalearth_cities):
+        countries = read_file(naturalearth_lowres)
+        cities = read_file(naturalearth_cities)
+        countries = countries[["geometry", "name"]].rename(columns={"name": "country"})
+
+        # Add first letter of country/city as an attribute column to be compared
+        countries["firstLetter"] = countries["country"].astype(str).str[0]
+        cities["firstLetter"] = cities["name"].astype(str).str[0]
+
+        result = sjoin(cities, countries, on_attribute="firstLetter")
+        assert (
+            result["country"].astype(str).str[0] == result["name"].astype(str).str[0]
+        ).all()
+        assert result.shape == (23, 5)
+
+    @pytest.mark.parametrize(
+        "attr1_key_change_dict, attr2_key_change_dict",
+        [
+            pytest.param(
+                {True: "merge", False: "no_merge"},
+                {True: "merge", False: "no_merge"},
+                id="merge on string attributes",
+            ),
+            pytest.param(
+                {True: 2, False: 1},
+                {True: 2, False: 1},
+                id="merge on integer attributes",
+            ),
+            pytest.param(
+                {True: True, False: False},
+                {True: True, False: False},
+                id="merge on boolean attributes",
+            ),
+            pytest.param(
+                {True: True, False: False},
+                {True: "merge", False: "no_merge"},
+                id="merge on mixed attributes",
+            ),
+        ],
+    )
+    def test_sjoin_multiple_attributes_datatypes(
+        self, dfs_shared_attribute, attr1_key_change_dict, attr2_key_change_dict
+    ):
+        left_gdf, right_gdf = dfs_shared_attribute
+        left_gdf["attr1"] = left_gdf["attr1"].map(attr1_key_change_dict)
+        left_gdf["attr2"] = left_gdf["attr2"].map(attr2_key_change_dict)
+        right_gdf["attr1"] = right_gdf["attr1"].map(attr1_key_change_dict)
+        right_gdf["attr2"] = right_gdf["attr2"].map(attr2_key_change_dict)
+
+        joined = sjoin(left_gdf, right_gdf, on_attribute=("attr1", "attr2"))
+        assert (["A", "B"] == joined["attr_tracker"].values).all()
+
+    def test_sjoin_multiple_attributes_check_header(self, dfs_shared_attribute):
+        left_gdf, right_gdf = dfs_shared_attribute
+        joined = sjoin(left_gdf, right_gdf, on_attribute=["attr1"])
+
+        assert (["A", "B", "E", "F"] == joined["attr_tracker"].values).all()
+        assert {"attr2_left", "attr2_right", "attr1"}.issubset(joined.columns)
+        assert "attr1_left" not in joined
+
+    def test_sjoin_error_column_does_not_exist(self, dfs_shared_attribute):
+        left_gdf, right_gdf = dfs_shared_attribute
+        right_gdf_dropped_attr = right_gdf.drop("attr1", axis=1)
+        left_gdf_dropped_attr = left_gdf.drop("attr1", axis=1)
+
+        with pytest.raises(
+            ValueError,
+            match="Expected column attr1 is missing from the right dataframe.",
+        ):
+            sjoin(left_gdf, right_gdf_dropped_attr, on_attribute="attr1")
+
+        with pytest.raises(
+            ValueError,
+            match="Expected column attr1 is missing from the left dataframe.",
+        ):
+            sjoin(left_gdf_dropped_attr, right_gdf, on_attribute="attr1")
+
+        with pytest.raises(
+            ValueError,
+            match="Expected column attr1 is missing from both of the dataframes.",
+        ):
+            sjoin(left_gdf_dropped_attr, right_gdf_dropped_attr, on_attribute="attr1")
+
+    def test_sjoin_error_use_geometry_column(self, dfs_shared_attribute):
+        left_gdf, right_gdf = dfs_shared_attribute
+        with pytest.raises(
+            ValueError,
+            match="Active geometry column cannot be used as an input for "
+            "on_attribute parameter.",
+        ):
+            sjoin(left_gdf, right_gdf, on_attribute="geometry")
+        with pytest.raises(
+            ValueError,
+            match="Active geometry column cannot be used as an input for "
+            "on_attribute parameter.",
+        ):
+            sjoin(left_gdf, right_gdf, on_attribute=["attr1", "geometry"])
+
 
 class TestIndexNames:
     @pytest.mark.parametrize("how", ["inner", "left", "right"])