3

My points table has 246,338 features and the polygon table has 43,957 features. I am relating the values of one attribute from the polygons to update a column in the points table.

This question is similar to Speed up point sampling with ST_Value function in PostGIS, raster/vector overlay but that question specifically addresses point sampling of raster, not polygons, using ST_Value, and the optimum raster tile storage size.

Both tables have had maintenance performed and have spatial indexes that are updated. The polygons are records of fertilizer application rates applied to a field and are fairly messy with plenty of overlaps:

close up of the fertilizer polygons:

The polygons table has been checked for invalid geometries and those have been fixed. Running ST_NPoints tells me there are an average of only 5 vertices per polygon so ST_Simplify is not necessary. The average polygon size is 112 m2 and range from 0 (only 60 of them) to 826 m2. I cannot dissolve (ST_Union) polygons on anything because I need the unique numeric fertilizer values from those individual polygons.

I've tried two queries in PostGIS, but have abandoned them after running them for a few hours:

Query 1:

update points
set fert_rate = rate
from points a, polys b
where st_intersects(a.geom, b.geom);

Query 2:

with pnt_sample as (
  select rate
  from points a, polys b
  where st_intersects(a.geom, b.geom)
)

update points
set fert_rate = pnt_sample.rate
from pnt_sample;

Having QGIS read the tables from the DB and run the 'point sampling' tool takes only a couple of minutes. However this is not the optimal solution for many reasons, particularly because the output must be saved as a shapefile, headers are truncated to 10 characters, and then must be imported back into the DB, header names fixed, etc.

My questions are:

  1. What is QGIS doing in its algorithm that I don't have in my query?
  2. How can I incorporate that advantage into my approach so I don't have to go outside of my database?

I would like to keep any suggestions/solutions within PostGIS/QGIS.

I'm using PostGIS 2.4.1 and QGIS 2.18.14 on a Windows 10 machine.

EDIT: In response to John Powell and Thingamubob's comments, I ran EXPLAIN ANALYZE on the two queries above. Both tables had vacuum/analyze performed just prior to running EXPLAIN ANALYZE. Here is the output for the result of the first query:

query 1:

And here is the output from the second query, the one using the CTE:

query 2:

Each of these EXPLAIN ANALYZE queries took about an hour to run.

  • For this kind of question, it is essential to include the output of EXPLAIN. Otherwise it is guesswork. It very much sounds like the spatial index is not being used, though, given the performance differences between QGIS and Postgis. – John Powell May 15 '18 at 16:56
  • In QGIS, Processing can output another format than shapefile. Especially in QGIS 3, it might be even quicker because of the port to C++ of Processing. – etrimaille May 16 '18 at 0:34
  • @Gustry, can you please explain how this is done? I know some tools allow for output of a temporary file. This would be ideal, as I could easily push that back to my postgis DB without the intermediate troubles of saving as a shapefile. This particular 'point sampling tool' allows for only one output option, which is a shapefile. – pdavis May 16 '18 at 15:05
  • your queries are creating multiple sets of records per point for each polygon it intersects and does not filter for a unique return value for the SET command to use; i.e. it creates exhaustive overhead with unoptimized records, plus wrong results. given the proper setup, I run similar updates on twice the rows on lesser equipped computers in less than a minute. I, too, assume the indexes are not used properly (did you e.g. run VACUUM ANALYZE <table> on both tables?). to be able to validate, though, you have to provide the EXPLAIN ANALYZE <query>result... – ThingumaBob May 17 '18 at 20:43
  • @JohnPowellakaBarça I've edited the question to include the output of EXPLAIN ANALYZE. Any insights would be greatly appreciated. – pdavis May 21 '18 at 14:25
1

Does it work any better with an inner join rather than a full outer join? I expect an inner join would make better use of the indexes, where the full outer join would assemble every possible combination and then evaluate the results.

(This was edited with a presumption that you have a key column called "id" with a single-field primary key or a unique index)

with pnt_sample as (
  select a.id, b.rate
  from points a
  JOIN polys b
  ON st_intersects(a.geom, b.geom)
)

update points p
set p.fert_rate = ps.rate
from pnt_sample ps
where p.id = ps.id;

Alternatively:

update points
set a.fert_rate = b.rate
from points a JOIN polys b
ON st_intersects(a.geom, b.geom);

Related: Acquiring ArcGIS-like speed in Postgis and https://stackoverflow.com/questions/17946221/sql-join-and-different-types-of-joins

If polygons overlap, then choosing the correct value may require some grouping.

  • Nate, your suggestion is still running after 2+ hours; however I noticed if I simply 'select' instead of 'update' using your query or either of mine, it returns the expected values in 11 seconds in all cases. So it seems the intersection itself is not the issue, but something to do with the update process? – pdavis May 11 '18 at 20:57
  • New update - it ran for 4 or 5 hours and finished correctly. What's strange is that it took this long, but the message when it finished said 'query returned successfully in 12 min.'. So I wonder why the discrepancy, but that's a whole different question which I'll ask the SE experts on Monday. – pdavis May 11 '18 at 23:02
  • this query will try to set the fert_rate column to each entry in pnt_sample; you need to add a filter to get only one row as the SET parameter. can you be sure ST_Intersects does always return one point per polygon in your Query 1? – ThingumaBob May 12 '18 at 8:11
  • @ThingumaBob. Good point. I had just grabbed one of the queries and added a join. I put the join on the other query since I'm not sure which tables fert_rate and rate are in. – Nate Wanner May 12 '18 at 23:19
  • @ThingumaBob and Nate, you both have good points. First, the 'fert_rate' is in the points file and 'rate' is the polygon value being assigned to the points. There are many points per polygon and each point needs the rate from the poly. Nate, in regards to your edit: Overlaps are not a concern because after point sampling I actually eliminate points from overlapping areas, having extracted and isolating those locations in a separate table. But I want to start with the full point set before filtering based on other criteria. I'm going to run the query in your edit and see what happens. – pdavis May 15 '18 at 14:56

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.