Cropping a Shapefile into a Rectangular Box with Python pyshp

Question

Using 'pyshp' I am hoping to reduce the set of US highway shapefiles to only those in Puerto Rico. The idea is very easy since I only need those pieces which fit into the window.

>> import shapefile
>> sf = shapefile.Reader("roadl_usa") 
>> len(sf.records())

193504

Not familiar with the writing object but hoping it's easy since I am not building from scratch. How do I select only those records in the bounding box [-67.5, -65.0] and [17.8, 18.6] and write a new shapefile?

What is the output like? I am wondering if .shx, .txt, .dbs also get exported and what they mean.

---

This script almost works except it leads to IndexError: list index out of range - maybe this is because I am deleting things in the process of iterating?

import shapefile
e = shapefile.Editor("roadl_usa")

prbbox = [-67.5, 17.8,-65.0, 18.6]

L = len(e.shapes())    
for t in range(L):
    if f( e.shape(t).bbox, prbbox) > 0:
        e.delete(t)

e.save('pr-roads')

Answer 1

Non-Python solution

ogr2ogr -f "ESRI Shapefile" pr-roads.shp roadl_usa.shp -clipsrc -67.5  17.8 -65.1  18.6

Answer 2

If you have ogr2ogr modules modules, you can actually run the line given by John Mangual in Python.

This link have same helpful answers regarding to that. Using ogr2ogr to convert GML to shapefile in Python?

Also you can run line with os module as in this link

FWTools ogr2ogr in Python

Answer 3

Your approach is very interesting. The Editor class in PyShp is still a work in progress but I'll look into it further.

Here's some code I just tested using your bounding box on a US roads shapefile. PyShp does export the SHX and the DBF files. Those are the minimum requirements of a shapefile. The SHP file is the geometry, the SHX file is a simple record index for the geometry, and the DBF file contains the attributes.

import shapefile

# Create a reader instance
r = shapefile.Reader("roadl_usa") 

# Create a writer instance
w = shapefile.Writer(r.shapeType)

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

# Our selection box
xmin = -67.5
xmax = -65.0
ymin = 17.8
ymax = 18.6

for road in r.iterShapeRecords():
    # Shape geometry, Database Attributes
    geom = road.shape
    rec = road.record
    # Get the bounding box of the shape
    sxmin, symin, sxmax, symax = geom.bbox
    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")