wassname/geopandas
1
0
Fork 0
mirror of https://github.com/wassname/geopandas.git synced 2026-09-24 13:30:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

ENH: Add clip function (#1128)

Browse Source 
Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>
Co-authored-by: Nathan Korinek <nako1890@colorado.edu>
Co-authored-by: Joris Van den Bossche <jorisvandenbossche@gmail.com>
This commit is contained in:
Leah Wasser
2020-02-07 15:39:26 +01:00
committed by GitHub
co-authored by Martin Fleischmann Nathan Korinek Joris Van den Bossche
parent 81272cbafd
commit 0dde1c7253
6 changed files with 701 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/source/reference.rst
+2 -1
View File
@@ -157,8 +157,9 @@ API Pages
GeoDataFrame
GeoSeries
overlay
read_file
sjoin
overlay
clip
tools.geocode
datasets.get_path
examples/plot_clip.py
+116
View File
@@ -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()
geopandas/__init__.py
+2
View File
@@ -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
geopandas/tools/__init__.py
+2
View File
@@ -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",
]
geopandas/tools/clip.py
+252
View File
@@ -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
geopandas/tools/tests/test_clip.py
+327
View File
@@ -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)