Repository: geopandas/geopandas

Issue:
ENH: Return make_valid to overlay and use to skip the buffer(0) calls
In https://github.com/geopandas/geopandas/pull/1553 `make_valid` was removed. However, in our code the unexpected calls to `buffer(0)` exposed this bug in GEOS https://github.com/Toblerity/Shapely/issues/1069. 

I propose `make_valid=True` is restored and the `buffer(0)` calls are skipped if `make_valid=False`. In this case documentation could warm the user to ensure that geometries are valid before the call.

Discussion:
bretttully: @martinfleis I'm getting strange behaviour with the `geom_almost_equals` function. 

I am running: `pytest --pdb -vvs geopandas/tests/test_overlay.py::test_overlay_strict` when it hits the assertion error:

```python
(Pdb) left.geom_almost_equals(right)
0    False
1    False
2    False
dtype: bool

(Pdb) left.geom_equals(right)
0    True
1    True
2    True
dtype: bool

(Pdb) left.at[0].wkt
'POLYGON ((3 1, 1 1, 1 3, 3 3, 3 1))'

(Pdb) right.at[0].wkt
'POLYGON ((1 1, 1 3, 3 3, 3 1, 1 1))'

(Pdb) left.at[0].equals(right.at[0])
True

(Pdb) left.at[0].almost_equals(right.at[0])
False
```

It seems wrong that `almost_equals` fails when `equals` succeeds... What am I missing?
martinfleis: This is a bit confusing, but drawing from shapely's documentation:
 
> `eqauls` Returns True if the set-theoretic boundary, interior, and exterior of the object coincide with those of the other.

> `almost_equals` Returns True if the object is approximately equal to the other at all points to specified decimal place precision.

Which means that `equals` does not care about the order of coordinates, while `almost_equals` does (they can be slightly off but needs to be in the same order). This is a behaviour we inherit from shapely.

In pygeos, this is more clear, because `almost_equals` is called `equals_exact`:

> Returns True if A and B are structurally equal.

> This method uses exact coordinate equality, which requires coordinates to be equal (within specified tolerance) and and in the same order for all components of a geometry. This is in contrast with the equals function which uses spatial (topological) equality.
bretttully: @martinfleis that is a bit confusing.

See my proposed solution in the commit above. Not sure if this is acceptable, or if I should be going through all of the shapes in the overlay test and updating them?
martinfleis: I think that it is in principle a correct solution. The point of `check_less_precise` in `assert_geodataframe_equal` is to be a bit more forgiving, which is retained in your solution. Can you just update docstrings of `assert_geoseries_equal` and `assert_geodataframe_equal` to make clear that if `check_less_precise=True` both `geom_equals ` and `geom_almost_equals` are used, not only one?

https://github.com/geopandas/geopandas/blob/79aeefacd71dccbe86fac741d1e4809fba0e3d4e/geopandas/testing.py#L85-L86

https://github.com/geopandas/geopandas/blob/79aeefacd71dccbe86fac741d1e4809fba0e3d4e/geopandas/testing.py#L164-L165
martinfleis: With changes merged in #1826 we can maybe revert changes to `assert_geoseries_equal` and use `normalize` keyword in overlay tests. That may clean-up the PR considerably.

Patch:
diff --git a/.gitignore b/.gitignore
index 8eefebdf..5e9c00e8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,6 +34,7 @@ nosetests.xml
 coverage.xml
 *.cover
 .hypothesis/
+result_images
 
 # Sphinx documentation
 doc/_build/
diff --git a/benchmarks/sindex.py b/benchmarks/sindex.py
index 7b3a8f13..d48bd299 100644
--- a/benchmarks/sindex.py
+++ b/benchmarks/sindex.py
@@ -1,8 +1,13 @@
-from geopandas import read_file, datasets
-from geopandas.sindex import VALID_QUERY_PREDICATES
+from shapely.geometry import Point
 
+from geopandas import read_file, datasets, GeoSeries
+
+
+# Derive list of valid query predicates based on underlying index backend;
+# we have to create a non-empty instance of the index to get these
+index = GeoSeries([Point(0, 0)]).sindex
+predicates = sorted(p for p in index.valid_query_predicates if p is not None)
 
-predicates = sorted(VALID_QUERY_PREDICATES, key=lambda x: (x is None, x))
 geom_types = ("mixed", "points", "polygons")
 
 
diff --git a/doc/environment.yml b/doc/environment.yml
index 584cbef9..11b55595 100644
--- a/doc/environment.yml
+++ b/doc/environment.yml
@@ -39,4 +39,4 @@ dependencies:
   - pygeos=0.9
   - pip
   - pip:
-      - sphinx-toggleprompt
\ No newline at end of file
+      - sphinx-toggleprompt
diff --git a/doc/source/docs/reference.rst b/doc/source/docs/reference.rst
index 473df638..66bda929 100644
--- a/doc/source/docs/reference.rst
+++ b/doc/source/docs/reference.rst
@@ -15,5 +15,6 @@ The API Reference provides an overview of all public objects, functions and meth
   GeoDataFrame <reference/geodataframe>
   Input/output <reference/io>
   Tools <reference/tools>
+  Spatial index <reference/sindex>
   Testing <reference/testing>
 
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
new file mode 100644
index 00000000..ffed1c98
--- /dev/null
+++ b/doc/source/docs/reference/sindex.rst
@@ -0,0 +1,50 @@
+=============
+Spatial index
+=============
+.. currentmodule:: geopandas
+
+GeoPandas offers built-in support for spatial indexing using an R-Tree algorithm.
+Depending on the ability to import ``pygeos``, GeoPandas will either use
+``pygeos.STRtree`` or ``rtree.index.Index``. The main interface for both is the
+same and follows the ``pygeos`` model.
+
+``GeoSeries.sindex`` creates a spatial index, which can use the methods and
+properties documented below.
+
+Constructor
+-----------
+.. autosummary::
+   :toctree: api/
+
+    GeoSeries.sindex
+
+PyGEOS STRtree
+--------------
+.. currentmodule:: geopandas.sindex.PyGEOSSTRTreeIndex
+.. autosummary::
+   :toctree: api/
+
+    intersection
+    is_empty
+    query
+    query_bulk
+    size
+    valid_query_predicates
+
+rtree Rtree
+-----------
+.. currentmodule:: geopandas.sindex.RTreeIndex
+.. autosummary::
+   :toctree: api/
+
+    intersection
+    is_empty
+    query
+    query_bulk
+    size
+    valid_query_predicates
+
+Furthermore, the ``rtree``-based spatial index offers full capability of
+``rtree.index.Index`` - see the full API in the `rtree documentation`_.
+
+.. _rtree documentation: https://rtree.readthedocs.io/en/stable/class.html
diff --git a/doc/source/docs/user_guide.rst b/doc/source/docs/user_guide.rst
index 02b0210e..1086f58b 100644
--- a/doc/source/docs/user_guide.rst
+++ b/doc/source/docs/user_guide.rst
@@ -13,7 +13,7 @@ Advanced topics can be found in the :doc:`Advanced Guide <advanced_guide>` and f
   Data Structures <user_guide/data_structures>
   Reading and Writing Files <user_guide/io>
   Indexing and Selecting Data <user_guide/indexing>
-  Making Maps <user_guide/mapping>
+  Making Maps and plots <user_guide/mapping>
   Managing Projections <user_guide/projections>
   Geometric Manipulations <user_guide/geometric_manipulations>
   Set Operations with overlay <user_guide/set_operations>
diff --git a/doc/source/docs/user_guide/mapping.rst b/doc/source/docs/user_guide/mapping.rst
index dbad330f..f2a96b5a 100644
--- a/doc/source/docs/user_guide/mapping.rst
+++ b/doc/source/docs/user_guide/mapping.rst
@@ -11,7 +11,7 @@
    plt.close('all')
 
 
-Mapping Tools
+Mapping and Plotting Tools
 =========================================
 
 
@@ -222,6 +222,41 @@ We can set the ``zorder`` for cities higher than for world to move it of top.
     @savefig zorder_set.png
     world.plot(ax=ax, zorder=1);
 
+
+Pandas Plots
+-----------------
+
+Plotting methods also allow for different plot styles from pandas
+along with the default ``geo`` plot. These methods can be accessed using
+the ``kind`` keyword argument in :meth:`~GeoDataFrame.plot`, and include:
+
+* ``geo`` for mapping
+* ``line`` for line plots
+* ``bar`` or ``barh`` for bar plots
+* ``hist`` for histogram
+* ``box`` for boxplot
+* ``kde`` or ``density`` for density plots
+* ``area``  for area plots
+* ``scatter`` for scatter plots
+* ``hexbin`` for hexagonal bin plots
+* ``pie`` for pie plots
+
+.. ipython:: python
+
+    gdf = world.head(10)
+    @savefig pandas_line_plot.png
+    gdf.plot(kind='scatter', x="pop_est", y="gdp_md_est")
+
+You can also create these other plots using the ``GeoDataFrame.plot.<kind>`` accessor methods instead of providing the ``kind`` keyword argument.
+
+.. ipython:: python
+
+    @savefig pandas_bar_plot.png
+    gdf.plot.bar()
+
+For more information check out the `pandas documentation <https://pandas.pydata.org/pandas-docs/stable/user_guide/visualization.html>`_.
+
+
 Other Resources
 -----------------
 Links to jupyter Notebooks for different mapping tasks:
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index 8690d90e..635aeccb 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -5,6 +5,7 @@ import os
 import warnings
 
 import pandas as pd
+import pyproj
 import shapely
 import shapely.geos
 
@@ -16,6 +17,7 @@ import shapely.geos
 PANDAS_GE_025 = str(pd.__version__) >= LooseVersion("0.25.0")
 PANDAS_GE_10 = str(pd.__version__) >= LooseVersion("1.0.0")
 PANDAS_GE_11 = str(pd.__version__) >= LooseVersion("1.1.0")
+PANDAS_GE_115 = str(pd.__version__) >= LooseVersion("1.1.5")
 PANDAS_GE_12 = str(pd.__version__) >= LooseVersion("1.2.0")
 
 
@@ -202,3 +204,9 @@ try:
     HAS_RTREE = True
 except ImportError:
     HAS_RTREE = False
+
+# -----------------------------------------------------------------------------
+# pyproj compat
+# -----------------------------------------------------------------------------
+
+PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
diff --git a/geopandas/array.py b/geopandas/array.py
index bc1e302e..a25e01ba 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -18,7 +18,7 @@ import shapely.geometry
 from shapely.geometry.base import BaseGeometry
 import shapely.ops
 import shapely.wkt
-from pyproj import CRS
+from pyproj import CRS, Transformer
 
 try:
     import pygeos
@@ -708,6 +708,168 @@ class GeometryArray(ExtensionArray):
             crs=self.crs,
         )
 
+    def to_crs(self, crs=None, epsg=None):
+        """Returns a ``GeometryArray`` with all geometries transformed to a new
+        coordinate reference system.
+
+        Transform all geometries in a GeometryArray to a different coordinate
+        reference system.  The ``crs`` attribute on the current GeometryArray must
+        be set.  Either ``crs`` or ``epsg`` may be specified for output.
+
+        This method will transform all points in all objects.  It has no notion
+        or projecting entire geometries.  All segments joining points are
+        assumed to be lines in the current projection, not geodesics.  Objects
+        crossing the dateline (or other projection boundary) will have
+        undesirable behavior.
+
+        Parameters
+        ----------
+        crs : pyproj.CRS, optional if `epsg` is specified
+            The value can be anything accepted
+            by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        epsg : int, optional if `crs` is specified
+            EPSG code specifying output projection.
+
+        Returns
+        -------
+        GeometryArray
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> from geopandas.array import from_shapely, to_wkt
+        >>> a = from_shapely([Point(1, 1), Point(2, 2), Point(3, 3)], crs=4326)
+        >>> to_wkt(a)
+        array(['POINT (1 1)', 'POINT (2 2)', 'POINT (3 3)'], dtype=object)
+        >>> a.crs  # doctest: +SKIP
+        <Geographic 2D CRS: EPSG:4326>
+        Name: WGS 84
+        Axis Info [ellipsoidal]:
+        - Lat[north]: Geodetic latitude (degree)
+        - Lon[east]: Geodetic longitude (degree)
+        Area of Use:
+        - name: World
+        - bounds: (-180.0, -90.0, 180.0, 90.0)
+        Datum: World Geodetic System 1984
+        - Ellipsoid: WGS 84
+        - Prime Meridian: Greenwich
+
+        >>> a = a.to_crs(3857)
+        >>> to_wkt(a)
+        array(['POINT (111319 111325)', 'POINT (222639 222684)',
+               'POINT (333958 334111)'], dtype=object)
+        >>> a.crs  # doctest: +SKIP
+        <Projected CRS: EPSG:3857>
+        Name: WGS 84 / Pseudo-Mercator
+        Axis Info [cartesian]:
+        - X[east]: Easting (metre)
+        - Y[north]: Northing (metre)
+        Area of Use:
+        - name: World - 85°S to 85°N
+        - bounds: (-180.0, -85.06, 180.0, 85.06)
+        Coordinate Operation:
+        - name: Popular Visualisation Pseudo-Mercator
+        - method: Popular Visualisation Pseudo Mercator
+        Datum: World Geodetic System 1984
+        - Ellipsoid: WGS 84
+        - Prime Meridian: Greenwich
+
+        """
+        if self.crs is None:
+            raise ValueError(
+                "Cannot transform naive geometries.  "
+                "Please set a crs on the object first."
+            )
+        if crs is not None:
+            crs = CRS.from_user_input(crs)
+        elif epsg is not None:
+            crs = CRS.from_epsg(epsg)
+        else:
+            raise ValueError("Must pass either crs or epsg.")
+
+        # skip if the input CRS and output CRS are the exact same
+        if self.crs.is_exact_same(crs):
+            return self
+
+        transformer = Transformer.from_crs(self.crs, crs, always_xy=True)
+
+        new_data = vectorized.transform(self.data, transformer.transform)
+        return GeometryArray(new_data, crs=crs)
+
+    def estimate_utm_crs(self, datum_name="WGS 84"):
+        """Returns the estimated UTM CRS based on the bounds of the dataset.
+
+        .. versionadded:: 0.9
+
+        .. note:: Requires pyproj 3+
+
+        Parameters
+        ----------
+        datum_name : str, optional
+            The name of the datum to use in the query. Default is WGS 84.
+
+        Returns
+        -------
+        pyproj.CRS
+
+        Examples
+        --------
+        >>> world = geopandas.read_file(
+        ...     geopandas.datasets.get_path("naturalearth_lowres")
+        ... )
+        >>> germany = world.loc[world.name == "Germany"]
+        >>> germany.geometry.values.estimate_utm_crs()  # doctest: +SKIP
+        <Projected CRS: EPSG:32632>
+        Name: WGS 84 / UTM zone 32N
+        Axis Info [cartesian]:
+        - E[east]: Easting (metre)
+        - N[north]: Northing (metre)
+        Area of Use:
+        - name: World - N hemisphere - 6°E to 12°E - by country
+        - bounds: (6.0, 0.0, 12.0, 84.0)
+        Coordinate Operation:
+        - name: UTM zone 32N
+        - method: Transverse Mercator
+        Datum: World Geodetic System 1984
+        - Ellipsoid: WGS 84
+        - Prime Meridian: Greenwich
+        """
+        try:
+            from pyproj.aoi import AreaOfInterest
+            from pyproj.database import query_utm_crs_info
+        except ImportError:
+            raise RuntimeError("pyproj 3+ required for estimate_utm_crs.")
+
+        if not self.crs:
+            raise RuntimeError("crs must be set to estimate UTM CRS.")
+
+        minx, miny, maxx, maxy = self.total_bounds
+        # ensure using geographic coordinates
+        if not self.crs.is_geographic:
+            lon, lat = Transformer.from_crs(
+                self.crs, "EPSG:4326", always_xy=True
+            ).transform((minx, maxx, minx, maxx), (miny, miny, maxy, maxy))
+            x_center = np.mean(lon)
+            y_center = np.mean(lat)
+        else:
+            x_center = np.mean([minx, maxx])
+            y_center = np.mean([miny, maxy])
+
+        utm_crs_list = query_utm_crs_info(
+            datum_name=datum_name,
+            area_of_interest=AreaOfInterest(
+                west_lon_degree=x_center,
+                south_lat_degree=y_center,
+                east_lon_degree=x_center,
+                north_lat_degree=y_center,
+            ),
+        )
+        try:
+            return CRS.from_epsg(utm_crs_list[0].code)
+        except IndexError:
+            raise RuntimeError("Unable to determine UTM CRS")
+
     #
     # Coordinate related properties
     #
@@ -915,6 +1077,41 @@ class GeometryArray(ExtensionArray):
     def nbytes(self):
         return self.data.nbytes
 
+    def shift(self, periods=1, fill_value=None):
+        """
+        Shift values by desired number.
+
+        Newly introduced missing values are filled with
+        ``self.dtype.na_value``.
+
+        Parameters
+        ----------
+        periods : int, default 1
+            The number of periods to shift. Negative values are allowed
+            for shifting backwards.
+
+        fill_value : object, optional (default None)
+            The scalar value to use for newly introduced missing values.
+            The default is ``self.dtype.na_value``.
+
+        Returns
+        -------
+        GeometryArray
+            Shifted.
+
+        Notes
+        -----
+        If ``self`` is empty or ``periods`` is 0, a copy of ``self`` is
+        returned.
+
+        If ``periods > len(self)``, then an array of size
+        len(self) is returned, with all values filled with
+        ``self.dtype.na_value``.
+        """
+        shifted = super(GeometryArray, self).shift(periods, fill_value)
+        shifted.crs = self.crs
+        return shifted
+
     # -------------------------------------------------------------------------
     # ExtensionArray specific
     # -------------------------------------------------------------------------
@@ -980,7 +1177,7 @@ class GeometryArray(ExtensionArray):
         pandas.factorize
         ExtensionArray.factorize
         """
-        return from_wkb(values)
+        return from_wkb(values, crs=original.crs)
 
     def _values_for_argsort(self):
         # type: () -> np.ndarray
diff --git a/geopandas/base.py b/geopandas/base.py
index a4edd738..878fa533 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2578,6 +2578,55 @@ GeometryCollection
 
     @property
     def sindex(self):
+        """Generate the spatial index
+
+        Creates R-tree spatial index based on ``pygeos.STRtree`` or
+        ``rtree.index.Index``.
+
+        Note that the  spatial index may not be fully
+        initialized until the first use.
+
+        Examples
+        --------
+        >>> from shapely.geometry import box
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(5), range(5)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        dtype: geometry
+
+        Query the spatial index with a single geometry based on the bounding box:
+
+        >>> s.sindex.query(box(1, 1, 3, 3))
+        array([1, 2, 3])
+
+        Query the spatial index with a single geometry based on the predicate:
+
+        >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
+        array([2])
+
+        Query the spatial index with an array of geometries based on the bounding
+        box:
+
+        >>> s2 = geopandas.GeoSeries([box(1, 1, 3, 3), box(4, 4, 5, 5)])
+        >>> s2
+        0    POLYGON ((3.00000 1.00000, 3.00000 3.00000, 1....
+        1    POLYGON ((5.00000 4.00000, 5.00000 5.00000, 4....
+        dtype: geometry
+
+        >>> s.sindex.query_bulk(s2)
+        array([[0, 0, 0, 1],
+               [1, 2, 3, 4]])
+
+        Query the spatial index with an array of geometries based on the predicate:
+
+        >>> s.sindex.query_bulk(s2, predicate="contains")
+        array([[0],
+               [2]])
+        """
         return self.geometry.values.sindex
 
     @property
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 78fc3a34..d3d18d77 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1312,20 +1312,6 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
 
         return self
 
-    def plot(self, *args, **kwargs):
-        """Generate a plot of the geometries in the ``GeoDataFrame``.
-
-        If the ``column`` parameter is given, colors plot according to values
-        in that column, otherwise calls ``GeoSeries.plot()`` on the
-        ``geometry`` column.
-
-        Wraps the ``plot_dataframe()`` function, and documentation is copied
-        from there.
-        """
-        return plot_dataframe(self, *args, **kwargs)
-
-    plot.__doc__ = plot_dataframe.__doc__
-
     def dissolve(self, by=None, aggfunc="first", as_index=True):
         """
         Dissolve geometries within `groupby` into single observation.
@@ -1617,6 +1603,24 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
         )
         return self.geometry.difference(other)
 
+    if compat.PANDAS_GE_025:
+        from pandas.core.accessor import CachedAccessor
+
+        plot = CachedAccessor("plot", geopandas.plotting.GeoplotAccessor)
+    else:
+
+        def plot(self, *args, **kwargs):
+            """Generate a plot of the geometries in the ``GeoDataFrame``.
+            If the ``column`` parameter is given, colors plot according to values
+            in that column, otherwise calls ``GeoSeries.plot()`` on the
+            ``geometry`` column.
+            Wraps the ``plot_dataframe()`` function, and documentation is copied
+            from there.
+            """
+            return plot_dataframe(self, *args, **kwargs)
+
+    plot.__doc__ = plot_dataframe.__doc__
+
 
 def _dataframe_set_geometry(self, col, drop=False, inplace=False, crs=None):
     if inplace:
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 82c6bd92..8444d978 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -6,15 +6,14 @@ import pandas as pd
 from pandas import Series, MultiIndex
 from pandas.core.internals import SingleBlockManager
 
-from pyproj import CRS, Transformer
+from pyproj import CRS
 from shapely.geometry.base import BaseGeometry
 
 from geopandas.base import GeoPandasBase, _delegate_property
 from geopandas.plotting import plot_series
 
-from .array import GeometryArray, GeometryDtype, from_shapely
+from .array import GeometryDtype, from_shapely
 from .base import is_geometry_type
-from . import _vectorized as vectorized
 from . import _compat as compat
 
 
@@ -201,7 +200,7 @@ class GeoSeries(GeoPandasBase, Series):
                 s = pd.Series(data, index=index, name=name, **kwargs)
             # prevent trying to convert non-geometry objects
             if s.dtype != object:
-                if s.empty:
+                if s.empty or data is None:
                     s = s.astype(object)
                 else:
                     warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
@@ -927,27 +926,8 @@ class GeoSeries(GeoPandasBase, Series):
         GeoSeries.set_crs : assign CRS
 
         """
-        if self.crs is None:
-            raise ValueError(
-                "Cannot transform naive geometries.  "
-                "Please set a crs on the object first."
-            )
-        if crs is not None:
-            crs = CRS.from_user_input(crs)
-        elif epsg is not None:
-            crs = CRS.from_epsg(epsg)
-        else:
-            raise ValueError("Must pass either crs or epsg.")
-
-        # skip if the input CRS and output CRS are the exact same
-        if self.crs.is_exact_same(crs):
-            return self
-
-        transformer = Transformer.from_crs(self.crs, crs, always_xy=True)
-
-        new_data = vectorized.transform(self.values.data, transformer.transform)
         return GeoSeries(
-            GeometryArray(new_data), crs=crs, index=self.index, name=self.name
+            self.values.to_crs(crs=crs, epsg=epsg), index=self.index, name=self.name
         )
 
     def estimate_utm_crs(self, datum_name="WGS 84"):
@@ -988,40 +968,7 @@ class GeoSeries(GeoPandasBase, Series):
         - Ellipsoid: WGS 84
         - Prime Meridian: Greenwich
         """
-        try:
-            from pyproj.aoi import AreaOfInterest
-            from pyproj.database import query_utm_crs_info
-        except ImportError:
-            raise RuntimeError("pyproj 3+ required for estimate_utm_crs.")
-
-        if not self.crs:
-            raise RuntimeError("crs must be set to estimate UTM CRS.")
-
-        minx, miny, maxx, maxy = self.total_bounds
-        # ensure using geographic coordinates
-        if not self.crs.is_geographic:
-            lon, lat = Transformer.from_crs(
-                self.crs, "EPSG:4326", always_xy=True
-            ).transform((minx, maxx, minx, maxx), (miny, miny, maxy, maxy))
-            x_center = np.mean(lon)
-            y_center = np.mean(lat)
-        else:
-            x_center = np.mean([minx, maxx])
-            y_center = np.mean([miny, maxy])
-
-        utm_crs_list = query_utm_crs_info(
-            datum_name=datum_name,
-            area_of_interest=AreaOfInterest(
-                west_lon_degree=x_center,
-                south_lat_degree=y_center,
-                east_lon_degree=x_center,
-                north_lat_degree=y_center,
-            ),
-        )
-        try:
-            return CRS.from_epsg(utm_crs_list[0].code)
-        except IndexError:
-            raise RuntimeError("Unable to determine UTM CRS")
+        return self.values.estimate_utm_crs(datum_name)
 
     def to_json(self, **kwargs):
         """
diff --git a/geopandas/io/sql.py b/geopandas/io/sql.py
index 250478ef..4aa3767d 100644
--- a/geopandas/io/sql.py
+++ b/geopandas/io/sql.py
@@ -306,7 +306,9 @@ def _psql_insert_copy(tbl, conn, keys, data_iter):
 
     dbapi_conn = conn.connection
     with dbapi_conn.cursor() as cur:
-        sql = "COPY {} ({}) FROM STDIN WITH CSV".format(tbl.table.fullname, columns)
+        sql = 'COPY "{}"."{}" ({}) FROM STDIN WITH CSV'.format(
+            tbl.table.schema, tbl.table.name, columns
+        )
         cur.copy_expert(sql=sql, file=s_buf)
 
 
@@ -399,14 +401,14 @@ def _write_postgis(
     # Convert geometries to EWKB
     gdf = _convert_to_ewkb(gdf, geom_name, srid)
 
+    if schema is not None:
+        schema_name = schema
+    else:
+        schema_name = "public"
+
     if if_exists == "append":
         # Check that the geometry srid matches with the current GeoDataFrame
         with _get_conn(con) as connection:
-            if schema is not None:
-                schema_name = schema
-            else:
-                schema_name = "public"
-
             # Only check SRID if table exists
             if connection.dialect.has_table(connection, name, schema):
                 target_srid = connection.execute(
@@ -429,7 +431,7 @@ def _write_postgis(
         gdf.to_sql(
             name,
             connection,
-            schema=schema,
+            schema=schema_name,
             if_exists=if_exists,
             index=index,
             index_label=index_label,
diff --git a/geopandas/plotting.py b/geopandas/plotting.py
index 9c395685..7159aa2d 100644
--- a/geopandas/plotting.py
+++ b/geopandas/plotting.py
@@ -115,7 +115,7 @@ def _PolygonPatch(polygon, **kwargs):
 
     path = Path.make_compound_path(
         Path(np.asarray(polygon.exterior.coords)[:, :2]),
-        *[Path(np.asarray(ring.coords)[:, :2]) for ring in polygon.interiors]
+        *[Path(np.asarray(ring.coords)[:, :2]) for ring in polygon.interiors],
     )
     return PathPatch(path, **kwargs)
 
@@ -254,7 +254,7 @@ def _plot_point_collection(
     vmax=None,
     marker="o",
     markersize=None,
-    **kwargs
+    **kwargs,
 ):
     """
     Plots a collection of Point and MultiPoint geometries to `ax`
@@ -488,7 +488,7 @@ def plot_dataframe(
     classification_kwds=None,
     missing_kwds=None,
     aspect="auto",
-    **style_kwds
+    **style_kwds,
 ):
     """
     Plot a GeoDataFrame.
@@ -508,6 +508,21 @@ def plot_dataframe(
         If np.array or pd.Series are used then it must have same length as
         dataframe. Values are used to color the plot. Ignored if `color` is
         also set.
+    kind: str
+        The kind of plots to produce:
+         - 'geo': Map (default)
+         Pandas Kinds
+         - 'line' : line plot
+         - 'bar' : vertical bar plot
+         - 'barh' : horizontal bar plot
+         - 'hist' : histogram
+         - 'box' : BoxPlot
+         - 'kde' : Kernel Density Estimation plot
+         - 'density' : same as 'kde'
+         - 'area' : area plot
+         - 'pie' : pie plot
+         - 'scatter' : scatter plot
+         - 'hexbin' : hexbin plot.
     cmap : str (default None)
         The name of a colormap recognized by matplotlib.
     color : str (default None)
@@ -683,7 +698,7 @@ GON (((-122.84000 49.00000, -120.0000...
             figsize=figsize,
             markersize=markersize,
             aspect=aspect,
-            **style_kwds
+            **style_kwds,
         )
 
     # To accept pd.Series and np.arrays as column
@@ -820,7 +835,7 @@ GON (((-122.84000 49.00000, -120.0000...
             vmax=mx,
             markersize=markersize,
             cmap=cmap,
-            **style_kwds
+            **style_kwds,
         )
 
     if missing_kwds is not None and not expl_series[nan_idx].empty:
@@ -899,6 +914,30 @@ GON (((-122.84000 49.00000, -120.0000...
     return ax
 
 
+if geopandas._compat.PANDAS_GE_025:
+    from pandas.plotting import PlotAccessor
+
+    class GeoplotAccessor(PlotAccessor):
+
+        __doc__ = plot_dataframe.__doc__
+        _pandas_kinds = PlotAccessor._all_kinds
+
+        def __call__(self, *args, **kwargs):
+            data = self._parent.copy()
+            kind = kwargs.pop("kind", "geo")
+            if kind == "geo":
+                return plot_dataframe(data, *args, **kwargs)
+            if kind in self._pandas_kinds:
+                # Access pandas plots
+                return PlotAccessor(data)(kind=kind, **kwargs)
+            else:
+                # raise error
+                raise ValueError(f"{kind} is not a valid plot kind")
+
+        def geo(self, *args, **kwargs):
+            return self(kind="geo", *args, **kwargs)
+
+
 def _mapclassify_choro(values, scheme, **classification_kwds):
     """
     Wrapper for choropleth schemes from mapclassify for use with plot_dataframe
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 0716bc98..03eab7c9 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -82,6 +82,14 @@ if compat.HAS_RTREE:
             -------
             set
                 Set of valid predicates for this spatial index.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries([Point(0, 0), Point(1, 1)])
+            >>> s.sindex.valid_query_predicates  # doctest: +SKIP
+            {'contains', 'crosses', 'intersects', 'within', 'touches', \
+'overlaps', None, 'covers', 'contains_properly'}
             """
             return {
                 None,
@@ -96,7 +104,10 @@ if compat.HAS_RTREE:
             }
 
         def query(self, geometry, predicate=None, sort=False):
-            """Compatibility layer for pygeos.query.
+            """Return the index of all geometries in the tree with extents that
+            intersect the envelope of the input geometry.
+
+            Compatibility layer for pygeos-based ``sindex.query``.
 
             This is not a vectorized function, if speed is important,
             please use PyGEOS.
@@ -121,6 +132,29 @@ if compat.HAS_RTREE:
             -------
             matches : ndarray of shape (n_results, )
                 Integer indices for matching geometries from the spatial index.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.query(box(1, 1, 3, 3))
+            array([1, 2, 3])
+
+            >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
+            array([2])
             """
 
             # handle invalid predicates
@@ -202,11 +236,16 @@ if compat.HAS_RTREE:
             return np.array(tree_idx, dtype=np.intp)
 
         def query_bulk(self, geometry, predicate=None, sort=False):
-            """Compatibility layer for pygeos.query_bulk.
+            """
+            Returns all combinations of each input geometry and geometries in
+            the tree where the envelope of each input geometry intersects with
+            the envelope of a tree geometry.
 
-            Iterates over `geometry` and queries index.
+            Compatibility layer for pygeos-based ``sindex.query_bulk``.
+
+            Iterates over ``geometry`` and queries index.
             This operation is not vectorized and may be slow.
-            Use PyGEOS with `query_bulk` for speed.
+            Use PyGEOS with ``query_bulk`` for speed.
 
             Parameters
             ----------
@@ -230,6 +269,36 @@ if compat.HAS_RTREE:
             ndarray with shape (2, n)
                 The first subarray contains input geometry integer indexes.
                 The second subarray contains tree geometry integer indexes.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+            >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
+            >>> s2
+            0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
+            1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
+            dtype: geometry
+
+            >>> s.sindex.query_bulk(s2)
+            array([[0, 0, 0, 1, 1],
+                   [2, 3, 4, 5, 6]])
+
+            >>> s.sindex.query_bulk(s2, predicate="contains")
+            array([[0],
+                   [3]])
             """
             # Iterates over geometry, applying func.
             tree_index = []
@@ -249,11 +318,61 @@ if compat.HAS_RTREE:
             coordinates : sequence or array
                 Sequence of the form (min_x, min_y, max_x, max_y)
                 to query a rectangle or (x, y) to query a point.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.intersection(box(1, 1, 3, 3).bounds)
+            array([1, 2, 3])
+
+            Alternatively, you can use ``query``:
+
+            >>> s.sindex.query(box(1, 1, 3, 3))
+            array([1, 2, 3])
+
             """
             return super().intersection(coordinates, objects=False)
 
         @property
         def size(self):
+            """Size of the spatial index
+
+            Number of leaves (input geometries) in the index.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.size
+            10
+            """
             if hasattr(self, "_size"):
                 size = self._size
             else:
@@ -266,6 +385,32 @@ if compat.HAS_RTREE:
 
         @property
         def is_empty(self):
+            """Check if the spatial index is empty
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.is_empty
+            False
+
+            >>> s2 = geopandas.GeoSeries()
+            >>> s2.sindex.is_empty
+            True
+            """
             return self.geometries.size == 0 or self.size == 0
 
         def __len__(self):
@@ -302,17 +447,29 @@ if compat.HAS_PYGEOS:
 
         @property
         def valid_query_predicates(self):
-            """Returns valid predicates for this spatial index.
+            """Returns valid predicates for the used spatial index.
 
             Returns
             -------
             set
                 Set of valid predicates for this spatial index.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries([Point(0, 0), Point(1, 1)])
+            >>> s.sindex.valid_query_predicates  # doctest: +SKIP
+            {'contains', 'crosses', 'covered_by', None, 'intersects', 'within', \
+'touches', 'overlaps', 'contains_properly', 'covers'}
             """
             return pygeos.strtree.VALID_PREDICATES | set([None])
 
         def query(self, geometry, predicate=None, sort=False):
-            """Wrapper for pygeos.query.
+            """
+            Return the index of all geometries in the tree with extents
+            that intersect the envelope of the input geometry.
+
+            Wrapper for pygeos.STRtree.query.
 
             This also ensures a deterministic (sorted) order for the results.
 
@@ -334,9 +491,32 @@ if compat.HAS_PYGEOS:
             matches : ndarray of shape (n_results, )
                 Integer indices for matching geometries from the spatial index.
 
-            See also
-            --------
+            Notes
+            -----
             See PyGEOS.strtree documentation for more information.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.query(box(1, 1, 3, 3))
+            array([1, 2, 3])
+
+            >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
+            array([2])
             """
 
             if predicate not in self.valid_query_predicates:
@@ -358,7 +538,19 @@ if compat.HAS_PYGEOS:
             return matches
 
         def query_bulk(self, geometry, predicate=None, sort=False):
-            """Wrapper to expose underlaying pygeos objects to pygeos.query_bulk.
+            """
+            Returns all combinations of each input geometry and geometries in
+            the tree where the envelope of each input geometry intersects with
+            the envelope of a tree geometry.
+
+            Wrapper to expose underlaying pygeos objects to pygeos.query_bulk.
+
+            In the context of a spatial join, input geometries are the “left”
+            geometries that determine the order of the results, and tree geometries
+            are “right” geometries that are joined against the left geometries.
+            This effectively performs an inner join, where only those combinations
+            of geometries that can be joined based on envelope overlap or optional
+            predicate are returned.
 
             This also allows a deterministic (sorted) order for the results.
 
@@ -385,9 +577,39 @@ if compat.HAS_PYGEOS:
                 The first subarray contains input geometry integer indexes.
                 The second subarray contains tree geometry integer indexes.
 
-            See also
-            --------
+            Notes
+            -----
             See PyGEOS.strtree documentation for more information.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+            >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
+            >>> s2
+            0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
+            1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
+            dtype: geometry
+
+            >>> s.sindex.query_bulk(s2)
+            array([[0, 0, 0, 1, 1],
+                   [2, 3, 4, 5, 6]])
+
+            >>> s.sindex.query_bulk(s2, predicate="contains")
+            array([[0],
+                   [3]])
             """
 
             if predicate not in self.valid_query_predicates:
@@ -416,11 +638,37 @@ if compat.HAS_PYGEOS:
         def intersection(self, coordinates):
             """Wrapper for pygeos.query that uses the RTree API.
 
+            Compatibility wrapper, use ``query`` instead.
+
             Parameters
             ----------
             coordinates : sequence or array
                 Sequence of the form (min_x, min_y, max_x, max_y)
                 to query a rectangle or (x, y) to query a point.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point, box
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.query(box(1, 1, 3, 3))
+            array([1, 2, 3])
+
+            >>> s.sindex.intersection(box(1, 1, 3, 3).bounds)
+            array([1, 2, 3])
+
             """
             # convert bounds to geometry
             # the old API uses tuples of bound, but pygeos uses geometries
@@ -452,8 +700,58 @@ if compat.HAS_PYGEOS:
 
         @property
         def size(self):
+            """Size of the spatial index
+
+            Number of leaves (input geometries) in the index.
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.size
+            10
+            """
             return len(self)
 
         @property
         def is_empty(self):
+            """Check if the spatial index is empty
+
+            Examples
+            --------
+            >>> from shapely.geometry import Point
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            3    POINT (3.00000 3.00000)
+            4    POINT (4.00000 4.00000)
+            5    POINT (5.00000 5.00000)
+            6    POINT (6.00000 6.00000)
+            7    POINT (7.00000 7.00000)
+            8    POINT (8.00000 8.00000)
+            9    POINT (9.00000 9.00000)
+            dtype: geometry
+
+            >>> s.sindex.is_empty
+            False
+
+            >>> s2 = geopandas.GeoSeries()
+            >>> s2.sindex.is_empty
+            True
+            """
             return len(self) == 0
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index 9aebbb39..5b95c922 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -137,7 +137,7 @@ def _overlay_union(df1, df2):
     return dfunion.reindex(columns=columns)
 
 
-def overlay(df1, df2, how="intersection", keep_geom_type=None):
+def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     """Perform spatial overlay between two GeoDataFrames.
 
     Currently only supports data GeoDataFrames with uniform geometry types,
@@ -159,6 +159,9 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
         if False, return all resulting geometries. Default is None,
         which will set keep_geom_type to True but warn upon dropping
         geometries.
+    make_valid : bool, default True
+        If True, any invalid input geometries are corrected with a call to `buffer(0)`,
+        if False, a `ValueError` is raised if any input geometries are invalid.
 
     Returns
     -------
@@ -179,9 +182,9 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
     >>> geopandas.overlay(df1, df2, how='union')
        df1_data  df2_data                                           geometry
     0       1.0       1.0  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1       2.0       1.0  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2       2.0       2.0  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
-    3       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    3       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
     5       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     6       NaN       2.0  POLYGON ((3.00000 4.00000, 3.00000 5.00000, 5....
@@ -189,27 +192,27 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
     >>> geopandas.overlay(df1, df2, how='intersection')
        df1_data  df2_data                                           geometry
     0         1         1  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1         2         1  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1         2         1  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2         2         2  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
 
     >>> geopandas.overlay(df1, df2, how='symmetric_difference')
        df1_data  df2_data                                           geometry
-    0       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    0       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     1       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
     2       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     3       NaN       2.0  POLYGON ((3.00000 4.00000, 3.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='difference')
                                                 geometry  df1_data
-    0  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....         1
+    0  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....         1
     1  MULTIPOLYGON (((2.00000 3.00000, 2.00000 4.000...         2
 
     >>> geopandas.overlay(df1, df2, how='identity')
        df1_data  df2_data                                           geometry
     0       1.0       1.0  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1       2.0       1.0  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2       2.0       2.0  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
-    3       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    3       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
 
     See also
@@ -262,12 +265,24 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
             )
 
     # Computations
-    df1 = df1.copy()
-    df2 = df2.copy()
-    if df1.geom_type.isin(polys).all():
-        df1[df1._geometry_column_name] = df1.geometry.buffer(0)
-    if df2.geom_type.isin(polys).all():
-        df2[df2._geometry_column_name] = df2.geometry.buffer(0)
+    def _make_valid(df):
+        df = df.copy()
+        if df.geom_type.isin(polys).all():
+            mask = ~df.geometry.is_valid
+            col = df._geometry_column_name
+            if make_valid:
+                df.loc[mask, col] = df.loc[mask, col].buffer(0)
+            elif mask.any():
+                raise ValueError(
+                    "You have passed make_valid=False along with "
+                    f"{mask.sum()} invalid input geometries. "
+                    "Use make_valid=True or make sure that all geometries "
+                    "are valid before using overlay."
+                )
+        return df
+
+    df1 = _make_valid(df1)
+    df2 = _make_valid(df2)
 
     with warnings.catch_warnings():  # CRS checked above, supress array-level warning
         warnings.filterwarnings("ignore", message="CRS mismatch between the CRS")


Test patch:
diff --git a/ci/envs/37-latest-conda-forge.yaml b/ci/envs/37-latest-conda-forge.yaml
index a0d148ac..1ca12938 100644
--- a/ci/envs/37-latest-conda-forge.yaml
+++ b/ci/envs/37-latest-conda-forge.yaml
@@ -18,6 +18,7 @@ dependencies:
   - rtree
   - matplotlib
   - mapclassify
+  - scipy
   - geopy
   - SQLalchemy
   - libspatialite
diff --git a/ci/envs/38-latest-conda-forge.yaml b/ci/envs/38-latest-conda-forge.yaml
index 5ebdc26a..a9ecb787 100644
--- a/ci/envs/38-latest-conda-forge.yaml
+++ b/ci/envs/38-latest-conda-forge.yaml
@@ -18,6 +18,7 @@ dependencies:
   - rtree
   - matplotlib
   - mapclassify
+  - scipy
   - geopy
   # installed in tests.yaml, because not available on windows
   # - postgis
diff --git a/ci/envs/39-latest-conda-forge.yaml b/ci/envs/39-latest-conda-forge.yaml
index 1f856524..bb8c8f35 100644
--- a/ci/envs/39-latest-conda-forge.yaml
+++ b/ci/envs/39-latest-conda-forge.yaml
@@ -19,6 +19,7 @@ dependencies:
   - matplotlib
   - descartes
   - mapclassify
+  - scipy
   - geopy
   # installed in tests.yaml, because not available on windows
   # - postgis
diff --git a/geopandas/io/tests/test_sql.py b/geopandas/io/tests/test_sql.py
index c9e64659..90e5c258 100644
--- a/geopandas/io/tests/test_sql.py
+++ b/geopandas/io/tests/test_sql.py
@@ -344,6 +344,21 @@ class TestIO:
         df = read_postgis(sql, engine, geom_col="geometry")
         validate_boro_df(df)
 
+    def test_write_postgis_uppercase_tablename(self, engine_postgis, df_nybb):
+        """Tests writing GeoDataFrame to PostGIS with uppercase tablename."""
+        engine = engine_postgis
+        table = "aTestTable"
+
+        # If table exists, delete it before trying to write with defaults
+        drop_table_if_exists(engine, table)
+
+        # Write to db
+        write_postgis(df_nybb, con=engine, name=table, if_exists="fail")
+        # Validate
+        sql = 'SELECT * FROM "{table}";'.format(table=table)
+        df = read_postgis(sql, engine, geom_col="geometry")
+        validate_boro_df(df)
+
     def test_write_postgis_sqlalchemy_connection(self, engine_postgis, df_nybb):
         """Tests that GeoDataFrame can be written to PostGIS with defaults."""
         with engine_postgis.begin() as con:
diff --git a/geopandas/testing.py b/geopandas/testing.py
index 437af0d3..7e0daa71 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -133,15 +133,58 @@ def assert_geoseries_equal(
     if not check_crs:
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", "CRS mismatch", UserWarning)
-            if check_less_precise:
-                assert geom_almost_equals(left, right)
-            else:
-                assert geom_equals(left, right)
+            _check_equality(left, right, check_less_precise)
+    else:
+        _check_equality(left, right, check_less_precise)
+
+
+def _truncated_string(geom):
+    """Truncated WKT repr of geom"""
+    s = str(geom)
+    if len(s) > 100:
+        return s[:100] + "..."
+    else:
+        return s
+
+
+def _check_equality(left, right, check_less_precise):
+    assert_error_message = (
+        "{0} out of {1} geometries are not {3}equal.\n"
+        "Indices where geometries are not {3}equal: {2} \n"
+        "The first not {3}equal geometry:\n"
+        "Left: {4}\n"
+        "Right: {5}\n"
+    )
+    if check_less_precise:
+        precise = "almost "
+        if not geom_almost_equals(left, right):
+            unequal_left_geoms = left[~left.geom_almost_equals(right)]
+            unequal_right_geoms = right[~left.geom_almost_equals(right)]
+            raise AssertionError(
+                assert_error_message.format(
+                    len(unequal_left_geoms),
+                    len(left),
+                    unequal_left_geoms.index.to_list(),
+                    precise,
+                    _truncated_string(unequal_left_geoms.iloc[0]),
+                    _truncated_string(unequal_right_geoms.iloc[0]),
+                )
+            )
     else:
-        if check_less_precise:
-            assert geom_almost_equals(left, right)
-        else:
-            assert geom_equals(left, right)
+        precise = ""
+        if not geom_equals(left, right):
+            unequal_left_geoms = left[~left.geom_almost_equals(right)]
+            unequal_right_geoms = right[~left.geom_almost_equals(right)]
+            raise AssertionError(
+                assert_error_message.format(
+                    len(unequal_left_geoms),
+                    len(left),
+                    unequal_left_geoms.index.to_list(),
+                    precise,
+                    _truncated_string(unequal_left_geoms.iloc[0]),
+                    _truncated_string(unequal_right_geoms.iloc[0]),
+                )
+            )
 
 
 def assert_geodataframe_equal(
@@ -212,10 +255,7 @@ def assert_geodataframe_equal(
         "GeoDataFrame shape mismatch, left: {lshape!r}, right: {rshape!r}.\n"
         "Left columns: {lcols!r}, right columns: {rcols!r}"
     ).format(
-        lshape=left.shape,
-        rshape=right.shape,
-        lcols=left.columns,
-        rcols=right.columns,
+        lshape=left.shape, rshape=right.shape, lcols=left.columns, rcols=right.columns
     )
 
     if check_like:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 4cc6e687..ef3cd143 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -4,6 +4,7 @@ import numpy as np
 import pandas as pd
 import six
 
+from pyproj import CRS
 import shapely
 import shapely.affinity
 import shapely.geometry
@@ -883,3 +884,55 @@ def test_isna_pdNA():
     t1 = T.copy()
     t1[0] = pd.NA
     assert t1[0] is None
+
+
+def test_shift_has_crs():
+    t = T.copy()
+    t.crs = 4326
+    assert t.shift(1).crs == t.crs
+    assert t.shift(0).crs == t.crs
+    assert t.shift(-1).crs == t.crs
+
+
+@pytest.mark.skipif(
+    not compat.PANDAS_GE_115, reason="crs only preserved in unique after pandas 1.1.5"
+)
+def test_unique_has_crs():
+    t = T.copy()
+    t.crs = 4326
+    assert t.unique().crs == t.crs
+
+
+class TestEstimateUtmCrs:
+    def setup_method(self):
+        self.esb = shapely.geometry.Point(-73.9847, 40.7484)
+        self.sol = shapely.geometry.Point(-74.0446, 40.6893)
+        self.landmarks = from_shapely([self.esb, self.sol], crs="epsg:4326")
+
+    def test_estimate_utm_crs__geographic(self):
+        if compat.PYPROJ_LT_3:
+            with pytest.raises(RuntimeError, match=r"pyproj 3\+ required"):
+                self.landmarks.estimate_utm_crs()
+        else:
+            assert self.landmarks.estimate_utm_crs() == CRS("EPSG:32618")
+            assert self.landmarks.estimate_utm_crs("NAD83") == CRS("EPSG:26918")
+
+    @pytest.mark.skipif(compat.PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__projected(self):
+        assert self.landmarks.to_crs("EPSG:3857").estimate_utm_crs() == CRS(
+            "EPSG:32618"
+        )
+
+    @pytest.mark.skipif(compat.PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__out_of_bounds(self):
+        with pytest.raises(RuntimeError, match="Unable to determine UTM CRS"):
+            from_shapely(
+                [shapely.geometry.Polygon([(0, 90), (1, 90), (2, 90)])], crs="EPSG:4326"
+            ).estimate_utm_crs()
+
+    @pytest.mark.skipif(compat.PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__missing_crs(self):
+        with pytest.raises(RuntimeError, match="crs must be set"):
+            from_shapely(
+                [shapely.geometry.Polygon([(0, 90), (1, 90), (2, 90)])]
+            ).estimate_utm_crs()
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index a804d4b7..14f4f719 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -392,7 +392,7 @@ class TestGeometryArrayCRS:
         assert df.column1.crs == self.wgs
         assert df.column1.values.crs == self.wgs
 
-    def test_to_crs(self):
+    def test_geoseries_to_crs(self):
         s = GeoSeries(self.geoms, crs=27700)
         s = s.to_crs(4326)
         assert s.crs == self.wgs
@@ -412,6 +412,11 @@ class TestGeometryArrayCRS:
         assert df.col1.crs == self.wgs
         assert df.col1.values.crs == self.wgs
 
+    def test_array_to_crs(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        arr = arr.to_crs(4326)
+        assert arr.crs == self.wgs
+
     def test_from_shapely(self):
         arr = from_shapely(self.geoms, crs=27700)
         assert arr.crs == self.osgb
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 1ffb6494..bdb44d54 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -1,4 +1,3 @@
-from distutils.version import LooseVersion
 import json
 import os
 import random
@@ -19,9 +18,9 @@ from shapely.geometry import (
     Polygon,
 )
 from shapely.geometry.base import BaseGeometry
-import pyproj
 
 from geopandas import GeoSeries, GeoDataFrame
+from geopandas._compat import PYPROJ_LT_3
 from geopandas.array import GeometryArray, GeometryDtype
 
 from geopandas.tests.util import geom_equals
@@ -29,9 +28,6 @@ from pandas.testing import assert_series_equal
 import pytest
 
 
-PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
-
-
 class TestSeries:
     def setup_method(self):
         self.tempdir = tempfile.mkdtemp()
@@ -394,6 +390,10 @@ class TestConstructor:
         s = GeoSeries()
         check_geoseries(s)
 
+    def test_data_is_none(self):
+        s = GeoSeries(index=range(3))
+        check_geoseries(s)
+
     def test_from_series(self):
         shapes = [
             Polygon([(random.random(), random.random()) for _ in range(3)])
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index b8281122..f0c095a6 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -2,7 +2,7 @@ import os
 
 import pandas as pd
 
-from shapely.geometry import Point, Polygon, LineString, GeometryCollection
+from shapely.geometry import Point, Polygon, LineString, GeometryCollection, box
 from fiona.errors import DriverError
 
 import geopandas
@@ -571,3 +571,24 @@ def test_keep_geom_type_geometry_collection():
     intersection = overlay(df1, df2, keep_geom_type=False)
     assert len(intersection) == 1
     assert (intersection.geom_type == "GeometryCollection").all()
+
+
+@pytest.mark.parametrize("make_valid", [True, False])
+def test_overlap_make_valid(make_valid):
+    bowtie = Polygon([(1, 1), (9, 9), (9, 1), (1, 9), (1, 1)])
+    assert not bowtie.is_valid
+    fixed_bowtie = bowtie.buffer(0)
+    assert fixed_bowtie.is_valid
+
+    df1 = GeoDataFrame({"col1": ["region"], "geometry": GeoSeries([box(0, 0, 10, 10)])})
+    df_bowtie = GeoDataFrame(
+        {"col1": ["invalid", "valid"], "geometry": GeoSeries([bowtie, fixed_bowtie])}
+    )
+
+    if make_valid:
+        df_overlay_bowtie = overlay(df1, df_bowtie, make_valid=make_valid)
+        assert df_overlay_bowtie.at[0, "geometry"].equals(fixed_bowtie)
+        assert df_overlay_bowtie.at[1, "geometry"].equals(fixed_bowtie)
+    else:
+        with pytest.raises(ValueError, match="1 invalid input geometries"):
+            overlay(df1, df_bowtie, make_valid=make_valid)
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index 2f59958f..884d64ee 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -29,6 +29,13 @@ matplotlib = pytest.importorskip("matplotlib")
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt  # noqa
 
+try:  # skipif and importorskip do not work for decorators
+    from matplotlib.testing.decorators import check_figures_equal
+
+    MPL_DECORATORS = True
+except ImportError:
+    MPL_DECORATORS = False
+
 
 @pytest.fixture(autouse=True)
 def close_figures(request):
@@ -42,6 +49,8 @@ try:
 except KeyError:
     MPL_DFT_COLOR = matplotlib.rcParams["axes.color_cycle"][0]
 
+plt.rcParams.update({"figure.max_open_warning": 0})
+
 
 class TestPointPlotting:
     def setup_method(self):
@@ -1466,6 +1475,79 @@ class TestPlotCollections:
         ax.cla()
 
 
+@pytest.mark.skipif(not compat.PANDAS_GE_025, reason="requires pandas > 0.24")
+class TestGeoplotAccessor:
+    def setup_method(self):
+        geometries = [Polygon([(0, 0), (1, 0), (1, 1)]), Point(1, 3)]
+        x = [1, 2]
+        y = [10, 20]
+        self.gdf = GeoDataFrame({"geometry": geometries, "x": x, "y": y})
+        self.df = pd.DataFrame({"x": x, "y": y})
+
+    def compare_figures(self, kind, fig_test, fig_ref, kwargs):
+        """Compare Figures."""
+        ax_pandas_1 = fig_test.subplots()
+        self.df.plot(kind=kind, ax=ax_pandas_1, **kwargs)
+        ax_geopandas_1 = fig_ref.subplots()
+        self.gdf.plot(kind=kind, ax=ax_geopandas_1, **kwargs)
+
+        ax_pandas_2 = fig_test.subplots()
+        getattr(self.df.plot, kind)(ax=ax_pandas_2, **kwargs)
+        ax_geopandas_2 = fig_ref.subplots()
+        getattr(self.gdf.plot, kind)(ax=ax_geopandas_2, **kwargs)
+
+    _pandas_kinds = []
+    if compat.PANDAS_GE_025:
+        from geopandas.plotting import GeoplotAccessor
+
+        _pandas_kinds = GeoplotAccessor._pandas_kinds
+
+    if MPL_DECORATORS:
+
+        @pytest.mark.parametrize("kind", _pandas_kinds)
+        @check_figures_equal(extensions=["png", "pdf"])
+        def test_pandas_kind(self, kind, fig_test, fig_ref):
+            """Test Pandas kind."""
+            import importlib
+
+            _scipy_dependent_kinds = ["kde", "density"]  # Needs scipy
+            _y_kinds = ["pie"]  # Needs y
+            _xy_kinds = ["scatter", "hexbin"]  # Needs x & y
+            kwargs = {}
+            if kind in _scipy_dependent_kinds:
+                if not importlib.util.find_spec("scipy"):
+                    with pytest.raises(
+                        ModuleNotFoundError, match="No module named 'scipy'"
+                    ):
+                        self.gdf.plot(kind=kind)
+            elif kind in _y_kinds:
+                kwargs = {"y": "y"}
+            elif kind in _xy_kinds:
+                kwargs = {"x": "x", "y": "y"}
+
+            self.compare_figures(kind, fig_test, fig_ref, kwargs)
+            plt.close("all")
+
+        @check_figures_equal(extensions=["png", "pdf"])
+        def test_geo_kind(self, fig_test, fig_ref):
+            """Test Geo kind."""
+            ax1 = fig_test.subplots()
+            self.gdf.plot(ax=ax1)
+            ax2 = fig_ref.subplots()
+            getattr(self.gdf.plot, "geo")(ax=ax2)
+            plt.close("all")
+
+    def test_invalid_kind(self):
+        """Test invalid kinds."""
+        with pytest.raises(ValueError, match="error is not a valid plot kind"):
+            self.gdf.plot(kind="error")
+        with pytest.raises(
+            AttributeError,
+            match="'GeoplotAccessor' object has no attribute 'error'",
+        ):
+            self.gdf.plot.error()
+
+
 def test_column_values():
     """
     Check that the dataframe plot method returns same values with an
diff --git a/geopandas/tests/test_testing.py b/geopandas/tests/test_testing.py
index 582c01b2..e3cec213 100644
--- a/geopandas/tests/test_testing.py
+++ b/geopandas/tests/test_testing.py
@@ -46,6 +46,14 @@ s4 = s1.copy()
 s4.crs = 4326
 s5 = s2.copy()
 s5.crs = 27700
+
+s6 = GeoSeries(
+    [
+        Polygon([(0, 3), (0, 0), (2, 0), (2, 2)]),
+        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
+    ]
+)
+
 df4 = GeoDataFrame(
     {"col1": [1, 2], "geometry": s1.copy(), "geom2": s4.copy(), "geom3": s5.copy()},
     crs=3857,
@@ -63,8 +71,15 @@ def test_geoseries():
     assert_geoseries_equal(s3, s2, check_series_type=False, check_dtype=False)
     assert_geoseries_equal(s1, s4, check_series_type=False)
 
-    with pytest.raises(AssertionError):
+    with pytest.raises(AssertionError) as error:
         assert_geoseries_equal(s1, s2, check_less_precise=True)
+    assert "1 out of 2 geometries are not almost equal" in str(error.value)
+    assert "not almost equal: [0]" in str(error.value)
+
+    with pytest.raises(AssertionError) as error:
+        assert_geoseries_equal(s2, s6, check_less_precise=False)
+    assert "1 out of 2 geometries are not equal" in str(error.value)
+    assert "not equal: [0]" in str(error.value)
 
 
 def test_geodataframe():