Clipping Shapefile in Pure Python

Question

I am trying to implement the code found at this site: http://geospatialpython.com/2015/05/clipping-shapefile-in-pure-python.html

The function I am trying to implement:

import shapefile

# Create a reader instance for our US Roads shapefile
r = shapefile.Reader("roadtrl020") 

# Create a writer instance copying the reader's shapefile type
w = shapefile.Writer(r.shapeType)

# Copy the database fields to the writer
w.fields = list(r.fields)

# Our selection box that contains Puerto Rico
xmin = -67.5
xmax = -65.0
ymin = 17.8
ymax = 18.6

# Iterate through the shapes and attributes at the same time
for road in r.iterShapeRecords():
    # Shape geometry
    geom = road.shape
    # Database attributes 
    rec = road.record
    # Get the bounding box of the shape (a single road)
    sxmin, symin, sxmax, symax = geom.bbox
    # Compare it to our Puerto Rico bounding box.
    # go to the next road as soon as a coordinate is outside the box
    if sxmin <  xmin: continue
    elif sxmax > xmax: continue
    elif symin < ymin: continue
    elif symax > ymax: continue
    # Road is inside our selection box.
    # Add it to the new shapefile
    w._shapes.append(geom)
    w.records.append(rec)
    
# Save the new shapefile! (.shp, .shx, .dbf)
w.save("Puerto_Rico_Roads")

I am getting the following error: AttributeError: 'Writer' object has no attribute '_shapes'

I think this issue may be connected to the issue discussed here: https://github.com/GeospatialPython/pyshp/issues/113

However when I convert the lines:

w._shapes.append(geom)
w.records.append(rec)

to something like:

w.shape(geom)
w.record(rec)

This does not help. And I am unsure of how to append multiple records or shapes to the writer.

Answer 1

I would suggest using the geopandas and shapely libraries for this operation, this will require much less code to achieve the same operation. You can create a Polygon for your BBOX to which your shapefile has to be clipped.

import geopandas as gpd
from shapely.geometry import Polygon

shp = gpd.read_file("roadtrl020.shp")

xmin, xmax, ymin, ymax = -67.5, -65.0, 17.8, 18.6

# Selection box
selection_box = Polygon([(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)])

clipped_shp = gpd.clip(shp, selection_box)

# Save clipped shapefile
clipped_shp.to_file("clipped_roads.shp", driver="ESRI Shapefile")

Answer 2

So i found that by making the conversion as I suggested above, but adding * to the last line (for unpacking the tuple), it seems to be working correctly now.

Corrected code:

# Create a reader instance for our US Roads shapefile
r = shapefile.Reader(loc)

# Create a writer instance copying the reader's shapefile type
save = saveloc + "\\" + "outputshp"
w = shapefile.Writer(save, r.shapeType)

# Copy the database fields to the writer
w.fields = list(r.fields)
print(r.fields)

# Our selection box that contains Puerto Rico
#personal definition of coordinates (top left and bottom right coords)
xmin, xmax, ymin, ymax = coord1[1], coord2[1], coord2[0], coord1[0]

# Iterate through the shapes and attributes at the same time
for road in r.iterShapeRecords():
    # Shape geometry
    geom = road.shape
    # Database attributes
    rec = road.record
    # Get the bounding box of the shape (a single road)
    sxmin, symin, sxmax, symax = geom.bbox
    # Compare it to our Puerto Rico bounding box.
    # go to the next road as soon as a coordinate is outside the box
    if sxmin < xmin:
        continue
    elif sxmax > xmax:
        continue
    elif symin < ymin:
        continue
    elif symax > ymax:
        continue
    # Road is inside our selection box.
    # Add it to the new shapefile

    w.shape(geom)
    w.record(*rec)

Thanks to this page for help: Creating new shapefile using pyshp