"""
earthpy.clip
============
A module to clip vector data using GeoPandas.
"""
import pandas as pd
import geopandas as gpd
def _clip_points(shp, clip_obj):
"""Clip point geometry to the clip_obj GeoDataFrame extent.
Clip an input point GeoDataFrame to the polygon extent of the clip_obj
parameter. Points that intersect the clip_obj geometry are extracted with
associated attributes and returned.
Parameters
----------
shp : GeoDataFrame
Composed of point geometry that is clipped to clip_obj.
clip_obj : GeoDataFrame
Reference polygon for clipping.
Returns
-------
GeoDataFrame
The returned GeoDataFrame is a subset of shp that intersects
with clip_obj.
"""
poly = clip_obj.geometry.unary_union
return shp[shp.geometry.intersects(poly)]
def _clip_multi_point(shp, clip_obj):
"""Clip multi point features to the clip_obj GeoDataFrame extent.
Clip an input multi point to the polygon extent of the clip_obj
parameter. Points that intersect the clip_obj geometry are
extracted with associated attributes returned.
Parameters
----------
shp : GeoDataFrame
multipoint geometry that is clipped to clip_obj.
clip_obj : GeoDataFrame
Reference polygon for clipping.
Returns
-------
GeoDataFrame
The returned GeoDataFrame is a clipped subset of shp
containing multi-point and point features.
"""
# Explode multi-point features when clipping then recreate geom
clipped = _clip_points(shp.explode().reset_index(level=[1]), clip_obj)
clipped = clipped.dissolve(by=[clipped.index]).drop(columns="level_1")[
shp.columns.tolist()
]
return clipped
def _clip_line_poly(shp, clip_obj):
"""Clip line and polygon geometry to the clip_obj GeoDataFrame extent.
Clip an input line or polygon to the polygon extent of the clip_obj
parameter. Lines or Polygons that intersect the clip_obj geometry are
extracted with associated attributes and returned.
Parameters
----------
shp : GeoDataFrame
Line or polygon geometry that is clipped to clip_obj.
clip_obj : GeoDataFrame
Reference polygon for clipping.
Returns
-------
GeoDataFrame
The returned GeoDataFrame is a clipped subset of shp
that intersects with clip_obj.
"""
# Create a single polygon object for clipping
poly = clip_obj.geometry.unary_union
spatial_index = shp.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))
shp_sub = shp.iloc[sidx]
# Clip the data - with these data
clipped = shp_sub.copy()
clipped["geometry"] = shp_sub.intersection(poly)
# Return the clipped layer with no null geometry values
return clipped[clipped.geometry.notnull()]
def _clip_multi_poly_line(shp, clip_obj):
"""Clip multi lines and polygons to the clip_obj GeoDataFrame extent.
Clip an input multi line or polygon to the polygon extent of the clip_obj
parameter. Lines or Polygons that intersect the clip_obj geometry are
extracted with associated attributes and returned.
Parameters
----------
shp : GeoDataFrame
multiLine or multipolygon geometry that is clipped to clip_obj.
clip_obj : GeoDataFrame
Reference polygon for clipping.
Returns
-------
GeoDataFrame
The returned GeoDataFrame is a clipped subset of shp
that intersects with clip_obj.
"""
# Clip multi polygons
clipped = _clip_line_poly(shp.explode().reset_index(level=[1]), clip_obj)
lines = clipped[
(clipped.geometry.type == "MultiLineString")
| (clipped.geometry.type == "LineString")
]
line_diss = lines.dissolve(by=[lines.index]).drop(columns="level_1")
polys = clipped[clipped.geometry.type == "Polygon"]
poly_diss = polys.dissolve(by=[polys.index]).drop(columns="level_1")
return gpd.GeoDataFrame(
pd.concat([poly_diss, line_diss], ignore_index=True)
)
[docs]def clip_shp(shp, clip_obj):
"""Clip points, lines, or polygon geometries to the clip_obj extent.
Both layers must be in the same Coordinate Reference System (CRS) and will
be clipped to the full extent of the clip object.
If there are multiple polygons in clip_obj,
data from shp will be clipped to the total boundary of
all polygons in clip_obj.
Parameters
----------
shp : GeoDataFrame
Vector layer (point, line, polygon) to be clipped to clip_obj.
clip_obj : GeoDataFrame
Polygon vector layer used to clip shp.
The clip_obj's geometry is dissolved into one geometric feature
and intersected with shp.
Returns
-------
GeoDataFrame
Vector data (points, lines, polygons) from shp clipped to
polygon boundary from clip_obj.
Examples
--------
Clipping points (glacier locations in the state of Colorado) with
a polygon (the boundary of Rocky Mountain National Park):
>>> import geopandas as gpd
>>> import earthpy.clip as cl
>>> from earthpy.io import path_to_example
>>> rmnp = gpd.read_file(path_to_example('rmnp.shp'))
>>> glaciers = gpd.read_file(path_to_example('colorado-glaciers.geojson'))
>>> glaciers.shape
(134, 2)
>>> rmnp_glaciers = cl.clip_shp(glaciers, rmnp)
>>> rmnp_glaciers.shape
(36, 2)
Clipping a line (the Continental Divide Trail) with a
polygon (the boundary of Rocky Mountain National Park):
>>> cdt = gpd.read_file(path_to_example('continental-div-trail.geojson'))
>>> rmnp_cdt_section = cl.clip_shp(cdt, rmnp)
>>> cdt['geometry'].length > rmnp_cdt_section['geometry'].length
0 True
dtype: bool
Clipping a polygon (Colorado counties) with another polygon
(the boundary of Rocky Mountain National Park):
>>> counties = gpd.read_file(path_to_example('colorado-counties.geojson'))
>>> counties.shape
(64, 13)
>>> rmnp_counties = cl.clip_shp(counties, rmnp)
>>> rmnp_counties.shape
(4, 13)
"""
try:
shp.geometry
clip_obj.geometry
except AttributeError:
raise AttributeError(
"Please make sure that your input and clip GeoDataFrames have a"
" valid geometry column"
)
if not any(shp.intersects(clip_obj.unary_union)):
raise ValueError("Shape and crop extent do not overlap.")
if any(shp.geometry.type == "MultiPoint"):
return _clip_multi_point(shp, clip_obj)
elif any(shp.geometry.type == "Point"):
return _clip_points(shp, clip_obj)
elif any(shp.geometry.type == "MultiPolygon") or any(
shp.geometry.type == "MultiLineString"
):
return _clip_multi_poly_line(shp, clip_obj)
else:
return _clip_line_poly(shp, clip_obj)