2

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.

1
  • You don't need elif when the if has a continue or break or return and your# Road is inside our selection box. comment is incorrect -- You have only proved that the envelopes overlap. This might be enough for your purpose, but it would include many, if not all, roads in PA, NJ, CT, MA, NH, and VT if you tried it with a NY polygon.
    – Vince
    Aug 1 '20 at 13:16
4

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")
3
  • This is much more succinct, thank you. However, the read_file() function seems to be much much slower than the shapefile library alternative? Does this sound right?
    – ushham
    Aug 2 '20 at 6:55
  • Yes this sounds right for large shapefiles. Aug 2 '20 at 9:50
  • Ok, thank you this makes sense. As my shapefile is large I am finding performance of my below solution is much better, though i acknowledge it is not pretty.
    – ushham
    Aug 2 '20 at 10:02
1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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