Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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_Transforms does not make a noticable difference. I see a performance of 1 point processed per 10-15 seconds.

share|improve this question

3 Answers 3

up vote 2 down vote accepted

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.

share|improve this answer
    
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?

share|improve this answer
1  
<#> 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
3  
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);
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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