Count number of points in multipolygon shapefile using Python

Question

I have a polygon shapefile of the U.S. made up of individual states as their attribute values. In addition, I have arrays storing latitude and longitude values of point events that I am also interested in. Essentially, I would like to 'spatial join' the points and polygons (or perform a check to see which polygon [i.e., state] each point is in), then sum the number of points in each state to find out which state has the most number of 'events'.

I believe the pseudocode would be something like:

Read in US.shp
Read in lat/lon points of events
Loop through each state in the shapefile and find number of points in each state
print 'Here is a list of the number of points in each state: '

Any libraries or syntax would be greatly appreciated.

Based on what I can tell, the OGR library is what I need, but I am having trouble with the syntax:

dsPolygons = ogr.Open('US.shp')  

polygonsLayer = dsPolygons.GetLayer()  


#Iterating all the polygons  
polygonFeature = polygonsLayer.GetNextFeature()  
k=0  
while polygonFeature:
    k = k + 1  
    print  "processing " + polygonFeature.GetField("STATE") + "-" + str(k) + " of " + str(polygonsLayer.GetFeatureCount())  

    geometry = polygonFeature.GetGeometryRef()          

    #Read in some points?
    geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(-122.33,47.09)
    point.AddPoint(-110.11,33.33)
    #geomcol.AddGeometry(point)
    print point.ExportToWkt()
    print point
    numCounts=0.0   
    while pointFeature:  
        if pointFeature.GetGeometryRef().Within(geometry):  
            numCounts = numCounts + 1  
        pointFeature = pointsLayer.GetNextFeature()
    polygonFeature = polygonsLayer.GetNextFeature()
    #Loop through to see how many events in each state

Answer 1

I like the question. I doubt I can give you the best answer, and definitely can't help with OGR, but FWIW I'll tell you what I'm doing right now.

I use GeoPandas, a geospatial extension of pandas. I recommend it — it's high-level and does a lot, giving you everything in Shapely and fiona for free. It is in active development by twitter/@kajord and others.

Here's a version of my working code. It assumes you have everything in shapefiles, but it's easy to generate a geopandas.GeoDataFrame from a list.

import geopandas as gpd

# Read the data.
polygons = gpd.GeoDataFrame.from_file('polygons.shp')
points = gpd.GeoDataFrame.from_file('points.shp')

# Make a copy because I'm going to drop points as I
# assign them to polys, to speed up subsequent search.
pts = points.copy() 

# We're going to keep a list of how many points we find.
pts_in_polys = []

# Loop over polygons with index i.
for i, poly in polygons.iterrows():

    # Keep a list of points in this poly
    pts_in_this_poly = []

    # Now loop over all points with index j.
    for j, pt in pts.iterrows():
        if poly.geometry.contains(pt.geometry):
            # Then it's a hit! Add it to the list,
            # and drop it so we have less hunting.
            pts_in_this_poly.append(pt.geometry)
            pts = pts.drop([j])

    # We could do all sorts, like grab a property of the
    # points, but let's just append the number of them.
    pts_in_polys.append(len(pts_in_this_poly))

# Add the number of points for each poly to the dataframe.
polygons['number of points'] = gpd.GeoSeries(pts_in_polys)

The developer tells me that spatial joins are 'new in the dev version', so if you feel like poking around in there, I'd love to hear how that goes! The main problem with my code is that it's slow.