Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
I have a 1 billion point CSV file and a shapefile with around 5,000 polygons. What would be the fastest way to spatially join points and polygons? For each point, I need to get the containing polygon id. (Polygons don't overlap.)
Usually, I'd load both data sets into PostGIS. Is there a faster way to get the work done?
I'm looking for an open-source solution.
If "fastest" includes the amount of your time that is spent, the solution will depend on what software you are comfortable with and can use expeditiously. The following remarks consequently focus on ideas for achieving the fastest possible computing times.
If you use a canned program, almost surely the best you can do is pre-process the polygons to set up a point-in-polygon data structure, such as a K-D tree or quadtree, whose performance will typically be O(log(V)*(N+V)) where V is the total number of vertices in the polygons and N is the number of points, because the data structure will take at least O(log(V)*V) effort to create and then will have to be probed for each point at a per-point cost O(log(V)).
You can do substantially better by first gridding the polygons, exploiting the assumption of no overlaps. Each grid cell is either entirely in a polygon interior (including the interior of the "universal polygon"), in which case label the cell with the polygon's id, or else it contains one or more polygon edges. The cost of this rasterization, equal to the number of grid cells referenced while rasterizing all the edges, is O(V/c) where c is the size of a cell, but the implicit constant in the big-O notation is small.
(One beauty of this approach is that you can exploit standard graphics routines. For example, if you have a system that (a) will draw the polygons on a virtual screen using (b) a distinct color for each polygon and (c) allows you to read the color of any pixel you care to address, you've got it made.)
With this grid in place, pre-screen the points by computing the cell containing each point (a O(1) operation requiring only a few clocks). Unless the points are clustered around the polygon boundaries, this will typically leave only about O(c) points with ambiguous results. The total cost of building the grid and pre-screening is therefore O(V/c + 1/c^2) + O(N). You have to use some other method (such as any of those recommended so far) to process the remaining points (that is, those which are close to polygon boundaries), at a cost of O(log(V) * N * c).
As c gets smaller, fewer and fewer points will be in the same grid cell with an edge and therefore fewer and fewer will require the subsequent O(log(V)) processing. Acting against this is the need to store O(1/c^2) grid cells and to spend O(V/c + 1/c^2) time rasterizing the polygons. Therefore there will be an optimal grid size c. Using it, the total computational cost is O(log(V) * N) but the implicit constant is typically way smaller than using the canned procedures, due to the O(N) speed of the pre-screening.
20 years ago I tested this approach (using uniformly spaced points throughout England and offshore and exploiting a relatively crude grid of around 400K cells offered by the video buffers of the time) and obtained two orders of magnitude speedup compared to the best published algorithm I could find. Even when the polygons are small and simple (like triangles), you are virtually assured of an order of magnitude speedup.
In my experience the computation was so fast that the entire operation was limited by data I/O speeds, not by the CPU. Anticipating that I/O might be the bottleneck, you would achieve the very fastest results by storing the points in as compressed a format as possible to minimize the data reading times. Also give some thought to how the results should be stored, so that you can limit disk writes.
For my part, I would probably load CSV data into a shp file and then write a python script using shapefile and shapely to get the containing polygon id and update the field value.
I don't know whether geotools and JTS is faster than shapefile/shapely ... Have no time to test it!
edit : By the way, the csv conversion to shapefile format is probably not required, since values could easily be formatted to be tested with spatial objects from your polygon shapefile.
I ended up converting the polygons to a raster and sampling it at the point positions. Since my polygons didn't overlap and high accuracy wasn't necessary (polygons represented land-use classes and their borders were considered rather uncertain anyway) this was the most time-efficient solution I could come up with.
I would quickly write a small java program based on the shapefile reader of geotools and the operation contains of JTS. I do not know how fast it can be...
Download the GUI. You can open both the Shapefile and CSV as virtual tables. This means that you don't actually import them into the database but they appear as tables and you can quicky join and query them any way you like.
You can do it fairly quickly using OGR in C/C++/Python (Python should be the slowest of the 3). Loop through all the polygons and set a filter on the points, loop through the filtered points and you'll know that each of the points that you loop through will belong to the current polygon. Here's sample code in python using OGR that will loop through the polygons and filter points accordingly. C/C++ code will look quite similar to this, and I would imagine you'll get a significant speed increase vs. python. You'll have to add a few lines of code to update the CSV as you go along:
from osgeo import ogr
from osgeo.gdalconst import *
inPolyDS = ogr.Open("winnipeg.shp", GA_ReadOnly)
inPolyLayer = inPolyDS.GetLayer(0)
inPointDS = ogr.Open("busstops.vrt", GA_ReadOnly)
inPointLayer = inPointDS.GetLayerByName("busstops")
inPolyFeat = inPolyLayer.GetNextFeature()
while inPolyFeat is not None:
inPtFeat = inPointLayer.GetNextFeature()
while inPtFeat is not None:
ptGeom = inPtFeat.GetGeometryRef()
# Do work here...
inPtFeat = inPointLayer.GetNextFeature()
inPolyFeat = inPolyLayer.GetNextFeature()
VRT File (busstops.vrt):
<OGRVRTDataSource>
<OGRVRTLayer name="busstops">
<SrcDataSource>busstops.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>WGS84</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y" reportSrcColumn="FALSE" />
</OGRVRTLayer>
</OGRVRTDataSource>
CSV File (busstops.csv):
FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey
CSVT File (busstops.csvt, OGR needs it to identify column types, otherwise it will not perform the spatial filter):
Integer,Real,Real,String
might try csv2shp csv2shp
curious to know what industry the billion point CSV is in?
