I need to attribute a geographical area (based on shapefile polygons) to 8 million GPS points. There are 405 polygons. I currently have the code below, and a quick extrapolation tells me it will take 84 hours to compute this.

I am new to GIS, but intuitively, I feel that there has to be a smarter way around this than, for each point, randomly testing each polygon with a point-in-polygon function until there is match. For example, simply putting the points and polygons in 2 groups (such as "East" and "West") simply based on their easternmost/westernmost longitude may divide by two the number of operations to perform by almost 2.

Are there indeed such much efficient algorithms? Are any of them easily implementable in python?

from shapely.geometry import shape, Point, Polygon
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

m = Basemap(projection='merc',llcrnrlat=15,urcrnrlat=55,llcrnrlon=70,urcrnrlon=140,lat_ts=20,resolution='c')
m.readshapefile(china_shapefile, 'prefectures', drawbounds=True)

print len(m.prefectures)

#get data



for coordinate in coordinates:
    for info,shape in zip(m.prefectures_info,m.prefectures):
        poly = Polygon(shape)
        #print info["NAME_1"],info["NAME_2"]
        xpt,ypt = m(coordinate[1],coordinate[0])
        point1 = Point(xpt,ypt)
        if point1.within(poly):
            print coordinate,info["NAME_1"],info["NAME_2"]#,info["NAME_3"]
  • 2
    Others will no doubt provide python code but as an algorithm I would think: For each polygon, Get extent, Filter sorted points to only include those that have a chance of being within this extent, Check if within, Next polygon. – David W Dec 1 '15 at 12:36
  • 2
    This would be simple and relatively fast if you were using GIS software, but you're just using Python. Sequential access isn't "random" (random would be much worse), but you're still pushing the short end of the lever -- you want to load the points in spatially indexed format, then perform 405 queries (not 405 sequential searches, but indexed queries) for point within polygon. Unfortunately, there are dozens of potential solutions, making this question too broad as written. – Vince Dec 1 '15 at 12:38
  • 5
    Any time you have a million things, a database would be an appropriate tool for manipulation, so I would recommend you start by loading your features in PostgreSQL and do a simple join -- it should take less than 84 hours to prepare to do this minutes-long query. – Vince Dec 1 '15 at 12:39
  • 1
    This previous question gives very good answers: gis.stackexchange.com/questions/102933/… Also see this post for more on rtree snorf.net/blog/2014/05/12/using-rtree-spatial-indexing-with-ogr – blord-castillo Dec 1 '15 at 14:11
  • 2
    If at all possible, and if you are looking for the most efficient way to do this, download QGIS(www.qgis.org) and then search on GIS SE how to achieve the results you want. What you are trying to accomplish without the use of GIS software is possible but as far as processing time and efficiency I can't think of anything better than actually using a piece of software designed for spatial analysis. – GeoJohn Dec 1 '15 at 14:53

As blord-castillo says, you can use a Spatial Index. They are many modules in Python for that (PiPy: spatial index, R-tree, Quad-Ttree, Kd-tree, Interval-tree, ...).

The most used method for indexing spatial data is the R-tree index (also used by QGIS or PyQGIS as stated by GeoJohn, see geoprocessing across multiple vector layers QGIS2 for example)

The fastest is Rtree (rtree python polygon index) because it is a simple Python wrapper of the C library libspatialindex. It is used in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc, for example (with Fiona, Shapely and Rtree).

Rather than using Pandas, you should use GeoPandas that includes optionally rtree as a dependency (spatial joins, Spatial indexes for geo data frames and series or More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc with GeoPandas)

This reference gives you 2 solutions for your problem (spatial join between points and polygons) with the limitations of the procedure:

  • Without an 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.

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.