I need to get a count of neighbors within 100m for each point in a set of 250 million points. The result would have the id of each point accompanied by the number of neighbors around them. Can this be done in postgis in a few hours (the table has an indexed geometry column)? How long should it take and what would be the most optimized query?

  • May I ask why is the timespan so critical? Is this query a one off, or are you running it frequently after modifying the points?
    – Guy
    Jun 9, 2015 at 9:54
  • Since by the current running time calculations, it could easily shoot to about 6-7 days. Which is not a reasonable time span. Also, I will be needing to run this and get results soon as part of a bigger analysis from time to time on different data sets.
    – picmate
    Jun 9, 2015 at 14:12

1 Answer 1


The query will take quite a while, for sure, w/ 250M points (perhaps too long). I'm assuming here that your points are in a metric coordinate system. Looks like this:

SELECT a.gid, count(*)
FROM pts a
JOIN pts b
ON ST_DWithin(a.geom, b.geom, 100)
WHERE a.gid != b.gid
GROUP BY a.gid;
  • 1
    maybe the ST_Envelope (on the start point) before the ST_DWithin postgis.net/docs/manual-1.4/ST_Envelope.html
    – Mapperz
    Jun 8, 2015 at 18:22
  • Thanks for the answer and the comment. I ran a similar query for about 2M points and it took close to 90minutes. Given 250M do you think this could be computed within a day? (I don't think it is going to scale linearly, even then it will be about 6 days correct?) If impossible, do you know of any other software I could use to get the job done?
    – picmate
    Jun 8, 2015 at 19:34
  • I think it probably will scale linearly, as essentially, a query like this will run by doing a table scan on pts a, and then using the spatial index with DWithin to find points from pts b. If you look at EXPLAIN you should see this confirmed. I have had to run even more complex queries on 200M plus polygons and then only way I could get it to run in anything approaching real time was to split up my area into 10 or so rectangular subsets (overlapping by 100m in your case), run these in parallel and then combine the results. Sort of a map-reduce approach, not very pure SQL, but it works. Jun 9, 2015 at 8:37
  • This kind of thing is n*log(n), so not quite linearly, but not quadratically either. I think @johnbarca's suggestion is best, just toss some more cores at it by running several areas in parallel. Jun 9, 2015 at 13:16
  • 1
    @PaulRamsey. You were missing a.gid in the group by so I took the liberty of editing your post. It is n*log(n) because the n is the sequence scan on one side of the self-join and the log(n) is the cost of the spatial index kicking in on ST_DWithin? Assuming that the points are reasonably randomly distributed, then the log(n) is more or less a constant and it scales with n, no? Jun 9, 2015 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.