Repository: geopandas/geopandas

Issue:
BUG: Clipping with empty shape produces cryptic error
- [x] I have checked that this issue has not already been reported.

- [x] I have confirmed this bug exists on the latest version of geopandas.

---

**Note**: Please read [this guide](https://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports) detailing how to provide the necessary information for us to reproduce your bug.

#### Code Sample, a copy-pastable example

```python
# Your code here
import geopandas as gpd
import numpy as np
from shapely.geometry import Point

gs = gpd.GeoSeries.from_xy(np.arange(4), np.arange(4)).buffer(1)
mask = Point(np.nan, np.nan).buffer(1)
gpd.clip(gs, mask)
```

#### Problem description

Passing an empty mask shape into `clip` produces a cryptic `IndexError`.

The above code results in the following error:

```python
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/USER/anaconda3/envs/rstools/lib/python3.8/site-packages/geopandas/tools/clip.py", line 174, in clip
    ((box_mask[0] <= box_gdf[2]) and (box_gdf[0] <= box_mask[2]))
IndexError: tuple index out of range
```

An error should be raised here, but the error that is raised is confusing and doesn't help the user understand what has happened. I think bad input data has been allowed to get too far into the machinery. 

#### Expected Output
I think that `clip` should raise a `ValueError` (or whatever error is most appropriate) indicating that the provided mask was empty. This is easy to check for with the `.is_empty` property.

#### Output of ``geopandas.show_versions()``

<details>

SYSTEM INFO
-----------
python     : 3.8.13 | packaged by conda-forge | (default, Mar 25 2022, 06:04:10)  [GCC 10.3.0]
executable : /home/USER/anaconda3/envs/rstools/bin/python
machine    : Linux-5.15.0-47-generic-x86_64-with-glibc2.10

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.11.0
GEOS lib   : /home/USER/anaconda3/envs/rstools/lib/libgeos_c.so
GDAL       : 3.5.2
GDAL data dir: /home/USER/anaconda3/envs/rstools/share/gdal
PROJ       : 9.0.1
PROJ data dir: /home/USER/anaconda3/envs/rstools/share/proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.11.1
numpy      : 1.23.3
pandas     : 1.5.0
pyproj     : 3.4.0
shapely    : 1.8.4
fiona      : 1.8.21
geoalchemy2: None
geopy      : None
matplotlib : 3.6.0
mapclassify: 2.4.3
pygeos     : 0.13
pyogrio    : v0.4.2
psycopg2   : None
pyarrow    : None
rtree      : 1.0.0

</details>


Discussion:
fbunt: I went ahead and added tests to check that empty geometries, series, and dataframes all produce an empty result.
fbunt: Well now I'm confused. It looks like the tests are failing because of NaN inputs to shapely geometries but as far as i can tell shapely itself is happy to accept them. Does anyone know what the critical difference is between the test environments?
jorisvandenbossche: Thanks a lot for the PR!

> Well now I'm confused. It looks like the tests are failing because of NaN inputs to shapely geometries but as far as i can tell shapely itself is happy to accept them. Does anyone know what the critical difference is between the test environments?

That seems to be an issue with older GEOS versions (older GEOS is failing to buffer a point if there are NaN coordinates). 

But I think we can avoid that by creating the empty geometry in a different way. For example, `Polygon()` should also create an empty geometry.

fbunt: Thanks for the help! I think it's good to go now. The failed test seems to be mamba related.
jorisvandenbossche: Thanks for the update!

> The failed test seems to be mamba related.

Yes, there was some flaky CI, I merged main, that should help now

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b1af58cc..8977f083 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@ for more details.
 New features and improvements:
 
 - Added ``normalize()`` method from shapely to GeoSeries/GeoDataframe (#2537)
+- Added ``make_valid()`` method from shapely to GeoSeries/GeoDataframe (#2539)
 - Add ``where`` filter to ``read_file`` (#2552)
 
 Deprecations and compatibility notes:
@@ -26,9 +27,11 @@ Deprecations and compatibility notes:
 - resolve ``matplotlib.cm`` warning in ``.explore()`` (#2596)
 
 Bug fixes:
+- Fix cryptic error message in ``geopandas.clip()`` when clipping with an empty geometry (#2589)
 - Accessing `gdf.geometry` where the active geometry column is missing, and a column named `"geometry"` is present 
   will now raise an `AttributeError`, rather than returning `gdf["geometry"]` (#2575) 
-
+- Combining GeoSeries/GeoDataFrames with ``pandas.concat`` will no longer silently
+  override CRS information if not all inputs have the same CRS (#2056).
 
 Version 0.11.1 (July 24, 2022)
 ------------------------------
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index e5736f2a..45d383dd 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -15,6 +15,7 @@ import shapely.geos
 import shapely.ops
 import shapely.wkb
 import shapely.wkt
+import shapely.validation
 
 from shapely.geometry.base import BaseGeometry
 
@@ -975,6 +976,26 @@ def normalize(data):
     return out
 
 
+def make_valid(data):
+    if compat.USE_SHAPELY_20:
+        return shapely.make_valid(data)
+    elif compat.USE_PYGEOS:
+        return pygeos.make_valid(data)
+    elif not compat.SHAPELY_GE_18:
+        raise NotImplementedError(
+            f"shapely >= 1.8 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+    else:
+        out = np.empty(len(data), dtype=object)
+        with compat.ignore_shapely2_warnings():
+            out[:] = [
+                shapely.validation.make_valid(geom) if geom is not None else None
+                for geom in data
+            ]
+    return out
+
+
 def project(data, other, normalized=False):
     if compat.USE_SHAPELY_20:
         return shapely.line_locate_point(data, other, normalized=normalized)
diff --git a/geopandas/array.py b/geopandas/array.py
index c2065565..779e7291 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -514,6 +514,9 @@ class GeometryArray(ExtensionArray):
     def normalize(self):
         return GeometryArray(vectorized.normalize(self.data), crs=self.crs)
 
+    def make_valid(self):
+        return GeometryArray(vectorized.make_valid(self.data), crs=self.crs)
+
     #
     # Binary predicates
     #
@@ -1317,7 +1320,7 @@ class GeometryArray(ExtensionArray):
         ExtensionArray
         """
         data = np.concatenate([ga.data for ga in to_concat])
-        return GeometryArray(data, crs=to_concat[0].crs)
+        return GeometryArray(data, crs=_get_common_crs(to_concat))
 
     def _reduce(self, name, skipna=True, **kwargs):
         # including the base class version here (that raises by default)
@@ -1388,3 +1391,26 @@ class GeometryArray(ExtensionArray):
             else:
                 return False
         return (self == item).any()
+
+
+def _get_common_crs(arr_seq):
+
+    crs_set = {arr.crs for arr in arr_seq}
+    crs_not_none = [crs for crs in crs_set if crs is not None]
+    names = [crs.name for crs in crs_not_none]
+
+    if len(crs_not_none) == 0:
+        return None
+    if len(crs_not_none) == 1:
+        if len(crs_set) != 1:
+            warnings.warn(
+                "CRS not set for some of the concatenation inputs. "
+                f"Setting output's CRS as {names[0]} "
+                "(the single non-null crs provided)."
+            )
+        return crs_not_none[0]
+
+    raise ValueError(
+        f"Cannot determine common CRS for concatenation inputs, got {names}. "
+        "Use `to_crs()` to transform geometries to the same CRS before merging."
+    )
diff --git a/geopandas/base.py b/geopandas/base.py
index 437a6119..8302b324 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -725,6 +725,45 @@ GeometryCollection
         """
         return _delegate_geo_method("normalize", self)
 
+    def make_valid(self):
+        """
+        Repairs invalid geometries.
+
+        Returns a ``GeoSeries`` with valid geometries.
+        If the input geometry is already valid, then it will be preserved.
+        In many cases, in order to create a valid geometry, the input
+        geometry must be split into multiple parts or multiple geometries.
+        If the geometry must be split into multiple parts of the same type
+        to be made valid, then a multi-part geometry will be returned
+        (e.g. a MultiPolygon).
+        If the geometry must be split into multiple parts of different types
+        to be made valid, then a GeometryCollection will be returned.
+
+        Examples
+        --------
+        >>> from shapely.geometry import MultiPolygon, Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)]),
+        ...         Polygon([(0, 2), (0, 1), (2, 0), (0, 0), (0, 2)]),
+        ...         LineString([(0, 0), (1, 1), (1, 0)]),
+        ...     ],
+        ...     crs='EPSG:3857',
+        ... )
+        >>> s
+        0    POLYGON ((0.000 0.000, 0.000 2.000, 1.000 1.00...
+        1    POLYGON ((0.000 2.000, 0.000 1.000, 2.000 0.00...
+        2    LINESTRING (0.000 0.000, 1.000 1.000, 1.000 0....
+        dtype: geometry
+
+        >>> s.make_valid()
+        0    MULTIPOLYGON (((1.000 1.000, 0.000 0.000, 0.00...
+        1    GEOMETRYCOLLECTION (POLYGON ((2.000 0.000, 0.0...
+        2    LINESTRING (0.000 0.000, 1.000 1.000, 1.000 0....
+        dtype: geometry
+        """
+        return _delegate_geo_method("make_valid", self)
+
     #
     # Reduction operations that return a Shapely geometry
     #
diff --git a/geopandas/tools/clip.py b/geopandas/tools/clip.py
index c10b12ba..a088db89 100644
--- a/geopandas/tools/clip.py
+++ b/geopandas/tools/clip.py
@@ -7,6 +7,7 @@ A module to clip vector data using GeoPandas.
 """
 import warnings
 
+import numpy as np
 import pandas.api.types
 from shapely.geometry import Polygon, MultiPolygon, box
 
@@ -168,7 +169,11 @@ def clip(gdf, mask, keep_geom_type=False):
     elif mask_is_list_like:
         box_mask = mask
     else:
-        box_mask = mask.bounds
+        # Avoid empty tuple returned by .bounds when geometry is empty. A tuple of
+        # all nan values is consistent with the behavior of
+        # {GeoSeries, GeoDataFrame}.total_bounds for empty geometries.
+        # TODO(shapely) can simpely use mask.bounds once relying on Shapely 2.0
+        box_mask = mask.bounds if not mask.is_empty else (np.nan,) * 4
     box_gdf = gdf.total_bounds
     if not (
         ((box_mask[0] <= box_gdf[2]) and (box_gdf[0] <= box_mask[2]))


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index d29932dc..606102e3 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -5,7 +5,16 @@ import numpy as np
 from numpy.testing import assert_array_equal
 from pandas import DataFrame, Index, MultiIndex, Series
 
-from shapely.geometry import LinearRing, LineString, MultiPoint, Point, Polygon
+import shapely
+
+from shapely.geometry import (
+    LinearRing,
+    LineString,
+    MultiPoint,
+    Point,
+    Polygon,
+    MultiPolygon,
+)
 from shapely.geometry.collection import GeometryCollection
 from shapely.ops import unary_union
 from shapely import wkt
@@ -682,6 +691,43 @@ class TestGeomMethods:
         expected = GeoSeries([polygon2, linestring, point])
         assert_geoseries_equal(series.normalize(), expected)
 
+    @pytest.mark.skipif(
+        not compat.SHAPELY_GE_18,
+        reason="make_valid keyword introduced in shapely 1.8.0",
+    )
+    def test_make_valid(self):
+        polygon1 = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
+        polygon2 = Polygon([(0, 2), (0, 1), (2, 0), (0, 0), (0, 2)])
+        linestring = LineString([(0, 0), (1, 1), (1, 0)])
+        series = GeoSeries([polygon1, polygon2, linestring])
+        out_polygon1 = MultiPolygon(
+            [
+                Polygon([(1, 1), (0, 0), (0, 2), (1, 1)]),
+                Polygon([(2, 0), (1, 1), (2, 2), (2, 0)]),
+            ]
+        )
+        out_polygon2 = GeometryCollection(
+            [Polygon([(2, 0), (0, 0), (0, 1), (2, 0)]), LineString([(0, 2), (0, 1)])]
+        )
+        expected = GeoSeries([out_polygon1, out_polygon2, linestring])
+        assert not series.is_valid.all()
+        result = series.make_valid()
+        assert_geoseries_equal(result, expected)
+        assert result.is_valid.all()
+
+    @pytest.mark.skipif(
+        compat.SHAPELY_GE_18,
+        reason="make_valid keyword introduced in shapely 1.8.0",
+    )
+    def test_make_valid_shapely_pre18(self):
+        s = GeoSeries([Point(1, 1)])
+        with pytest.raises(
+            NotImplementedError,
+            match=f"shapely >= 1.8 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed",
+        ):
+            s.make_valid()
+
     def test_convex_hull(self):
         # the convex hull of a square should be the same as the square
         squares = GeoSeries([self.sq for i in range(3)])
diff --git a/geopandas/tests/test_merge.py b/geopandas/tests/test_merge.py
index 5794b6ca..3f65346a 100644
--- a/geopandas/tests/test_merge.py
+++ b/geopandas/tests/test_merge.py
@@ -54,11 +54,6 @@ class TestMerging:
         self._check_metadata(res)
         exp = GeoDataFrame(pd.concat([pd.DataFrame(self.gdf), pd.DataFrame(self.gdf)]))
         assert_geodataframe_equal(exp, res)
-        # check metadata comes from first gdf
-        res4 = pd.concat([self.gdf.set_crs("epsg:4326"), self.gdf], axis=0)
-        # Note: this behaviour potentially does not make sense. If geom cols are
-        # concatenated but have different CRS, then the CRS will be overridden.
-        self._check_metadata(res4, crs="epsg:4326")
 
         # series
         res = pd.concat([self.gdf.geometry, self.gdf.geometry])
@@ -66,6 +61,48 @@ class TestMerging:
         assert isinstance(res, GeoSeries)
         assert isinstance(res.geometry, GeoSeries)
 
+    def test_concat_axis0_crs(self):
+
+        # CRS not set for both GeoDataFrame
+        res = pd.concat([self.gdf, self.gdf])
+        self._check_metadata(res)
+
+        # CRS set for both GeoDataFrame, same CRS
+        res1 = pd.concat([self.gdf.set_crs("epsg:4326"), self.gdf.set_crs("epsg:4326")])
+        self._check_metadata(res1, crs="epsg:4326")
+
+        # CRS not set for one GeoDataFrame, but set for the other GeoDataFrame
+        with pytest.warns(
+            UserWarning, match=r"CRS not set for some of the concatenation inputs.*"
+        ):
+            res2 = pd.concat([self.gdf, self.gdf.set_crs("epsg:4326")])
+            self._check_metadata(res2, crs="epsg:4326")
+
+        # CRS set for both GeoDataFrame, different CRS
+        with pytest.raises(
+            ValueError, match=r"Cannot determine common CRS for concatenation inputs.*"
+        ):
+            pd.concat([self.gdf.set_crs("epsg:4326"), self.gdf.set_crs("epsg:4327")])
+
+        # CRS not set for one GeoDataFrame, but set for the other GeoDataFrames,
+        # same CRS
+        with pytest.warns(
+            UserWarning, match=r"CRS not set for some of the concatenation inputs.*"
+        ):
+            res3 = pd.concat(
+                [self.gdf, self.gdf.set_crs("epsg:4326"), self.gdf.set_crs("epsg:4326")]
+            )
+            self._check_metadata(res3, crs="epsg:4326")
+
+        # CRS not set for one GeoDataFrame, but set for the other GeoDataFrames,
+        # different CRS
+        with pytest.raises(
+            ValueError, match=r"Cannot determine common CRS for concatenation inputs.*"
+        ):
+            pd.concat(
+                [self.gdf, self.gdf.set_crs("epsg:4326"), self.gdf.set_crs("epsg:4327")]
+            )
+
     def test_concat_axis1(self):
 
         res = pd.concat([self.gdf, self.df], axis=1)
diff --git a/geopandas/tools/tests/test_clip.py b/geopandas/tools/tests/test_clip.py
index c643d5ac..29006fda 100644
--- a/geopandas/tools/tests/test_clip.py
+++ b/geopandas/tools/tests/test_clip.py
@@ -450,3 +450,28 @@ def test_clip_single_multipoly_no_extra_geoms(
     multi = buffered_locations.dissolve(by="type").reset_index()
     clipped = clip(multi, masks)
     assert clipped.geom_type[0] == "Polygon"
+
+
+@pytest.mark.filterwarnings("ignore:All-NaN slice encountered")
+@pytest.mark.parametrize(
+    "mask",
+    [
+        Polygon(),
+        (np.nan,) * 4,
+        (np.nan, 0, np.nan, 1),
+        GeoSeries([Polygon(), Polygon()], crs="EPSG:3857"),
+        GeoSeries([Polygon(), Polygon()], crs="EPSG:3857").to_frame(),
+        GeoSeries([], crs="EPSG:3857"),
+        GeoSeries([], crs="EPSG:3857").to_frame(),
+    ],
+)
+def test_clip_empty_mask(buffered_locations, mask):
+    """Test that clipping with empty mask returns an empty result."""
+    clipped = clip(buffered_locations, mask)
+    assert_geodataframe_equal(
+        clipped,
+        GeoDataFrame([], columns=["geometry", "type"], crs="EPSG:3857"),
+        check_index_type=False,
+    )
+    clipped = clip(buffered_locations.geometry, mask)
+    assert_geoseries_equal(clipped, GeoSeries([], crs="EPSG:3857"))