4

I am trying to create a new table of linestrings with their existing attribution that are "different" to some polygons and am struggling to get a working PostGIS query.

I have a table of linestrings (~130,000) that intersect polygons (~13,0000 which are valid), I created this table using ST_Intersects.

Orignal Data

The query I originally tried to use was

CREATE TABLE new_lines AS SELECT
l.id,
l.value
ST_Difference(l.geom, p.geom)::geom(LineString, 27700) as geom
FROM
lines l,
polygons p;

This query ran for several hours and then failed due to GeometryCollection or MultiStrings, so had to amend my query to:

CREATE TABLE new_lines AS SELECT
l.id,
l.value,
ST_CollectionExtract(ST_Multi(ST_Difference(t.wkb_geometry, p.geom)), 2)::geometry(MultiLineString, 27700) as geom
FROM
lines l,
polygons p;

Whilst that was running I thought I would try the identical task in QGIS using "difference" and that successfully completes in 42 minutes giving me the results I expected.

QGIS Difference

The problem is that the Shapefile it creates is over the 2Gb limit hence needing to get a working PostGIS query.

I have also tried the other SQL found here: Trouble using ST_Difference to remove overlapping features and here Difference between two layers in PostGIS. So i tried this:

CREATE TABLE new_lines AS SELECT COALESCE(ST_Difference(l.geom, p.geom), l.geom) As geom FROM lines l LEFT JOIN polygons p;

But both give me some good results but they do not seem to apply ST_difference to all the lines which is really confusing.

Failed PostGIS query

Am sure this should be easy to do but cannot understand why my first query is taking so long and two why these other solutions only partly work.

EDIT1:

So I amended my query to the following:

CREATE TABLE new_lines AS SELECT
row_number() over() AS gid, g.*
FROM
(SELECT 
l.value 
ST_CollectionExtract(ST_Multi(ST_Difference(a.geom, b.geom)), 2)::geometry(MultiLineString, 27700) as geom
FROM lines a, polygons b
where st_intersects(a.geom,b.geom)) g;

Which successfully ran in around 20 minutes, however, not all the lines are clipped to the polygon boundary. Some are some are not.

I have checked that all geometries are valid and that there are no duplicate features in either dataset.

So am very confused.

Improve this question
12
  • First, I assume you have spatial indexes in place? Second, I have never tried coalesce in that context, though it is a good idea, but you have reversed dbaston's approach, the point of which is to select the geometry if there is no result from difference.
    – John Powell
    May 11, 2016 at 19:13
  • Yes got indexes added.
    – tjmgis
    May 11, 2016 at 19:19
  • 1
    Where is your ST_Intersects(l.geom, p.geom) clause? Without that, you have basically a cartesian join between the two tables, which would certainly explain the speed, if not the results.
    – John Powell
    May 11, 2016 at 19:25
  • 2
    Are you saying that all 130,000 lines intersect all 13,000 polygons? Or simply that every line intersects at least one polygon. If the latter, you really need ST_Intersects.
    – John Powell
    May 11, 2016 at 19:38
  • 1
    @tjmgis They are clipped at polygon boundary, the thing is that if a linestring intersects with two polygons you gonna get two linestrings, one clipped at polygon A, one clipped at polygon B.
    – Jakub Kania
    May 14, 2016 at 17:58

1 Answer 1

Reset to default
2

The query should use ST_UNION to make sure that the linestrings are clipped at all the points they intersect and a lateral join so it would use the proper indexes.

SELECT row_number() over() AS gid,
ST_CollectionExtract(ST_Multi(ST_Difference(a.geom, b.geom)), 2)::geometry(MultiLineString, 27700) as geom
FROM lines a
  ,LATERAL (
    SELECT ST_UNION(polygons.geom)
    FROM polygons
    WHERE ST_Intersects(a.geom,polygons.geom)
  ) AS b
Improve this answer
3
  • thanks so much. Yes the missing thing was to include an ST_Union. However, using a LATERAL seemed to take longer than just a standard query. Thanks again
    – tjmgis
    May 15, 2016 at 18:27
  • Hmm...Well, postgis indexes are very moody. I'd have to see the explain analyze and the list of indexes to say why it was slower. But oh well, it was just a shot in the dark.
    – Jakub Kania
    May 16, 2016 at 6:44
  • I would use ST_Collect instead of ST_Union, which is much faster.
    – caiohamamura
    Aug 4, 2017 at 21:45

Your Answer

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

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