19

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.

Improve this question
16

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.

Improve this answer
1
  • 6
    Very good point on time spent realizing the solution vs. computing time. Taking a long time to arrive at an optimal solution is only beneficial if you realize those savings through the optimization (esp. from an employer's point of view). – Sasa Ivetic Jul 25 '11 at 21:20
5

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.

Improve this answer
2
  • 4
    I'd directly load the data using a csv reader and populate an Rtree spatial index. The combination of Rtree and Shapely have an impressive performance (much better than PostGIS; I can't compare to JTS as I don't know Java). – Mike T Jul 25 '11 at 22:59
  • 2
    Good idea provided you don't need to store all 1b points in memory at once. At a minimum of 16 bytes per point (X/Y), you're looking at 16GB worth of data. If Rtree will build the index on local storage, then it will definitely improve performance. Importing 1b points to a single shapefile will not work either. OGR specs state shapefiles are limited to 8GB (4GB recommended). A single point shape uses 20 bytes. – Sasa Ivetic Jul 26 '11 at 13:34
4

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.

Improve this answer
3

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...

Improve this answer
1
  • 1
    If you have the data in PostGIS then GeoTools can make use of gist indexes etc. – Ian Turton♦ Jul 26 '11 at 10:46
3

Use Spatialite.

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.

Improve this answer
3

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
Improve this answer
2
  • 2
    Doesn't that loop through 1 bn points 5000 times (once for each polygon)? – underdark♦ Jul 25 '11 at 20:32
  • A spatial index is an absolute must. I mentioned Rtree before, and I'll mention it again! – Mike T Jul 26 '11 at 0:52
-1

might try csv2shp csv2shp

curious to know what industry the billion point CSV is in?

Improve this answer
1

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.