Source code for earthpy.clip

"""
earthpy.clip
============

A module to clip vector data using GeoPandas.

"""

# TODO: Clip poly should use OVERLAY not spatial indexing + intersects


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_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()]


[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.") # Multipolys / point / line don't clip properly if "Multi" in str(clip_obj.geom_type) or "Multi" in str(shp.geom_type): raise ValueError( """Clip doesn't currently support multipart geometries. Consider using .explode to create unique features in your GeoDataFrame""" ) if shp["geometry"].iloc[0].type == "Point": return _clip_points(shp, clip_obj) else: return _clip_line_poly(shp, clip_obj)