From 0dde1c7253ee4755f2ab8656c3bfa74c617e8d3a Mon Sep 17 00:00:00 2001 From: Leah Wasser Date: Fri, 7 Feb 2020 07:39:26 -0700 Subject: [PATCH] ENH: Add clip function (#1128) Co-authored-by: Martin Fleischmann Co-authored-by: Nathan Korinek Co-authored-by: Joris Van den Bossche --- doc/source/reference.rst | 3 +- examples/plot_clip.py | 116 ++++++++++ geopandas/__init__.py | 2 + geopandas/tools/__init__.py | 2 + geopandas/tools/clip.py | 252 ++++++++++++++++++++++ geopandas/tools/tests/test_clip.py | 327 +++++++++++++++++++++++++++++ 6 files changed, 701 insertions(+), 1 deletion(-) create mode 100644 examples/plot_clip.py create mode 100644 geopandas/tools/clip.py create mode 100644 geopandas/tools/tests/test_clip.py diff --git a/doc/source/reference.rst b/doc/source/reference.rst index 62f4403..a17b30a 100644 --- a/doc/source/reference.rst +++ b/doc/source/reference.rst @@ -157,8 +157,9 @@ API Pages GeoDataFrame GeoSeries - overlay read_file sjoin + overlay + clip tools.geocode datasets.get_path diff --git a/examples/plot_clip.py b/examples/plot_clip.py new file mode 100644 index 0000000..9a4b89d --- /dev/null +++ b/examples/plot_clip.py @@ -0,0 +1,116 @@ +""" +Clip Vector Data with GeoPandas +================================================================== + +Learn how to clip geometries to the boundary of a polygon geometry +using GeoPandas. + +.. currentmodule:: geopandas + +""" + +############################################################################### +# +# The example below shows you how to clip a set of vector geometries +# to the spatial extent / shape of another vector object. Both sets of geometries +# must be opened with GeoPandas as GeoDataFrames and be in the same Coordinate +# Reference System (CRS) for the :func:`clip` function in GeoPandas to work. +# +# This example uses GeoPandas example data ``'naturalearth_cities'`` and +# ``'naturalearth_lowres'``, alongside a custom rectangle geometry made with +# shapely and then turned into a GeoDataFrame. +# +# .. note:: +# The object to be clipped will be clipped to the full extent of the clip +# object. If there are multiple polygons in clip object, the input data will +# be clipped to the total boundary of all polygons in clip object. + +############################################################################### +# Import Packages +# --------------- +# +# To begin, import the needed packages. + +import matplotlib.pyplot as plt +import geopandas +from shapely.geometry import Polygon + +############################################################################### +# Get or Create Example Data +# -------------------------- +# +# Below, the example GeoPandas data is imported and opened as a GeoDataFrame. +# Additionally, a polygon is created with shapely and then converted into a +# GeoDataFrame with the same CRS as the GeoPandas world dataset. + +capitals = geopandas.read_file(geopandas.datasets.get_path("naturalearth_cities")) +world = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres")) + +# Create a subset of the world data that is just the South American continent +south_america = world[world["continent"] == "South America"] + +# Create a custom polygon +polygon = Polygon([(0, 0), (0, 90), (180, 90), (180, 0), (0, 0)]) +poly_gdf = geopandas.GeoDataFrame([1], geometry=[polygon], crs=world.crs) + +############################################################################### +# Plot the Unclipped Data +# ----------------------- + +fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) +world.plot(ax=ax1) +poly_gdf.boundary.plot(ax=ax1, color="red") +south_america.boundary.plot(ax=ax2, color="green") +capitals.plot(ax=ax2, color="purple") +ax1.set_title("All Unclipped World Data", fontsize=20) +ax2.set_title("All Unclipped Capital Data", fontsize=20) +ax1.set_axis_off() +ax2.set_axis_off() +plt.show() + +############################################################################### +# Clip the Data +# -------------- +# +# When you call :func:`clip`, the first object called is the object that will +# be clipped. The second object called is the clip extent. The returned output +# will be a new clipped GeoDataframe. All of the attributes for each returned +# geometry will be retained when you clip. +# +# .. note:: +# Recall that the data must be in the same CRS in order to use the +# :func:`clip` function. If the data are not in the same CRS, be sure to use +# the GeoPandas :meth:`~GeoDataFrame.to_crs` method to ensure both datasets +# are in the same CRS. + +############################################################################### +# Clip the World Data +# -------------------- + +world_clipped = geopandas.clip(world, polygon) + +# Plot the clipped data +# The plot below shows the results of the clip function applied to the world +# sphinx_gallery_thumbnail_number = 2 +fig, ax = plt.subplots(figsize=(12, 8)) +world_clipped.plot(ax=ax, color="purple") +world.boundary.plot(ax=ax) +poly_gdf.boundary.plot(ax=ax, color="red") +ax.set_title("World Clipped", fontsize=20) +ax.set_axis_off() +plt.show() + +############################################################################### +# Clip the Capitals Data +# ---------------------- + +capitals_clipped = geopandas.clip(capitals, south_america) + +# Plot the clipped data +# The plot below shows the results of the clip function applied to the capital cities +fig, ax = plt.subplots(figsize=(12, 8)) +capitals_clipped.plot(ax=ax, color="purple") +south_america.boundary.plot(ax=ax, color="green") +ax.set_title("Capitals Clipped", fontsize=20) +ax.set_axis_off() +plt.show() diff --git a/geopandas/__init__.py b/geopandas/__init__.py index 314a1f9..91c5e6d 100644 --- a/geopandas/__init__.py +++ b/geopandas/__init__.py @@ -7,6 +7,8 @@ from geopandas.io.sql import read_postgis # noqa from geopandas.tools import sjoin # noqa from geopandas.tools import overlay # noqa from geopandas.tools._show_versions import show_versions # noqa +from geopandas.tools import clip # noqa + import geopandas.datasets # noqa diff --git a/geopandas/tools/__init__.py b/geopandas/tools/__init__.py index 0b5b149..6c5fa61 100644 --- a/geopandas/tools/__init__.py +++ b/geopandas/tools/__init__.py @@ -3,6 +3,7 @@ from .geocoding import geocode, reverse_geocode from .overlay import overlay from .sjoin import sjoin from .util import collect +from .clip import clip __all__ = [ "collect", @@ -11,4 +12,5 @@ __all__ = [ "overlay", "reverse_geocode", "sjoin", + "clip", ] diff --git a/geopandas/tools/clip.py b/geopandas/tools/clip.py new file mode 100644 index 0000000..da022b8 --- /dev/null +++ b/geopandas/tools/clip.py @@ -0,0 +1,252 @@ +""" +geopandas.clip +============== + +A module to clip vector data using GeoPandas. + +""" +import warnings + +import numpy as np +import pandas as pd + +from shapely.geometry import Polygon, MultiPolygon + +from geopandas import GeoDataFrame, GeoSeries + + +def _clip_points(gdf, poly): + """Clip point geometry to the polygon extent. + + Clip an input point GeoDataFrame to the polygon extent of the poly + parameter. Points that intersect the poly geometry are extracted with + associated attributes and returned. + + Parameters + ---------- + gdf : GeoDataFrame, GeoSeries + Composed of point geometry that will be clipped to the poly. + + poly : (Multi)Polygon + Reference geometry used to spatially clip the data. + + Returns + ------- + GeoDataFrame + The returned GeoDataFrame is a subset of gdf that intersects + with poly. + """ + spatial_index = gdf.sindex + bbox = poly.bounds + sidx = list(spatial_index.intersection(bbox)) + gdf_sub = gdf.iloc[sidx] + + return gdf_sub[gdf_sub.geometry.intersects(poly)] + + +def _clip_line_poly(gdf, poly): + """Clip line and polygon geometry to the polygon extent. + + Clip an input line or polygon to the polygon extent of the poly + parameter. Parts of Lines or Polygons that intersect the poly geometry are + extracted with associated attributes and returned. + + Parameters + ---------- + gdf : GeoDataFrame, GeoSeries + Line or polygon geometry that is clipped to poly. + + poly : (Multi)Polygon + Reference polygon for clipping. + + Returns + ------- + GeoDataFrame + The returned GeoDataFrame is a clipped subset of gdf + that intersects with poly. + """ + spatial_index = gdf.sindex + + # Create a box for the initial intersection + bbox = poly.bounds + # Get a list of id's for each object that overlaps the bounding box and + # subset the data to just those lines + sidx = list(spatial_index.intersection(bbox)) + gdf_sub = gdf.iloc[sidx] + + # Clip the data with the polygon + if isinstance(gdf_sub, GeoDataFrame): + clipped = gdf_sub.copy() + clipped["geometry"] = gdf_sub.intersection(poly) + + # Return the clipped layer with no null geometry values or empty geometries + return clipped[~clipped.geometry.is_empty & clipped.geometry.notnull()] + else: + # GeoSeries + clipped = gdf_sub.intersection(poly) + return clipped[~clipped.is_empty & clipped.notnull()] + + +def clip(gdf, mask, keep_geom_type=False): + """Clip points, lines, or polygon geometries to the mask extent. + + Both layers must be in the same Coordinate Reference System (CRS). + The `gdf` will be clipped to the full extent of the clip object. + + If there are multiple polygons in mask, data from `gdf` will be + clipped to the total boundary of all polygons in mask. + + Parameters + ---------- + gdf : GeoDataFrame or GeoSeries + Vector layer (point, line, polygon) to be clipped to mask. + mask : GeoDataFrame, GeoSeries, (Multi)Polygon + Polygon vector layer used to clip `gdf`. + The mask's geometry is dissolved into one geometric feature + and intersected with `gdf`. + keep_geom_type : boolean, default False + If True, return only geometries of original type in case of intersection + resulting in multiple geometry types or GeometryCollections. + If False, return all resulting geometries (potentially mixed-types). + + Returns + ------- + GeoDataFrame or GeoSeries + Vector data (points, lines, polygons) from `gdf` clipped to + polygon boundary from mask. + + Examples + -------- + Clip points (global cities) with a polygon (the South American continent): + + >>> import geopandas + >>> path = + >>> world = geopandas.read_file( + ... geopandas.datasets.get_path('naturalearth_lowres')) + >>> south_america = world[world['continent'] == "South America"] + >>> capitals = geopandas.read_file( + ... geopandas.datasets.get_path('naturalearth_cities')) + >>> capitals.shape + (202, 2) + >>> sa_capitals = geopandas.clip(capitals, south_america) + >>> sa_capitals.shape + (12, 2) + """ + if not isinstance(gdf, (GeoDataFrame, GeoSeries)): + raise TypeError( + "'gdf' should be GeoDataFrame or GeoSeries, got {}".format(type(gdf)) + ) + + if not isinstance(mask, (GeoDataFrame, GeoSeries, Polygon, MultiPolygon)): + raise TypeError( + "'mask' should be GeoDataFrame, GeoSeries or" + "(Multi)Polygon, got {}".format(type(gdf)) + ) + + if isinstance(mask, (GeoDataFrame, GeoSeries)): + box_mask = mask.total_bounds + else: + box_mask = mask.bounds + box_gdf = gdf.total_bounds + if not ( + ((box_mask[0] <= box_gdf[2]) and (box_gdf[0] <= box_mask[2])) + and ((box_mask[1] <= box_gdf[3]) and (box_gdf[1] <= box_mask[3])) + ): + return GeoDataFrame(columns=gdf.columns, crs=gdf.crs) + + if isinstance(mask, (GeoDataFrame, GeoSeries)): + poly = mask.geometry.unary_union + else: + poly = mask + + geom_types = gdf.geometry.type + poly_idx = np.asarray((geom_types == "Polygon") | (geom_types == "MultiPolygon")) + line_idx = np.asarray( + (geom_types == "LineString") + | (geom_types == "LinearRing") + | (geom_types == "MultiLineString") + ) + point_idx = np.asarray((geom_types == "Point") | (geom_types == "MultiPoint")) + geomcoll_idx = np.asarray((geom_types == "GeometryCollection")) + + if point_idx.any(): + point_gdf = _clip_points(gdf[point_idx], poly) + else: + point_gdf = None + + if poly_idx.any(): + poly_gdf = _clip_line_poly(gdf[poly_idx], poly) + else: + poly_gdf = None + + if line_idx.any(): + line_gdf = _clip_line_poly(gdf[line_idx], poly) + else: + line_gdf = None + + if geomcoll_idx.any(): + geomcoll_gdf = _clip_line_poly(gdf[geomcoll_idx], poly) + else: + geomcoll_gdf = None + + order = pd.Series(range(len(gdf)), index=gdf.index) + concat = pd.concat([point_gdf, line_gdf, poly_gdf, geomcoll_gdf]) + + if keep_geom_type: + geomcoll_concat = (concat.geom_type == "GeometryCollection").any() + geomcoll_orig = geomcoll_idx.any() + + new_collection = geomcoll_concat and not geomcoll_orig + + if geomcoll_orig: + warnings.warn( + "keep_geom_type can not be called on a " + "GeoDataFrame with GeometryCollection." + ) + else: + polys = ["Polygon", "MultiPolygon"] + lines = ["LineString", "MultiLineString", "LinearRing"] + points = ["Point", "MultiPoint"] + + # Check that the gdf for multiple geom types (points, lines and/or polys) + orig_types_total = sum( + [ + gdf.geom_type.isin(polys).any(), + gdf.geom_type.isin(lines).any(), + gdf.geom_type.isin(points).any(), + ] + ) + + # Check how many geometry types are in the clipped GeoDataFrame + clip_types_total = sum( + [ + concat.geom_type.isin(polys).any(), + concat.geom_type.isin(lines).any(), + concat.geom_type.isin(points).any(), + ] + ) + + # Check there aren't any new geom types in the clipped GeoDataFrame + more_types = orig_types_total < clip_types_total + + if orig_types_total > 1: + warnings.warn( + "keep_geom_type can not be called on a mixed type GeoDataFrame." + ) + elif new_collection or more_types: + orig_type = gdf.geom_type.iloc[0] + if new_collection: + concat = concat.explode() + if orig_type in polys: + concat = concat.loc[concat.geom_type.isin(polys)] + elif orig_type in lines: + concat = concat.loc[concat.geom_type.isin(lines)] + + # preserve the original order of the input + if isinstance(concat, GeoDataFrame): + concat["_order"] = order + return concat.sort_values(by="_order").drop(columns="_order") + else: + concat = GeoDataFrame(geometry=concat) + concat["_order"] = order + return concat.sort_values(by="_order").geometry diff --git a/geopandas/tools/tests/test_clip.py b/geopandas/tools/tests/test_clip.py new file mode 100644 index 0000000..64aea51 --- /dev/null +++ b/geopandas/tools/tests/test_clip.py @@ -0,0 +1,327 @@ +"""Tests for the clip module.""" + +import warnings + +import numpy as np + +import shapely +from shapely.geometry import Polygon, Point, LineString, LinearRing, GeometryCollection + +import geopandas +from geopandas import GeoDataFrame, GeoSeries, clip +from geopandas.testing import assert_geodataframe_equal + +import pytest + + +@pytest.fixture +def point_gdf(): + """Create a point GeoDataFrame.""" + pts = np.array([[2, 2], [3, 4], [9, 8], [-12, -15]]) + gdf = GeoDataFrame([Point(xy) for xy in pts], columns=["geometry"], crs="EPSG:4326") + return gdf + + +@pytest.fixture +def single_rectangle_gdf(): + """Create a single rectangle for clipping.""" + poly_inters = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]) + gdf = GeoDataFrame([1], geometry=[poly_inters], crs="EPSG:4326") + gdf["attr2"] = "site-boundary" + return gdf + + +@pytest.fixture +def larger_single_rectangle_gdf(): + """Create a slightly larger rectangle for clipping. + The smaller single rectangle is used to test the edge case where slivers + are returned when you clip polygons. This fixture is larger which + eliminates the slivers in the clip return. + """ + poly_inters = Polygon([(-5, -5), (-5, 15), (15, 15), (15, -5), (-5, -5)]) + gdf = GeoDataFrame([1], geometry=[poly_inters], crs="EPSG:4326") + gdf["attr2"] = ["study area"] + return gdf + + +@pytest.fixture +def buffered_locations(point_gdf): + """Buffer points to create a multi-polygon.""" + buffered_locs = point_gdf + buffered_locs["geometry"] = buffered_locs.buffer(4) + buffered_locs["type"] = "plot" + return buffered_locs + + +@pytest.fixture +def donut_geometry(buffered_locations, single_rectangle_gdf): + """Make a geometry with a hole in the middle (a donut).""" + donut = geopandas.overlay( + buffered_locations, single_rectangle_gdf, how="symmetric_difference" + ) + return donut + + +@pytest.fixture +def two_line_gdf(): + """Create Line Objects For Testing""" + linea = LineString([(1, 1), (2, 2), (3, 2), (5, 3)]) + lineb = LineString([(3, 4), (5, 7), (12, 2), (10, 5), (9, 7.5)]) + gdf = GeoDataFrame([1, 2], geometry=[linea, lineb], crs="EPSG:4326") + return gdf + + +@pytest.fixture +def multi_poly_gdf(donut_geometry): + """Create a multi-polygon GeoDataFrame.""" + multi_poly = donut_geometry.unary_union + out_df = GeoDataFrame(geometry=GeoSeries(multi_poly), crs="EPSG:4326") + out_df["attr"] = ["pool"] + return out_df + + +@pytest.fixture +def multi_line(two_line_gdf): + """Create a multi-line GeoDataFrame. + This GDF has one multiline and one regular line.""" + # Create a single and multi line object + multiline_feat = two_line_gdf.unary_union + linec = LineString([(2, 1), (3, 1), (4, 1), (5, 2)]) + out_df = GeoDataFrame(geometry=GeoSeries([multiline_feat, linec]), crs="EPSG:4326") + out_df["attr"] = ["road", "stream"] + return out_df + + +@pytest.fixture +def multi_point(point_gdf): + """Create a multi-point GeoDataFrame.""" + multi_point = point_gdf.unary_union + out_df = GeoDataFrame( + geometry=GeoSeries( + [multi_point, Point(2, 5), Point(-11, -14), Point(-10, -12)] + ), + crs="EPSG:4326", + ) + out_df["attr"] = ["tree", "another tree", "shrub", "berries"] + return out_df + + +@pytest.fixture +def mixed_gdf(): + """Create a Mixed Polygon and LineString For Testing""" + point = Point([(2, 3), (11, 4), (7, 2), (8, 9), (1, 13)]) + line = LineString([(1, 1), (2, 2), (3, 2), (5, 3), (12, 1)]) + poly = Polygon([(3, 4), (5, 2), (12, 2), (10, 5), (9, 7.5)]) + ring = LinearRing([(1, 1), (2, 2), (3, 2), (5, 3), (12, 1)]) + gdf = GeoDataFrame( + [1, 2, 3, 4], geometry=[point, poly, line, ring], crs="EPSG:4326" + ) + return gdf + + +@pytest.fixture +def geomcol_gdf(): + """Create a Mixed Polygon and LineString For Testing""" + point = Point([(2, 3), (11, 4), (7, 2), (8, 9), (1, 13)]) + poly = Polygon([(3, 4), (5, 2), (12, 2), (10, 5), (9, 7.5)]) + coll = GeometryCollection([point, poly]) + gdf = GeoDataFrame([1], geometry=[coll], crs="EPSG:4326") + return gdf + + +@pytest.fixture +def sliver_line(): + """Create a line that will create a point when clipped.""" + linea = LineString([(10, 5), (13, 5), (15, 5)]) + lineb = LineString([(1, 1), (2, 2), (3, 2), (5, 3), (12, 1)]) + gdf = GeoDataFrame([1, 2], geometry=[linea, lineb], crs="EPSG:4326") + return gdf + + +def test_not_gdf(single_rectangle_gdf): + """Non-GeoDataFrame inputs raise attribute errors.""" + with pytest.raises(TypeError): + clip((2, 3), single_rectangle_gdf) + with pytest.raises(TypeError): + clip(single_rectangle_gdf, (2, 3)) + + +def test_returns_gdf(point_gdf, single_rectangle_gdf): + """Test that function returns a GeoDataFrame (or GDF-like) object.""" + out = clip(point_gdf, single_rectangle_gdf) + assert isinstance(out, GeoDataFrame) + + +def test_returns_series(point_gdf, single_rectangle_gdf): + """Test that function returns a GeoSeries if GeoSeries is passed.""" + out = clip(point_gdf.geometry, single_rectangle_gdf) + assert isinstance(out, GeoSeries) + + +def test_non_overlapping_geoms(): + """Test that a bounding box returns error if the extents don't overlap""" + unit_box = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]) + unit_gdf = GeoDataFrame([1], geometry=[unit_box], crs="EPSG:4326") + non_overlapping_gdf = unit_gdf.copy() + non_overlapping_gdf = non_overlapping_gdf.geometry.apply( + lambda x: shapely.affinity.translate(x, xoff=20) + ) + out = clip(unit_gdf, non_overlapping_gdf) + assert_geodataframe_equal( + out, GeoDataFrame(columns=unit_gdf.columns, crs=unit_gdf.crs) + ) + + +def test_clip_points(point_gdf, single_rectangle_gdf): + """Test clipping a points GDF with a generic polygon geometry.""" + clip_pts = clip(point_gdf, single_rectangle_gdf) + pts = np.array([[2, 2], [3, 4], [9, 8]]) + exp = GeoDataFrame([Point(xy) for xy in pts], columns=["geometry"], crs="EPSG:4326") + assert_geodataframe_equal(clip_pts, exp) + + +def test_clip_poly(buffered_locations, single_rectangle_gdf): + """Test clipping a polygon GDF with a generic polygon geometry.""" + clipped_poly = clip(buffered_locations, single_rectangle_gdf) + assert len(clipped_poly.geometry) == 3 + assert all(clipped_poly.geom_type == "Polygon") + + +def test_clip_poly_series(buffered_locations, single_rectangle_gdf): + """Test clipping a polygon GDF with a generic polygon geometry.""" + clipped_poly = clip(buffered_locations.geometry, single_rectangle_gdf) + assert len(clipped_poly) == 3 + assert all(clipped_poly.geom_type == "Polygon") + + +def test_clip_multipoly_keep_slivers(multi_poly_gdf, single_rectangle_gdf): + """Test a multi poly object where the return includes a sliver. + Also the bounds of the object should == the bounds of the clip object + if they fully overlap (as they do in these fixtures).""" + clipped = clip(multi_poly_gdf, single_rectangle_gdf) + assert np.array_equal(clipped.total_bounds, single_rectangle_gdf.total_bounds) + # Assert returned data is a geometry collection given sliver geoms + assert "GeometryCollection" in clipped.geom_type[0] + + +def test_clip_multipoly_keep_geom_type(multi_poly_gdf, single_rectangle_gdf): + """Test a multi poly object where the return includes a sliver. + Also the bounds of the object should == the bounds of the clip object + if they fully overlap (as they do in these fixtures).""" + clipped = clip(multi_poly_gdf, single_rectangle_gdf, keep_geom_type=True) + assert np.array_equal(clipped.total_bounds, single_rectangle_gdf.total_bounds) + # Assert returned data is a not geometry collection + assert (clipped.geom_type == "Polygon").any() + + +def test_clip_single_multipoly_no_extra_geoms( + buffered_locations, larger_single_rectangle_gdf +): + """When clipping a multi-polygon feature, no additional geom types + should be returned.""" + multi = buffered_locations.dissolve(by="type").reset_index() + clipped = clip(multi, larger_single_rectangle_gdf) + assert clipped.geom_type[0] == "Polygon" + + +def test_clip_multiline(multi_line, single_rectangle_gdf): + """Test that clipping a multiline feature with a poly returns expected output.""" + clipped = clip(multi_line, single_rectangle_gdf) + assert clipped.geom_type[0] == "MultiLineString" + + +def test_clip_multipoint(single_rectangle_gdf, multi_point): + """Clipping a multipoint feature with a polygon works as expected. + should return a geodataframe with a single multi point feature""" + clipped = clip(multi_point, single_rectangle_gdf) + assert clipped.geom_type[0] == "MultiPoint" + assert hasattr(clipped, "attr") + # All points should intersect the clip geom + assert all(clipped.intersects(single_rectangle_gdf.unary_union)) + + +def test_clip_lines(two_line_gdf, single_rectangle_gdf): + """Test what happens when you give the clip_extent a line GDF.""" + clip_line = clip(two_line_gdf, single_rectangle_gdf) + assert len(clip_line.geometry) == 2 + + +def test_clip_with_multipolygon(buffered_locations, single_rectangle_gdf): + """Test clipping a polygon with a multipolygon.""" + multi = buffered_locations.dissolve(by="type").reset_index() + clipped = clip(single_rectangle_gdf, multi) + assert clipped.geom_type[0] == "Polygon" + + +def test_mixed_geom(mixed_gdf, single_rectangle_gdf): + """Test clipping a mixed GeoDataFrame""" + clipped = clip(mixed_gdf, single_rectangle_gdf) + assert ( + clipped.geom_type[0] == "Point" + and clipped.geom_type[1] == "Polygon" + and clipped.geom_type[2] == "LineString" + ) + + +def test_mixed_series(mixed_gdf, single_rectangle_gdf): + """Test clipping a mixed GeoSeries""" + clipped = clip(mixed_gdf.geometry, single_rectangle_gdf) + assert ( + clipped.geom_type[0] == "Point" + and clipped.geom_type[1] == "Polygon" + and clipped.geom_type[2] == "LineString" + ) + + +def test_clip_warning_no_extra_geoms(buffered_locations, single_rectangle_gdf): + """Test a user warning is provided if no new geometry types are found.""" + with pytest.warns(UserWarning): + clip(buffered_locations, single_rectangle_gdf, True) + warnings.warn( + "keep_geom_type was called when no extra geometry types existed.", + UserWarning, + ) + + +def test_clip_with_polygon(single_rectangle_gdf): + """Test clip when using a shapely object""" + polygon = Polygon([(0, 0), (5, 12), (10, 0), (0, 0)]) + clipped = clip(single_rectangle_gdf, polygon) + exp_poly = polygon.intersection( + Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]) + ) + exp = GeoDataFrame([1], geometry=[exp_poly], crs="EPSG:4326") + exp["attr2"] = "site-boundary" + assert_geodataframe_equal(clipped, exp) + + +def test_clip_with_line_extra_geom(single_rectangle_gdf, sliver_line): + """When the output of a clipped line returns a geom collection, + and keep_geom_type is True, no geometry collections should be returned.""" + clipped = clip(sliver_line, single_rectangle_gdf, keep_geom_type=True) + assert len(clipped.geometry) == 1 + # Assert returned data is a not geometry collection + assert not (clipped.geom_type == "GeometryCollection").any() + + +def test_clip_line_keep_slivers(single_rectangle_gdf, sliver_line): + """Test the correct output if a point is returned + from a line only geometry type.""" + clipped = clip(sliver_line, single_rectangle_gdf) + # Assert returned data is a geometry collection given sliver geoms + assert "Point" == clipped.geom_type[0] + assert "LineString" == clipped.geom_type[1] + + +def test_warning_extra_geoms_mixed(single_rectangle_gdf, mixed_gdf): + """Test the correct warnings are raised if keep_geom_type is + called on a mixed GDF""" + with pytest.warns(UserWarning): + clip(mixed_gdf, single_rectangle_gdf, keep_geom_type=True) + + +def test_warning_geomcoll(single_rectangle_gdf, geomcol_gdf): + """Test the correct warnings are raised if keep_geom_type is + called on a GDF with GeometryCollection""" + with pytest.warns(UserWarning): + clip(geomcol_gdf, single_rectangle_gdf, keep_geom_type=True)