Problem definition is quite simple: I need to achieve point-in-polygon detection for:

  • billions of lat/lon points
  • over thousands of polygons

I need to be able, for every point, to determine if it is in one or several of the polygons available (in the thousands). Currently we are using PostGIS with ST_Contains but the whole things needs 3 good entire days (roughly 80 hours) to compute from start to finish.

Is there any solution (free or not) that would offer significant speed improvement for this workload over PostGIS? It can be anything from an analytical database (Exasol is one example that comes to mind, however although hugely faster than Postgres as a database, I'm not sure if it's any faster for geographical computations) to a dedicated C/C++/Rust piece of code.

  • 5
    sub divide the polygons – Ian Turton Jun 5 '20 at 16:19
  • 2
    And I would ignore the worry about rounding errors (no one cares about nanometers) – Ian Turton Jun 5 '20 at 19:03
  • 2
    Splitting complex polygons with many vertices into smaller ones with less vertices is a relevant approach everywhere so consider PostGIS and ST_Subdivide as an example. Read a study with SpatiaLite gaia-gis.it/spatialite-3.0.0-BETA1/WorldBorders.pdf – user30184 Jun 8 '20 at 7:23
  • 1
    Sounds interesting. This one is quite fast as well youtu.be/_r4IqjGqGEY but the methods are not open. – user30184 Jun 8 '20 at 19:01
  • 1
    Raster operations are sometimes much faster than playing with geometries. If you could convert your polygons into a raster in a clever way you migh be able to do point-in-polygon check by reading the pixel value. You should use binary encoding: first polygon=1, next=2, then 4, 8 and so on. Increase the value for overlapping polygons. If the pixel value is 10 you know that it comes from 2+8. – user30184 Jun 9 '20 at 7:57

You need some kind of precomputed index structure.

In this case, a grid might be useful. They are usually not used for general-purpose applications because it would be hard to choose the proper grid cell size, but your polygons should have known characteristics.

First, choose a fixed grid. For each grid cell, store two lists: polygons that completely cover the grid cell, and polygons that partically cover the grid cell. Make the grid cells small enough so that almost all of them are completely inside or outside all polygons.

To check a point, compute the grid cell number from its coordinates, and then read that cell's lists. You need to do raycasting only for partially covered cells (and for complex polygons, it is useful to precompute the intersection with the cell).

You need to write some code yourself, but at least the initial computation of the grid index can be done with PostGIS.

  • Brilliant man. Bravo. – Jivan Jun 9 '20 at 15:43

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.