Overlay Analysis in QGIS (or Preferably Python)

Question

I would like to perform an overlay analysis in Python or QGIS to output a list of each shapefile that falls within certain boundaries. I have unsuccessfully tried to use the Spatial Join, Intersect, and Union commands that were described in this similar GIS Exchange post.

For each parcel described in this shapefile (links to download), I would like to see whether it exists within the regions of this shapefile.

Would you please point me in the right direction on how to do this? If there is a Python library that I can use, I would prefer that method.

Answer 1

1) The easiest solution is to use the processing module in the QGIS Python console:

import processing
processing.runalg("qgis:joinattributesbylocation","BKMapPLUTO.shp","DCP_nyc_freshzoning.shp","['intersects']",0,"sum,mean,min,max,median",0,'result.shp')

2) Without a GIS, you can use Fiona (read and write shapefiles as Python dictionaries) and Shapely (geometries) -> look at More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc and replace within by intersects. I follow here the same procedure.

a) The trivial solution: iterating through all the geometries (very longtime with your shapefiles...).

from shapely.geometry import shape
import fiona
with fiona.open('BKMapPLUTO.shp') as input:
    schema = input.schema
    schema['properties']['NAME'] = 'str:30'
    meta = input.meta
    with fiona.open('DCP_nyc_freshzoning.shp') as regions:
       with fiona.open ('result2.shp', 'w', **meta) as output:
           for feat in input:
               poly1 = shape(feat['geometry'])
               for poly2 in regions:
                   if poly1.intersects(shape(poly2['geometry'])):
                        feat['properties']['NAME']=  poly2['properties']['NAME']
                        output.write(feat)

2) To accelerate the procedure, you can use a boundary spatial index (rtree here)

regions = [pol for pol in fiona.open('DCP_nyc_freshzoning.shp')]
# creation of the spatial index
from rtree import index
idx = index.Index()
for pos, poly in enumerate(regions):
    idx.insert(pos, shape(poly['geometry']).bounds)

with fiona.open('BKMapPLUTO.shp') as input:
    schema = input.schema
    schema['properties']['NAME'] = 'str:30'
    meta = input.meta
    with fiona.open ('result2.shp', 'w', **meta) as output:
        for feat in input:
             poly = shape(feat['geometry'])
             for j in idx.intersection(poly.bounds):
                  if poly.intersects(shape(regions[j]['geometry'])):
                     feat['properties']['NAME']=  regions[j]['properties']['NAME']
                     output.write(feat) 

3) You can do the same with QgsSpatialIndex() in the Python console of QGIS

index = QgsSpatialIndex()
for pos, poly in enumerate(regions):
    index.insertFeature(poly)