Spatial join in Python without QGIS, ArcGIS, PostGIS

Question

I want to follow that post More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

My target is to count points in polygons and export to new shapefiles the polygons that have a specific number of points on the count or score.

import fiona
from shapely.geometry import shape
polygons = [pol for pol in fiona.open('poly.shp')]
points = [pt for pt in fiona.open('point.shp')]
# attributes of the polygons
for poly in polygons:
   print poly['properties']
for pt in points:
    print i['properties']
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     #iterate through polygons
     for j, poly in enumerate(polygons):
        if point.within(shape(poly['geometry'])):
             # sum of attributes values
             polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

If I run this code with my shapefiles I take keyerror 'score', that is because on my files I don't have a score field.

How to fix that without changing my schema from the input shapefiles?

I don't want to import new field in my files. Any ideas?

score 3 · Answer 1 · edited Apr 13 '17 at 12:34

The script is for a specific shapefile with a score field. In your case, why not simply use a dictionary of polygons as in How to find out which polygon contains the most points using Python OGR?

import fiona
from shapely.geometry import shape
polygons = [pol for pol in fiona.open('polys.shp')]
points = [pt for pt in fiona.open('points.shp')]

 totpts = {}
 for i, poly in enumerate(polygons) :
     totpts[i] = 0
     for j, pt in enumerate(points):
           point = shape(pt['geometry'])
           if shape(poly['geometry']).contains(point):
               totpts[i] += 1

print totpts
{0: 3, 1: 2, 2: 6}

If you open the polygons shapefile now

polys = fiona.open('polys.shp')
# show the third polygon
polys[totpts.keys()[2]]
{'geometry': {'type': 'Polygon', 'coordinates': [[(151.05805243445693, -185.31460674157304), (167.9943820224719, -167.48689138576782), (215.68352059925093, -180.85767790262173), (251.33895131086143, -218.29588014981277), (249.11048689138573, -244.14606741573033), (189.83333333333331, -253.95131086142322), (144.37265917602994, -221.8614232209738), (151.05805243445693, -185.31460674157304)]]}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'id', 3)])}

Spatial join in Python without QGIS, ArcGIS, PostGIS

1 Answers1