Repository: geopandas/geopandas
Pull request #2224 (2021-11-09 21:53:12)

BUG: keep colors aligned with shapes, even when empty
_plot_polygon_collection eliminates empty shapes, but doesn't eliminate corresponding colors.  Operations like intersection that perform cropping can accidentally change colors of remaining shapes.

Here's a simple demonstration of the problem:

```python
import geopandas as gpd
from shapely.geometry import box

s = gpd.GeoSeries([
    box(0,0,1,1),
    box(2,2,3,3),
], index=["r", "b"])

# correct:
s.plot(color=["red", "blue"])

# incorrect: should show blue, but shows red
s2 = s.intersection(box(2,2,4,4))
s2.plot(color=["red", "blue"])
```

martinfleis: Thanks for the catch and a fix! A few notes:

1. Can we filter the input array of colours to remove those that match empty geometries rather than creating and empty Patch? It seems like a bit dirty solution to me as it is now.
1. This will likely affect other geom types as well. Can you check if it is the case for LineStrings and Points and fix those in the same way if needed?
2. Can you add some tests to `test_plotting.py` to ensure we won't break it in the future?
tylerharter: Thanks for the feedback!

1. I was originally trying to do something like that, but it was turning out to be clunky/buggy, probably because I'm not familiar with the code.  Happy to revisit (in the next couple weeks) if you have some pointers for an elegant way.  Do you think _flatten_multi_geoms would be the place to attempt that?  Should I be using is_list_like to determine which attributes to filter (to match filtering of empty geoms)?
2. good call.  Points was fine.  LineStrings had a different bug, which I also fixed in the last commit (can't do a [:,:2] slice on an empty array).
3. sounds good, just added them
AkiyamaKunka: Awesome catch!
martinfleis: I have refactored this PR to skip empty or missing geometries instead of creating empty patches. It is done in the `_flatten_multi_geoms` which has been renamed to `_sanitize_geoms` for this reason.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8dd90834..53fb4d2f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,9 +2,20 @@
 
 ## Development version
 
+New features and improvements:
+
+- Added ``minimum_bounding_circle()`` method from shapely to GeoSeries/GeoDataframe (#2621).
+
+Deprecations and compatibility notes:
+
+- Added warning that ``unary_union`` will return ``'GEOMETRYCOLLECTION EMPTY'`` instead
+  of None for all-None GeoSeries. (#2618)
+
 Bug fixes:
 
 - Ensure that GeoDataFrame created from DataFrame is a copy, not a view (#2667)
+- Fix mismatch between geometries and colors in ``plot()`` if an empty or missing
+  geometry is present (#2224)
 
 ## Version 0.12.2 (December 10, 2022)
 
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index bf7f4ec2..7d7317b9 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -86,6 +86,7 @@ Constructive methods and attributes
    GeoSeries.centroid
    GeoSeries.convex_hull
    GeoSeries.envelope
+   GeoSeries.minimum_bounding_circle
    GeoSeries.simplify
    GeoSeries.normalize
 
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index d1f132b4..266da762 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -20,6 +20,7 @@ PANDAS_GE_115 = Version(pd.__version__) >= Version("1.1.5")
 PANDAS_GE_12 = Version(pd.__version__) >= Version("1.2.0")
 PANDAS_GE_13 = Version(pd.__version__) >= Version("1.3.0")
 PANDAS_GE_14 = Version(pd.__version__) >= Version("1.4.0rc0")
+PANDAS_GE_15 = Version(pd.__version__) >= Version("1.5.0")
 PANDAS_GE_20 = Version(pd.__version__) >= Version("2.0.0.dev0")
 
 
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 20aa7332..36408804 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -677,6 +677,18 @@ def representative_point(data):
         return out
 
 
+def minimum_bounding_circle(data):
+    if compat.USE_SHAPELY_20:
+        return shapely.minimum_bounding_circle(data)
+    elif compat.USE_PYGEOS:
+        return pygeos.minimum_bounding_circle(data)
+    else:
+        raise NotImplementedError(
+            f"shapely >= 2.0 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+
+
 #
 # Binary predicates
 #
@@ -1021,21 +1033,41 @@ def relate(data, other):
 
 
 def unary_union(data):
+    warning_msg = (
+        "`unary_union` returned None due to all-None GeoSeries. In future, "
+        "`unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead."
+    )
+
     if compat.USE_SHAPELY_20:
         data = shapely.union_all(data)
-        if data is None:  # shapely 2.0a1
-            return None
-        if data.is_empty:  # shapely 2.0
+        if data is None or data.is_empty:  # shapely 2.0a1 and 2.0
+            warnings.warn(
+                warning_msg,
+                FutureWarning,
+                stacklevel=4,
+            )
             return None
         else:
             return data
     elif compat.USE_PYGEOS:
-        return _pygeos_to_shapely(pygeos.union_all(data))
+        result = _pygeos_to_shapely(pygeos.union_all(data))
+        if result is None:
+            warnings.warn(
+                warning_msg,
+                FutureWarning,
+                stacklevel=4,
+            )
+        return result
     else:
         data = [g for g in data if g is not None]
         if data:
             return shapely.ops.unary_union(data)
         else:
+            warnings.warn(
+                warning_msg,
+                FutureWarning,
+                stacklevel=4,
+            )
             return None
 
 
diff --git a/geopandas/array.py b/geopandas/array.py
index 5f472333..22fceb36 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -511,6 +511,11 @@ class GeometryArray(ExtensionArray):
     def representative_point(self):
         return GeometryArray(vectorized.representative_point(self.data), crs=self.crs)
 
+    def minimum_bounding_circle(self):
+        return GeometryArray(
+            vectorized.minimum_bounding_circle(self.data), crs=self.crs
+        )
+
     def normalize(self):
         return GeometryArray(vectorized.normalize(self.data), crs=self.crs)
 
@@ -1329,9 +1334,8 @@ class GeometryArray(ExtensionArray):
             # TODO(pygeos)
             return getattr(to_shapely(self), name)()
         raise TypeError(
-            "cannot perform {name} with type {dtype}".format(
-                name=name, dtype=self.dtype
-            )
+            f"'{type(self).__name__}' with dtype {self.dtype} "
+            f"does not support reduction '{name}'"
         )
 
     def __array__(self, dtype=None):
diff --git a/geopandas/base.py b/geopandas/base.py
index ada4c79d..b4c6a908 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -692,6 +692,39 @@ GeometryCollection
         """
         return _delegate_geo_method("representative_point", self)
 
+    def minimum_bounding_circle(self):
+        """Returns a ``GeoSeries`` of geometries representing the minimum bounding
+        circle that encloses each geometry.
+
+        Examples
+        --------
+
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1), (0, 1), (0, 0)]),
+        ...         LineString([(0, 0), (1, 1), (1, 0)]),
+        ...         Point(0, 0),
+        ...     ]
+        ... )
+        >>> s
+        0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
+        1    LINESTRING (0.00000 0.00000, 1.00000 1.00000, ...
+        2                              POINT (0.00000 0.00000)
+        dtype: geometry
+
+        >>> s.minimum_bounding_circle()
+        0    POLYGON ((1.20711 0.50000, 1.19352 0.36205, 1....
+        1    POLYGON ((1.20711 0.50000, 1.19352 0.36205, 1....
+        2                              POINT (0.00000 0.00000)
+        dtype: geometry
+
+        See also
+        --------
+        GeoSeries.convex_hull : convex hull geometry
+        """
+        return _delegate_geo_method("minimum_bounding_circle", self)
+
     def normalize(self):
         """Returns a ``GeoSeries`` of normalized
         geometries to normal form (or canonical form).
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 7be42d28..6b5eb009 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1581,6 +1581,7 @@ individually so that features may have different properties
         sort=True,
         observed=False,
         dropna=True,
+        **kwargs,
     ):
         """
         Dissolve geometries within `groupby` into single observation.
@@ -1633,7 +1634,12 @@ individually so that features may have different properties
             if a non-default value is given for this parameter.
 
             .. versionadded:: 0.9.0
+        **kwargs :
+            Keyword arguments to be passed to the pandas `DataFrameGroupby.agg` method
+            which is used by `dissolve`. In particular, `numeric_only` may be
+            supplied, which will be required in pandas 2.0 for certain aggfuncs.
 
+            .. versionadded:: 0.13.0
         Returns
         -------
         GeoDataFrame
@@ -1679,7 +1685,26 @@ individually so that features may have different properties
 
         # Process non-spatial component
         data = self.drop(labels=self.geometry.name, axis=1)
-        aggregated_data = data.groupby(**groupby_kwargs).agg(aggfunc)
+        with warnings.catch_warnings(record=True) as record:
+            aggregated_data = data.groupby(**groupby_kwargs).agg(aggfunc, **kwargs)
+        for w in record:
+            if str(w.message).startswith("The default value of numeric_only"):
+                msg = (
+                    f"The default value of numeric_only in aggfunc='{aggfunc}' "
+                    "within pandas.DataFrameGroupBy.agg used in dissolve is "
+                    "deprecated. In pandas 2.0, numeric_only will default to False. "
+                    "Either specify numeric_only as additional argument in dissolve() "
+                    "or select only columns which should be valid for the function."
+                )
+                warnings.warn(msg, FutureWarning, stacklevel=2)
+            else:
+                # Only want to capture specific warning,
+                # other warnings from pandas should be passed through
+                # TODO this is not an ideal approach
+                warnings.showwarning(
+                    w.message, w.category, w.filename, w.lineno, w.file, w.line
+                )
+
         aggregated_data.columns = aggregated_data.columns.to_flat_index()
 
         # Process spatial component
diff --git a/geopandas/plotting.py b/geopandas/plotting.py
index 4145307f..91a9f49d 100644
--- a/geopandas/plotting.py
+++ b/geopandas/plotting.py
@@ -26,11 +26,12 @@ def deprecated(new, warning_type=FutureWarning):
     return old
 
 
-def _flatten_multi_geoms(geoms, prefix="Multi"):
+def _sanitize_geoms(geoms, prefix="Multi"):
     """
     Returns Series like geoms and index, except that any Multi geometries
     are split into their components and indices are repeated for all component
-    in the same Multi geometry.  Maintains 1:1 matching of geometry to value.
+    in the same Multi geometry. At the same time, empty or missing geometries are
+    filtered out.  Maintains 1:1 matching of geometry to value.
 
     Prefix specifies type of geometry to be flatten. 'Multi' for MultiPoint and similar,
     "Geom" for GeometryCollection.
@@ -42,9 +43,15 @@ def _flatten_multi_geoms(geoms, prefix="Multi"):
     component_index : index array
         indices are repeated for all components in the same Multi geometry
     """
+    # TODO(shapely) look into simplifying this with
+    # shapely.get_parts(geoms, return_index=True) from shapely 2.0
     components, component_index = [], []
 
-    if not geoms.geom_type.str.startswith(prefix).any():
+    if (
+        not geoms.geom_type.str.startswith(prefix).any()
+        and not geoms.is_empty.any()
+        and not geoms.isna().any()
+    ):
         return geoms, np.arange(len(geoms))
 
     for ix, geom in enumerate(geoms):
@@ -52,6 +59,8 @@ def _flatten_multi_geoms(geoms, prefix="Multi"):
             for poly in geom.geoms:
                 components.append(poly)
                 component_index.append(ix)
+        elif geom is None or geom.is_empty:
+            continue
         else:
             components.append(geom)
             component_index.append(ix)
@@ -154,7 +163,7 @@ def _plot_polygon_collection(
     """
     from matplotlib.collections import PatchCollection
 
-    geoms, multiindex = _flatten_multi_geoms(geoms)
+    geoms, multiindex = _sanitize_geoms(geoms)
     if values is not None:
         values = np.take(values, multiindex, axis=0)
 
@@ -171,9 +180,7 @@ def _plot_polygon_collection(
 
     _expand_kwargs(kwargs, multiindex)
 
-    collection = PatchCollection(
-        [_PolygonPatch(poly) for poly in geoms if not poly.is_empty], **kwargs
-    )
+    collection = PatchCollection([_PolygonPatch(poly) for poly in geoms], **kwargs)
 
     if values is not None:
         collection.set_array(np.asarray(values))
@@ -213,7 +220,7 @@ def _plot_linestring_collection(
     """
     from matplotlib.collections import LineCollection
 
-    geoms, multiindex = _flatten_multi_geoms(geoms)
+    geoms, multiindex = _sanitize_geoms(geoms)
     if values is not None:
         values = np.take(values, multiindex, axis=0)
 
@@ -283,7 +290,7 @@ def _plot_point_collection(
     if values is not None and color is not None:
         raise ValueError("Can only specify one of 'values' and 'color' kwargs")
 
-    geoms, multiindex = _flatten_multi_geoms(geoms)
+    geoms, multiindex = _sanitize_geoms(geoms)
     # values are expanded below as kwargs["c"]
 
     x = [p.x if not p.is_empty else None for p in geoms]
@@ -427,7 +434,7 @@ def plot_series(
         style_kwds["vmax"] = style_kwds.get("vmax", values.max())
 
     # decompose GeometryCollections
-    geoms, multiindex = _flatten_multi_geoms(s.geometry, prefix="Geom")
+    geoms, multiindex = _sanitize_geoms(s.geometry, prefix="Geom")
     values = np.take(values, multiindex, axis=0) if cmap else None
     # ensure indexes are consistent
     if color_given and isinstance(color, pd.Series):
@@ -829,7 +836,7 @@ GON (((-122.84000 49.00000, -120.0000...
     mx = values[~np.isnan(values)].max() if vmax is None else vmax
 
     # decompose GeometryCollections
-    geoms, multiindex = _flatten_multi_geoms(df.geometry, prefix="Geom")
+    geoms, multiindex = _sanitize_geoms(df.geometry, prefix="Geom")
     values = np.take(values, multiindex, axis=0)
     nan_idx = np.take(nan_idx, multiindex, axis=0)
     expl_series = geopandas.GeoSeries(geoms)


Test patch:
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index dd0eafc9..6d9b740e 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -10,6 +10,7 @@ from geopandas import _compat as compat
 from pandas.testing import assert_frame_equal
 import pytest
 
+from geopandas._compat import PANDAS_GE_15, PANDAS_GE_20
 from geopandas.testing import assert_geodataframe_equal
 
 
@@ -82,11 +83,52 @@ def test_first_dissolve(nybb_polydf, first):
 
 
 def test_mean_dissolve(nybb_polydf, first, expected_mean):
-    test = nybb_polydf.dissolve("manhattan_bronx", aggfunc="mean")
-    assert_frame_equal(expected_mean, test, check_column_type=False)
+    if not PANDAS_GE_15:
+        test = nybb_polydf.dissolve("manhattan_bronx", aggfunc="mean")
+        test2 = nybb_polydf.dissolve("manhattan_bronx", aggfunc=np.mean)
+    elif PANDAS_GE_15 and not PANDAS_GE_20:
+        with pytest.warns(FutureWarning, match=".*used in dissolve is deprecated.*"):
+            test = nybb_polydf.dissolve("manhattan_bronx", aggfunc="mean")
+            test2 = nybb_polydf.dissolve("manhattan_bronx", aggfunc=np.mean)
+    else:  # pandas 2.0
+        test = nybb_polydf.dissolve(
+            "manhattan_bronx", aggfunc="mean", numeric_only=True
+        )
+        # for non pandas "mean", numeric only cannot be applied. Drop columns manually
+        test2 = nybb_polydf.drop(columns=["BoroName"]).dissolve(
+            "manhattan_bronx", aggfunc=np.mean
+        )
 
-    test = nybb_polydf.dissolve("manhattan_bronx", aggfunc=np.mean)
     assert_frame_equal(expected_mean, test, check_column_type=False)
+    assert_frame_equal(expected_mean, test2, check_column_type=False)
+
+
+@pytest.mark.skipif(not PANDAS_GE_15 or PANDAS_GE_20, reason="warning for pandas 1.5.x")
+def test_mean_dissolve_warning_capture(nybb_polydf, first, expected_mean):
+    with pytest.warns(
+        FutureWarning,
+        match=".*used in dissolve is deprecated.*",
+    ):
+        nybb_polydf.dissolve("manhattan_bronx", aggfunc="mean")
+
+    # test no warning for aggfunc first which doesn't have numeric only semantics
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        nybb_polydf.dissolve("manhattan_bronx", aggfunc="first")
+
+
+def test_dissolve_emits_other_warnings(nybb_polydf):
+    # we only do something special for pandas 1.5.x, but expect this
+    # test to be true on any version
+    def sum_and_warn(group):
+        warnings.warn("foo")
+        if PANDAS_GE_20:
+            return group.sum(numeric_only=False)
+        else:
+            return group.sum()
+
+    with pytest.warns(UserWarning, match="foo"):
+        nybb_polydf.dissolve("manhattan_bronx", aggfunc=sum_and_warn)
 
 
 def test_multicolumn_dissolve(nybb_polydf, first):
@@ -123,7 +165,7 @@ def test_dissolve_none(nybb_polydf):
 
 
 def test_dissolve_none_mean(nybb_polydf):
-    test = nybb_polydf.dissolve(aggfunc="mean")
+    test = nybb_polydf.dissolve(aggfunc="mean", numeric_only=True)
     expected = GeoDataFrame(
         {
             nybb_polydf.geometry.name: [nybb_polydf.geometry.unary_union],
diff --git a/geopandas/tests/test_explore.py b/geopandas/tests/test_explore.py
index 3ccd49ed..7e4a8544 100644
--- a/geopandas/tests/test_explore.py
+++ b/geopandas/tests/test_explore.py
@@ -14,6 +14,7 @@ import matplotlib.colors as colors  # noqa
 from branca.colormap import StepColormap  # noqa
 
 BRANCA_05 = Version(branca.__version__) > Version("0.4.2")
+FOLIUM_G_014 = Version(folium.__version__) > Version("0.14.0")
 
 
 class TestExplore:
@@ -524,8 +525,14 @@ class TestExplore:
         df2["values"] = df2["BoroCode"] * 10.0
         m = df2[df2["values"] >= 30].explore("values", vmin=0)
         out_str = self._fetch_map_string(m)
-        assert 'case"1":return{"color":"#7ad151","fillColor":"#7ad151"' in out_str
-        assert 'case"2":return{"color":"#22a884","fillColor":"#22a884"' in out_str
+        if FOLIUM_G_014:
+            assert 'case"0":return{"color":"#fde725","fillColor":"#fde725"' in out_str
+            assert 'case"1":return{"color":"#7ad151","fillColor":"#7ad151"' in out_str
+            assert 'default:return{"color":"#22a884","fillColor":"#22a884"' in out_str
+        else:
+            assert 'case"1":return{"color":"#7ad151","fillColor":"#7ad151"' in out_str
+            assert 'case"2":return{"color":"#22a884","fillColor":"#22a884"' in out_str
+            assert 'default:return{"color":"#fde725","fillColor":"#fde725"' in out_str
 
         df2["values_negative"] = df2["BoroCode"] * -10.0
         m = df2[df2["values_negative"] <= 30].explore("values_negative", vmax=0)
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 606102e3..eccd5c00 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -385,8 +385,9 @@ class TestGeomMethods:
         g2 = GeoSeries([p1, None])
         self._test_unary_topological("unary_union", p1, g2)
 
-        g3 = GeoSeries([None, None])
-        assert g3.unary_union is None
+        with pytest.warns(FutureWarning, match="`unary_union` returned None"):
+            g3 = GeoSeries([None, None])
+            assert g3.unary_union is None
 
     def test_cascaded_union_deprecated(self):
         p1 = self.t1
@@ -927,6 +928,21 @@ class TestGeomMethods:
         assert isinstance(e, GeoSeries)
         assert self.g3.crs == e.crs
 
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="minimum_bounding_circle is only implemented for pygeos, not shapely",
+    )
+    def test_minimum_bounding_circle(self):
+        mbc = self.g1.minimum_bounding_circle()
+        centers = GeoSeries([Point(0.5, 0.5)] * 2)
+        assert np.all(mbc.centroid.geom_equals_exact(centers, 0.001))
+        assert_series_equal(
+            mbc.area,
+            Series([1.560723, 1.560723]),
+        )
+        assert isinstance(mbc, GeoSeries)
+        assert self.g1.crs == mbc.crs
+
     def test_total_bounds(self):
         bbox = self.sol.x, self.sol.y, self.esb.x, self.esb.y
         assert isinstance(self.landmarks.total_bounds, np.ndarray)
diff --git a/geopandas/tests/test_op_output_types.py b/geopandas/tests/test_op_output_types.py
index 5a447f39..789982b9 100644
--- a/geopandas/tests/test_op_output_types.py
+++ b/geopandas/tests/test_op_output_types.py
@@ -348,5 +348,5 @@ def test_constructor_sliced_in_pandas_methods(df2):
     # drop the secondary geometry columns as not hashable
     hashable_test_df = df2.drop(columns=["geometry2", "geometry3"])
     assert type(hashable_test_df.duplicated()) == pd.Series
-    assert type(df2.quantile()) == pd.Series
+    assert type(df2.quantile(numeric_only=True)) == pd.Series
     assert type(df2.memory_usage()) == pd.Series
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index 69f34071..e92f04c7 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -9,6 +9,7 @@ from shapely.geometry import Point, Polygon, LineString, GeometryCollection, box
 import geopandas
 from geopandas import GeoDataFrame, GeoSeries, overlay, read_file
 from geopandas import _compat
+from geopandas._compat import PANDAS_GE_20
 
 from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
 import pytest
@@ -755,13 +756,18 @@ def test_non_overlapping(how):
     result = overlay(df1, df2, how=how)
 
     if how == "intersection":
+        if PANDAS_GE_20:
+            index = None
+        else:
+            index = pd.Index([], dtype="object")
+
         expected = GeoDataFrame(
             {
                 "col1": np.array([], dtype="int64"),
                 "col2": np.array([], dtype="int64"),
                 "geometry": [],
             },
-            index=pd.Index([], dtype="object"),
+            index=index,
         )
     elif how == "union":
         expected = GeoDataFrame(
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 27d45f1d..3875adaf 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -326,7 +326,11 @@ def test_to_csv(df):
 def test_numerical_operations(s, df):
     # df methods ignore the geometry column
     exp = pd.Series([3, 4], index=["value1", "value2"])
-    assert_series_equal(df.sum(), exp)
+    if not compat.PANDAS_GE_20:
+        res = df.sum()
+    else:
+        res = df.sum(numeric_only=True)
+    assert_series_equal(res, exp)
 
     # series methods raise error (not supported for geometry)
     with pytest.raises(TypeError):
@@ -343,8 +347,7 @@ def test_numerical_operations(s, df):
     with pytest.raises(TypeError):
         df + 1
 
-    with pytest.raises((TypeError, AssertionError)):
-        # TODO(pandas 0.23) remove AssertionError -> raised in 0.23
+    with pytest.raises(TypeError):
         s + 1
 
     # boolean comparisons work
@@ -536,7 +539,10 @@ def test_groupby(df):
     assert_frame_equal(res, exp)
 
     # reductions ignore geometry column
-    res = df.groupby("value2").sum()
+    if not compat.PANDAS_GE_20:
+        res = df.groupby("value2").sum()
+    else:
+        res = df.groupby("value2").sum(numeric_only=True)
     exp = pd.DataFrame({"value1": [2, 1], "value2": [1, 2]}, dtype="int64").set_index(
         "value2"
     )
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index 771eac82..2a4041eb 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -335,6 +335,36 @@ class TestPointPlotting:
         ax = gdf.plot()
         assert len(ax.collections) == 1
 
+    @pytest.mark.parametrize(
+        "geoms",
+        [
+            [
+                box(0, 0, 1, 1),
+                box(7, 7, 8, 8),
+            ],
+            [
+                LineString([(1, 1), (1, 2)]),
+                LineString([(7, 1), (7, 2)]),
+            ],
+            [
+                Point(1, 1),
+                Point(7, 7),
+            ],
+        ],
+    )
+    def test_empty_geometry_colors(self, geoms):
+        s = GeoSeries(
+            geoms,
+            index=["r", "b"],
+        )
+        s2 = s.intersection(box(5, 0, 10, 10))
+        ax = s2.plot(color=["red", "blue"])
+        blue = np.array([0.0, 0.0, 1.0, 1.0])
+        if s.geom_type["r"] == "LineString":
+            np.testing.assert_array_equal(ax.get_children()[0].get_edgecolor()[0], blue)
+        else:
+            np.testing.assert_array_equal(ax.get_children()[0].get_facecolor()[0], blue)
+
     def test_multipoints(self):
 
         # MultiPoints