6

I've to compare two polygon layers, 'land_use' and 'land_parcels', in PostGIS. My aim is to create linestrings, where boundaries of 'land_use' polygons do not match boundaries of 'land_parcel' polygons.

enter image description here enter image description here enter image description here

Please note, that a simmilar question has already been asked (see 'How to create lines to visualize differences between polygon features in PostGIS?'). Can anyone help me adapting the query posted by Mike T?

SELECT row_number() OVER () rn, *
FROM (
  SELECT (ST_Dump(ST_SymDifference(ST_Snap(luse, lparcels, tol), ST_Snap(lparcels, luse, tol)))).*
  FROM (
    SELECT ST_Union(DISTINCT land_use_only.geom) luse, ST_Union(DISTINCT land_parcels_only.geom) lparcels, 1e-5 tol
    FROM (
      SELECT ST_Boundary(geom) geom FROM land_use
      EXCEPT SELECT ST_Boundary(geom) geom FROM land_parcels
    ) land_use_only,
    (
      SELECT ST_Boundary(geom) geom FROM land_parcels
      EXCEPT SELECT ST_Boundary(geom) geom FROM land_use
    ) land_parcels_only
  ) s
) s;
1

Wow, that long query is really necessary? What's wrong with:

SELECT 1 gid, ST_Difference(ST_Union(lparcels.geom), ST_Union(luse.geom)) geom
FROM boundary_test.land_use luse,  boundary_test.land_parcels lparcels

Or if you really need single part geometry:

SELECT (dump_geom).path[1],(dump_geom).geom geom
FROM (SELECT (ST_Dump(ST_Difference(ST_Union(lparcels.geom), ST_Union(luse.geom)))) dump_geom
FROM boundary_test.land_use luse,  boundary_test.land_parcels lparcels) t

If you want to snap vertices, you shouldn't ST_Snap twice, you should only snap to the most accurate reference, if land_parcels is more accurate you should:

ST_Snap(luse.geom, lparcels.geom, tolerance)

And shouldn't ST_Snap(lparcels.geom, luse.geom, tolerance). If you snap land_parcels to land_use you will actually include misplace error to land_parcels, assuming land_parcels is more accurate.

EDIT:

If you need to return LineString just do:

SELECT (dump_geom).path[1],(dump_geom).geom geom
FROM (SELECT (ST_Dump(ST_Boundary(ST_Difference(ST_Union(lparcels.geom), ST_Union(luse.geom))))) dump_geom
FROM boundary_test.land_use luse,  boundary_test.land_parcels lparcels) t

Snapping

You should not notice visually the error included by incorrect snap 'cause it may enter very tiny difference polygons. See what I'm talking about snapping.

If you snap A->B:

enter image description here

Or snap B->A:

enter image description here

If you snap both (like you are in the query you posted) you will just switch places:

enter image description here

The right approach is to first snap features to the most accurate mapping and then snap the most accurate mapping to the already snapped features. So if A is more accurate, snap B to A (generate a snapped B), then snap A to snapped B.

Difference vs Symmetric difference

You should also be aware about what each one does.

Symmetric difference doesn't care about the order of operands, because its result is what A has but B doesn't + what B has but A doesn't.

Difference will care about the order of operands. If it is A - B you will get what A has but B doesn't.

Be careful on code you blindly copy and paste. You have to understand what's going on or you may get results you think are ok, but in fact they are not.

FINALLY

I think this is what you wanted:

SELECT 1 gid, ST_Difference(lu.geom, lp.geom) geom
FROM 
    (SELECT ST_Union(ST_Boundary(geom)) geom FROM land_use) lu, 
    (SELECT ST_Union(ST_Boundary(geom)) geom FROM land_parcels) lp

You can ST_Dump to break linestrings.

  • The query posted above returns polygon features. I want to create linestrings where 'land_use' boundaries do not math 'land_parcel' boundaries. Mike's query works properly after replacing ST_SymDifference with ST_Difference. – eclipsed_by_the_moon Jul 15 '16 at 8:54
  • It is really easy, just convert it by using ST_Boundary, I've edited the answer. – caiohamamura Jul 15 '16 at 15:25
  • Thanks, but the query does not return the results I need. Please take a closer look at the question. – eclipsed_by_the_moon Jul 16 '16 at 17:47
  • Ah, now I've got it, you want the "new" boundaries included by land_use which didn't exist before in land_parcels, right? I've updated above. – caiohamamura Jul 16 '16 at 18:18
  • Thanks, this one seems to return the expected results. But Mike's query is much faster on large data sets (e.g. 2.5 sec vs 109 sec for a few hundred features). – eclipsed_by_the_moon Jul 18 '16 at 13:58
0

Couldn't you just run an erase and end up with any stray line segments or add an erase to your script. Could you just covert your land use polygon to a polyline and run an intersect to see which areas actually overlap?

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.