3

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.

3

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)
1
  • In the QGIS console, I tried the first solution by providing the absolute paths to the directories. It did not output any parcels. I can now get the Intersection and other built-in functions in QGIS to work (before I wasn't letting it load sufficiently, rookie mistake!!). Do you have any thoughts as to why it isn't working. I'm using the same exact code as before, except now I'm using a larger shapefile (representative of all NYC). Jun 30 '15 at 14:57

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.