python - Intersect Shapefiles using Shapely

What I'd like to try and do is intersect some point buffers with census tracts, maintain their ID's and the tracts attributes.

I found a sample from http://www.macwright.org/2012/10/31/gis-with-python-shapely-fiona.html but can't exactly substitute the cascaded union for an intersect. I'm also having trouble getting the intersection function to import from shapely.ops.

from shapely.geometry import Point, mapping, shape

from fiona import collection
#from shapely.ops import intersection
import shapely.ops

bufSHP = 'data/h1_buf.shp'
intSHP = 'data/h1_buf_int_ct.shp'
ctSHP  = 'data/nyct2010.shp'
with collection(bufSHP, "r") as input:
    schema = input.schema.copy()
    with collection(intSHP, "w", "ESRI Shapefile", schema) as output:

        shapes = []
        for f in input:
            shapes.append(shape(f['geometry']))
        merged = shapes.intersection(ctSHP)
        output.write({
            'properties': {
                'uid': point['properties']['uid']
                },
            'geometry': mapping(merged)
            })

I've also got the code and sample shapefiles up here https://github.com/nygeog/questions/tree/master/shapely_intersect

---

UPDATE:

So for some reason using Shapely, the larger buffers (1 km) I lose some of the features (census tracts) in the intersect.

Also, since the Shapely responder said the OGR command is good to use I started looking for into that. Anyway, I was able to test cd'ing to the directory and it works great. However, I've been trying to experiment with it so I can loop through the many sized buffers and dif. census years. However, why is it that I cannot get the following script to work. Can I not substitute paths for the folder 'data' in the ogr2ogr response below?

import os
wp = '/Volumes/Echo/GIS/projects/naas/tasks/201411_geoprocessing/data/processing/'
folder = wp+'test'
theGDALcmd = 'ogr2ogr -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM nyct2010 A, hr1000m B WHERE ST_Intersects(A.geometry, B.geometry)" -dialect SQLITE '+folder+' '+folder+' -nln hr1000m_int_ct10_ogr_pytest'

os.system(theGDALcmd)

The error I get is

sh: ogr2ogr: command not found

However, when I go to the terminal, I have access to ogr2ogr.

---

UPDATE2:

Mike T include some code so I can run this ogr2ogr command in python.

Answer

While I'm a big user of both shapely and fiona, I wouldn't go this approach. This is a task of writing an effective SQL statement.

Using ogr2ogr with an SQLITE dialect, you can process this from a command line. Change directory to one before the shapefiles, so that all of the shapefiles are in one directory called data. OGR treats directories of shapefiles as datasets. Now try this command:

ogr2ogr -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM nyct2010 A, h1_buf B WHERE ST_Intersects(A.geometry, B.geometry)" -dialect SQLITE data data -nln h1_buf_int_ct2

This makes h1_buf_int_ct2.shp with the same attributes as A and B, showing the intersections. You could also have done the buffering operation somewhere in there too.

The same SQL statement works with PostGIS too, so you know.

---

Here's how to do the same from within Python using ogr:

from osgeo import ogr
ogr.UseExceptions()

ogr_ds = ogr.Open('/path/to/data', True)  # Windows: r'C:\path\to\data'
SQL = """\
    SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
    FROM nyct2010 A, h1_buf B
    WHERE ST_Intersects(A.geometry, B.geometry);
"""
layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
# copy result back to datasource as a new shapefile
layer2 = ogr_ds.CopyLayer(layer, 'h1_buf_int_ct3')
# save, close

layer = layer2 = ogr_ds = None