This web page shows how to use a Bounding Box point-in-polygon search before the more expensive Within spatial query of Shapely.

http://rexdouglass.com/fast-spatial-joins-in-python-with-a-spatial-index/

Use Rtree as an index to perform the much faster joins, then Shapely to do the spatial predicates to determine if a point is actually within a polygon. If done properly, this can be faster than most other GISes.

See examples here or here.

The second part of your question concerning 'SUM', use a dict object to accumulate populations using a polygon id as the key. Although, this type of thing is done much more nicely with PostGIS.

Fiona returns Python dictionaries and you can not use poly['properties']['score']) += point['properties']['score']) with a dictionary.

Example of summing attributes using the references given by Mike T:

# read the shapefiles 
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'] 
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])

# attributes of the points
for pt in points:
    print i['properties']
 OrderedDict([(u'score', 1)]) 
 .... # (same for the 8 points)

Now, we can use two methods, with or without a spatial index:

1. without

   iterate through points

   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']

2. with a R-tree index (you can use pyrtree or rtree)

   Create the R-tree index and store the features in it (bounding box)

   from rtree import index idx = index.Index() for pos, poly in enumerate(polygons): idx.insert(pos, shape(poly['geometry']).bounds)

   #iterate through points for i,pt in enumerate(points): point = shape(pt['geometry'])

   iterate through spatial index

   for j in idx.intersection(point.coords[0]): if point.within(shape(polygons[j]['geometry'])): polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

Result with the two solutions:

for poly in polygons:
   print poly['properties']    
 OrderedDict([(u'score', 2)]) # 2 points in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon

What is the difference ?

• Without the index, you must iterate through all the geometries (polygons and points).
• With a bounding spatial index (Spatial Index RTree), you iterate only through the geometries which have a chance to intersect with your current geometry ('filter' which can save a considerable amount of calculations and time...)
• but a Spatial Index is not a magic wand. When a very large part of the dataset has to be retrieved, a Spatial Index cannot give any speed benefit.

After:

schema = fiona.open('poly.shp').schema
with fiona.open ('output.shp', 'w', 'ESRI Shapefile', schema) as output:
    for poly in polygons:
        output.write(poly)

To go further, look at Using Rtree Spatial Indexing With OGR, Shapely, Fiona

Additionally - geopandas now optionally includes rtree as a dependency, see the github repo

So, instead of following all the (very nice) code above, you could simply do something like:

import geopandas
from geopandas.tools import sjoin
point = geopandas.GeoDataFrame.from_file('point.shp') # or geojson etc
poly = geopandas.GeoDataFrame.from_file('poly.shp')
pointInPolys = sjoin(point, poly, how='left')
pointSumByPoly = pointInPolys.groupby('PolyGroupByField')['fields', 'in', 'grouped', 'output'].agg(['sum'])

To get this snazzy functionality be sure to install the C-library libspatialindex first

EDIT: corrected package imports