备注
此页面是从 `gallery/plot_clip.ipynb`_ _.
使用GeoPandas裁剪矢量数据#
了解如何使用GeoPandas将几何图形剪裁到多边形几何图形的边界。
下面的示例显示如何将一组矢量几何图形剪裁到另一个矢量对象的空间范围/形状。这两组几何必须以GeoPandas作为GeoDataFrame打开,并且位于相同的坐标参考系(CRS)中 clip
功能在GeoPandas中工作。
此示例使用GeoPandas示例数据 'naturalearth_cities'
和 'naturalearth_lowres'
,旁边是使用Shapely制作的自定义矩形几何体,然后转换为GeoDataFrame。
注意事项
要剪裁的对象将被剪裁到剪裁对象的全部范围。如果裁剪对象中有多个多边形,则输入数据将被裁剪到裁剪对象中所有多边形的总边界。
导入包#
首先,导入所需的包。
[1]:
import matplotlib.pyplot as plt
import geopandas
from shapely.geometry import Polygon
获取或创建示例数据#
下面,示例GeoPandas数据被导入并作为GeoDataFrame打开。此外,使用Shapely创建一个多边形,然后将其转换为具有与GeoPandas世界数据集相同的CRS的GeoDataFrame。
[2]:
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)
绘制未剪裁的数据#
[3]:
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
是将被剪裁的对象。您传递的对象是剪裁范围。返回的输出将是新裁剪的GeoDataframe。剪裁时,每个返回的几何体的所有属性都将保留。
注意事项
回想一下,数据必须位于相同的CRS中才能使用 clip
方法。如果数据不在同一CRS中,请确保使用GeoPandas GeoDataFrame.to_crs
方法以确保两个数据集位于相同的CRS中。
剪裁世界数据#
[4]:
world_clipped = world.clip(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()
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 world_clipped = world.clip(polygon)
3 # Plot the clipped data
4 # The plot below shows the results of the clip function applied to the world
5 # sphinx_gallery_thumbnail_number = 2
6 fig, ax = plt.subplots(figsize=(12, 8))
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/geodataframe.py:2155, in GeoDataFrame.clip(self, mask, keep_geom_type)
2105 def clip(self, mask, keep_geom_type=False):
2106 """Clip points, lines, or polygon geometries to the mask extent.
2107
2108 Both layers must be in the same Coordinate Reference System (CRS).
(...)
2153 (12, 2)
2154 """
-> 2155 return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type)
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/tools/clip.py:184, in clip(gdf, mask, keep_geom_type)
181 else:
182 combined_mask = mask
--> 184 clipped = _clip_gdf_with_mask(gdf, combined_mask)
186 if keep_geom_type:
187 geomcoll_concat = (clipped.geom_type == "GeometryCollection").any()
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/tools/clip.py:49, in _clip_gdf_with_mask(gdf, mask)
46 else:
47 intersection_polygon = mask
---> 49 gdf_sub = gdf.iloc[gdf.sindex.query(intersection_polygon, predicate="intersects")]
51 # For performance reasons points don't need to be intersected with poly
52 non_point_mask = gdf_sub.geom_type != "Point"
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/base.py:2706, in GeoPandasBase.sindex(self)
2655 @property
2656 def sindex(self):
2657 """Generate the spatial index
2658
2659 Creates R-tree spatial index based on ``pygeos.STRtree`` or
(...)
2704 [2]])
2705 """
-> 2706 return self.geometry.values.sindex
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/array.py:291, in GeometryArray.sindex(self)
288 @property
289 def sindex(self):
290 if self._sindex is None:
--> 291 self._sindex = _get_sindex_class()(self.data)
292 return self._sindex
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/sindex.py:21, in _get_sindex_class()
19 if compat.HAS_RTREE:
20 return RTreeIndex
---> 21 raise ImportError(
22 "Spatial indexes require either `rtree` or `pygeos`. "
23 "See installation instructions at https://geopandas.org/install.html"
24 )
ImportError: Spatial indexes require either `rtree` or `pygeos`. See installation instructions at https://geopandas.org/install.html
注意事项
由于历史原因,裁剪方法也可以作为顶级函数使用 geopandas.clip
。建议使用该方法,因为该函数在将来可能会被弃用。
剪裁Capitals数据#
[5]:
capitals_clipped = capitals.clip(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()
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 capitals_clipped = capitals.clip(south_america)
3 # Plot the clipped data
4 # The plot below shows the results of the clip function applied to the capital cities
5 fig, ax = plt.subplots(figsize=(12, 8))
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/geodataframe.py:2155, in GeoDataFrame.clip(self, mask, keep_geom_type)
2105 def clip(self, mask, keep_geom_type=False):
2106 """Clip points, lines, or polygon geometries to the mask extent.
2107
2108 Both layers must be in the same Coordinate Reference System (CRS).
(...)
2153 (12, 2)
2154 """
-> 2155 return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type)
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/tools/clip.py:184, in clip(gdf, mask, keep_geom_type)
181 else:
182 combined_mask = mask
--> 184 clipped = _clip_gdf_with_mask(gdf, combined_mask)
186 if keep_geom_type:
187 geomcoll_concat = (clipped.geom_type == "GeometryCollection").any()
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/tools/clip.py:49, in _clip_gdf_with_mask(gdf, mask)
46 else:
47 intersection_polygon = mask
---> 49 gdf_sub = gdf.iloc[gdf.sindex.query(intersection_polygon, predicate="intersects")]
51 # For performance reasons points don't need to be intersected with poly
52 non_point_mask = gdf_sub.geom_type != "Point"
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/base.py:2706, in GeoPandasBase.sindex(self)
2655 @property
2656 def sindex(self):
2657 """Generate the spatial index
2658
2659 Creates R-tree spatial index based on ``pygeos.STRtree`` or
(...)
2704 [2]])
2705 """
-> 2706 return self.geometry.values.sindex
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/array.py:291, in GeometryArray.sindex(self)
288 @property
289 def sindex(self):
290 if self._sindex is None:
--> 291 self._sindex = _get_sindex_class()(self.data)
292 return self._sindex
File /usr/local/lib/python3.10/dist-packages/geopandas-0.10.2+79.g3abc6a7-py3.10.egg/geopandas/sindex.py:21, in _get_sindex_class()
19 if compat.HAS_RTREE:
20 return RTreeIndex
---> 21 raise ImportError(
22 "Spatial indexes require either `rtree` or `pygeos`. "
23 "See installation instructions at https://geopandas.org/install.html"
24 )
ImportError: Spatial indexes require either `rtree` or `pygeos`. See installation instructions at https://geopandas.org/install.html