3

I have the following datasets stored in PostGIS 9.5 with all the relevent indexes built.

  1. Collected Points - over 4 million points
  2. Polygons - almost 5000 multi polygons (and they are multi polygons)
  3. Road Lines - over 5 millions linestrings
  4. Access Points - 68k point features

For each Collected Point I need to test whether it is in a Polygon, if it is then select the nearest Access Point, if it is not in the Polygon, then select the closest point from the lines.

I have tried a number of methods and the closest I can get to the right SQL was this:

CREATE TABLE test AS 
WITH polygon AS (
    SELECT ST_Union(geom) as geom FROM polygons
)

SELECT
cp.*,
CASE
    WHEN ST_Intersects(cp.geom, polygon.geom) THEN (Select ap.geom FROM access_points ap ORDER BY ap.geom <-> cp.geom LIMIT 1) 
    ELSE ST_CLosestPoint(l.geom, cp.geom)
END AS closest_geom
FROM
polygon,
collected_points cp,
lines l

this code is running but even on a sample of 150 Collected Points the query has been running for a long time.

Maybe I am trying to do too much in a query or need to use more subqueries or CTS statements.

Can someone please suggest a better solution

  • this is going to likely be real slow as bounding box on a linestring may not be useful for closest point type ops and that function won't do any good. – Evan Carroll Jan 22 '18 at 21:07
  • For instance if [-100,+100] is the same <-> to 0,0 as -1,+1 despite both of the points in [-100,+100] being massively further away. – Evan Carroll Jan 22 '18 at 21:10
  • 1
    By unioning all polygons you will get no help from the index. Take a look at this answer for a better approach to find non intersecting gis.stackexchange.com/questions/4502/… – Nicklas Avén Jan 22 '18 at 22:19
1

First, to get decent performance, decompose your linestring to points.

CREATE TABLE pts AS
  SELECT (dp).geom AS lnpt, l.geom
  FROM lines AS l
  CROSS JOIN LATERAL ST_DumpPoints(geom) AS dp;

CREATE INDEX ON pts USING gist (lnpt);

Then run something like this,

SELECT cp.geom, COALESCE(ap.geom, l.geom) AS wanted
FROM collected_points AS cp
LEFT OUTER JOIN LATERAL (
  SELECT ap.geom
  FROM access_points
  WHERE EXISTS (
    SELECT 1
    FROM polygon AS p
    WHERE ST_Intersects( cp.geom, p.geom )
  )
  ORDER BY ap.geom <=> cp.geom
  OFFSET 0 ROWS
  FETCH FIRST 1 ROW ONLY
) AS ap
LEFT OUTER JOIN (
  SELECT l.geom
  FROM pts
  ORDER BY ap.geom <=> pts.lnpt
  OFFSET 0 ROWS
  FETCH FIRST 1 ROW ONLY
) AS l
  ON ap IS NULL;
  • thanks for your help. Just trying it out but getting a syntax error ERROR: syntax error at or near ")" LINE 14: ) AS ap, have tried editing but unsure what is causing it – tjmgis Jan 22 '18 at 21:43
  • @tjmgis fixed now – Evan Carroll Jan 22 '18 at 22:03
  • 1
    Be aware that you will get a different answer if you decompose the points. Tehe closest point on the line might be between the vertex points. – Nicklas Avén Jan 22 '18 at 22:14
  • If that's what your looking for you're going to have to calculate the closest point on every line and see which one is closer for every point in the result set. That is to say, it's simply going to be slow (afaik) – Evan Carroll Jan 22 '18 at 23:47
  • @EvanCarroll thank you for you help. Never used a JOIN LATERAL and it is fantastic. – tjmgis Jan 23 '18 at 9:38

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.