3

I would like to use Shapely to find the vertices that are located within the boundaries of a Shapely Polygon. So far, I could only manage to get it running by converting the LineString and the Polygon to a GeoPandas DataFrame and use the Geopandas within(...) function. Ideally, the function that I am looking for would return a list of Shapely Points that were either kept or are the ones that were within the Polygon.

Here a little code example of what I have done so far:

fault_ls and interface_ls are Shapely Linestrings.

# Creating a buffer around the fault trace
fault_polygon = fault_ls.buffer(radius)

# Creating GeoDataFrame from Polygon
fault_polygon_gdf = gpd.GeoDataFrame({'geometry': [fault_polygon]}, crs=crs)

# Create lists with X and Y coordinates from LineString
x = [i[0] for i in interface_ls.coords]
y = [i[1] for i in interface_ls.coords]

# Creating GeoDataFrame from LineString
interface_ls_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x, y, crs=crs))

vertices_in = interface_ls_gdf[interface_ls_gdf.within(fault_polygon_gdf.loc[0, 'geometry'])]

vertices_out = interface_ls_gdf[~interface_ls_gdf.within(fault_polygon_gdf.loc[0, 'geometry'])]

vertices_in are the vertices within the buffer, vertices_out the ones outside that I want to keep.

1 Answer 1

2

Why not fully taking advantage of GeoPandas in that case?

Let do the imports first:

import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon, MultiPoint
import geopandas as gpd

Then read some dummy Shapefiles, one containing a single 'fault' line and the other containing a network of let's say pipelines close to the fault:

fault_file = 'fault.shp'
pipeline_file = 'pipeline.shp'
# Load with fiona
fault_gdf = gpd.read_file(fault_file)
pipeline_gdf = gpd.read_file(pipeline_file)

Check what's in there:

>>> fault_gdf
   id       name                                           geometry
0   1  the_fault  LINESTRING (-50.34740 64.08437, -50.34060 64.0...

>>> fault_gdf.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> pipeline_gdf
   id       name                                           geometry
0   1    primary  LINESTRING (-50.34277 64.09243, -50.33218 64.0...
1   2  secondary  LINESTRING (-50.27735 64.09024, -50.27442 64.0...
2   3  secondary  LINESTRING (-50.30269 64.08963, -50.30723 64.0...
3   4   tertiary  LINESTRING (-50.27442 64.08850, -50.27802 64.0...
4   5   tertiary  LINESTRING (-50.26885 64.08594, -50.26402 64.0...
5   6   tertiary  LINESTRING (-50.27064 64.08735, -50.27357 64.0...

Then build your buffer (warning, WGS84, buffer size in angle, not precise, you'd better work in a cartesian system):

fault_gdf['buffer'] = fault_gdf['geometry'].buffer(0.0032)
poly = fault_gdf.iloc[0]['buffer'] # Get the shape out for convenience 

Then build a function to extract vertices of these linestrings:

def get_vertice_from_line(gdf):
    gdf['points'] = gdf['geometry'].apply(lambda x: \
                                          MultiPoint([y for y in x.coords]))

    return gdf

And one function to get inside and outside vertices (which get written in two new columns as MultiPoints):

def get_points_in_and_out(gdf_lines, polygon):
    gdf_lines['in'] = gdf_lines['points'].intersection(poly)
    gdf_lines['out'] = gdf_lines['points'].difference(
        gdf_lines['points'].intersection(poly)
    )
    # Just in case; replace empty value by None to fix a bug when plotting
    gdf_lines['in'] = gdf_lines['in'].apply(lambda x: None if x.is_empty else x)
    gdf_lines['out'] = gdf_lines['out'].apply(lambda x: None if x.is_empty else x)
    gdf_lines.set_geometry('points', inplace=True)

    return gdf_lines

Call these two guys in a row:

pipeline_gdf = get_vertice_from_line(pipeline_gdf)
pipeline_gdf = get_points_in_and_out(pipeline_gdf, poly)

Check:

>>> pipeline_gdf
   id  ...                                                out
0   1  ...  MULTIPOINT (-50.34277 64.09243, -50.33218 64.0...
1   2  ...  MULTIPOINT (-50.27735 64.09024, -50.27442 64.0...
2   3  ...  MULTIPOINT (-50.34126 64.08875, -50.33757 64.0...
3   4  ...                         POINT (-50.27442 64.08850)
4   5  ...  MULTIPOINT (-50.25750 64.08607, -50.25240 64.0...
5   6  ...                         POINT (-50.27064 64.08735)

[6 rows x 6 columns]

Hmm, maybe better with an image:

Check

Plot for convenience:

def plot_gdf(fault_gdf, pipeline_gdf):
    fig, ax = plt.subplots(figsize = (20,16))
    fault_gdf.plot(color='maroon',ax=ax)
    fault_gdf['buffer'].plot(color='sandybrown', alpha=0.4, ax=ax, zorder=-9)
    pipeline_gdf['in'].plot(color='limegreen', markersize=500, ax=ax)
    pipeline_gdf['out'].plot(color='lightcoral', markersize=500, ax=ax)
    pipeline_gdf.plot(color='blue', marker='+', markersize=200, ax=ax)
    plt.grid()

plot_gdf(fault_gdf, pipeline_gdf)

And there we go:

Results

Compared with the original shapes in QGIS:

QGIS shapes

You can then export your two sets of points to a database for example, keeping whatever original attributes they have. Don't worry, shapely is still used under the hood, but you add the advantages of dealing with dataframes.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.