# In PostGIS, how to count nearby points efficiently?

I have one dataset A with around 6K point geometries and another dataset B with around 600K point geometries.

I am trying to come up with an efficient way to update dataset A with an integer representing the number of objects from dataset B that are within a certain distance (100 meters) from the point in A.

The strategy I came up with did not complete overnight, so I am thinking I could be doing this more efficiently. (Or maybe this is just that expensive an operation?).

I am running this on my X220 i5 laptop with SSD and 8GB RAM, Ubuntu 12.10, PostgreSQL 9.1, Postgis 2.0.

This is the function I use. `nbi.testpoints` is the small dataset, `nbi.bridges` is the bigger dataset. Both tables have spatial indexes on them.

``````CREATE OR REPLACE FUNCTION getnearbycounts() RETURNS VOID
AS
\$\$
DECLARE
node nbi.testpoints%rowtype;
BEGIN
FOR node IN SELECT * FROM nbi.testpoints
LOOP
UPDATE nbi.testpoints SET closenbicount = (SELECT count(1) FROM nbi.bridges WHERE
ST_Transform(node.geom, 3785) <#> ST_Transform(nbi.bridges.geom, 3785) < 100)
WHERE nbi.testpoints.id = node.id;
RAISE NOTICE 'node % done', node.id;
END LOOP;
END;
\$\$
LANGUAGE PLPGSQL;
``````

I am doing the transform to a meter-based coordinate system to be able to give a distance threshold in meters, but replacing that with an approximation in degrees and doing away with the two `ST_Transform`s does not make a noticable difference. I see a performance of 1 point processed per 10-15 seconds.

-

If I understand correctly, the biggest problem you're facing in your current approach is that you're running `st_distance` on a huge dataset. You are running `st_distance` 600K * 6K = 3.6 * 10^9 times (you're computing distances between all features from table A against all features from table B).

Assuming 1ms per query that's 41 days.

Things I would try (in the order of diminished returns):

• build a bounding box 101 meters in each direction around your point of interest. Use `st_within` (or `st_contains` - whichever proves to be faster) to select close enough candidates in relation to this box for distance computation. This avoids running `st_distance` on the cartesian product between your tables and should make things orders of magnitude faster.

• add a new column to both tables with geometry in 3785 projection. If the first option works (and I think it should), this will not bring much improvement.

• if the above options fail consider moving the computation at runtime.

I had a similar project that required computing the status of 1 polygon in relation with 1 million plus other polygons -> although I got the computation down to less then 100 milliseconds updating the entire table still took too much time. Instead of computing for each feature in the table I moved the computation at runtime and increased the answer time by 100 milliseconds which was acceptable. This may or not be an option for you -> I don't know much enough about your project.

Finally, don't forget that you're using a database -> I always make this mistake. PLPGSQL will only get you so far - usually adding a new column (such as the bbox described at point 1), creating a spatial index for it, creating new tables or temporary tables will always run faster then intricate PLPGSQL code.

-
That is great feedback, thanks @diciu. Your first option is what I ended up going with, see my – mvexel Mar 16 '13 at 19:05

It turns out I need to use `ST_DWithin` to narrow down the scope of the search and prevent the query from performing a sequential scan on the entire bigger table on every iteration.

This is what I ended up using and takes only around 20 seconds to run on the entire table:

``````UPDATE nbi.testpoints SET closenbicount = (SELECT count(1) FROM nbi.bridges WHERE
ST_DWithin(nbi.bridges.geom, nbi.testpoints.geom, 0.01) AND
nbi.testpoints.geom <#> nbi.bridges.geom < 0.001);
``````

Much simpler and much more elegant if you ask me. Makes me wonder why <#> does not have this behavior built in?

-
<#> computes distances. Your query first computes the entire result set of distances and then applies a filter (<100) on the result set. What you're asking for is a smart(er) execution plan and not a more optimized <#> operator. – diciu Mar 15 '13 at 16:54
Drop the <#> operator entirely, it's duplicative to what you're getting from ST_DWithin and helps not at all in this context. It's only using in a "nearest N things" context in an ORDER BY clause. – Paul Ramsey Mar 15 '13 at 19:47
Thanks for the comments - I figured out soon enough that indeed operator choice was not the problem here, but the execution plan. – mvexel Mar 16 '13 at 19:00
@PaulRamsey my thought was to narrow down the search results first with ST_DWithin and then perform the final comparison on the subset with <#>, but you are saying it doesn't work like that? Hm. I will need to go back to <#> and ST_DWithin documentation and better understand what I am doing. Thanks for the feedback. – mvexel Mar 16 '13 at 19:02

Maybe a buffer and intersect? Something like

``````SELECT COUNT(*) FROM nbi.bridges WHERE
ST_Intersects(ST_Transform(node.geom, 3785),
ST_Buffer(ST_Transform(nbi.testpoints.geom, 3785), 100);
``````
-