1

Below, is a PostGIS SQL query where I have a table called nullss (polygon and line geometry) and a table called zones(polygon and line geometry as well) where the zones table spatially surrounds the nullss features.

PS** nullss does not mean the geometry is null, maybe its a bad name for a table but for me it makes sense

My Goal in this query is to find which zone type in the zones table shares the highest percentage of a nullss feature boundary, so that nullss feature will receive that zone type with the highest pct of feature that touches it.

In the query the first CTE is nullss, second CTE is zones,

third CTE calculates the sum length intersection of the nullss line with the line boundary/linestring zone

Problem

In some cases the intersection of the linestrings yields wrong calculations. Like the picture below

Picture Explanation AS you can see the NULLSS feature in this example is in pink and the majority of this feature is surrounded by the protection zone, and some by the conservation zone. when I run my query on this feature it yields wrong results

enter image description here

with nullss as(select objectid nid, shape as null_shape,st_boundary(shape) null_boundary,st_length((st_dump(st_boundary(shape))).geom) null_length,final_zone
            from null_test
            ),
zones as (select id zid, shape as zones_shape, st_boundary(shape) zone_boundary,final_zone
        from zones_test
        )

select nid, case 
        when sum(st_length(st_multi(st_collectionextract(st_intersection(null_boundary,st_snap(null_boundary,zone_boundary,0.001)),2))::geometry(MultiLineString,3424)))/null_length is null then 0
        else sum(st_length(st_multi(st_collectionextract(st_intersection(null_boundary,st_snap(null_boundary,zone_boundary,0.001)),2))::geometry(MultiLineString,3424)))/null_length
        end as pct,
            sum(st_length(st_multi(st_collectionextract(st_intersection(null_boundary,st_snap(null_boundary,zone_boundary,0.001)),2))::geometry(MultiLineString,3424))) intersect_length,
            st_union(st_multi(st_collectionextract(st_intersection(null_boundary,st_snap(null_boundary,zone_boundary,0.001)),2))::geometry(MultiLineString,3424)) line,
            null_length,zones.final_zone,nullss.null_shape 
            from nullss left join zones 
            on st_dwithin(null_boundary,zone_boundary,1)
            group by null_length,zones.final_zone,nullss.null_shape,nid

enter image description here

as you can see in the picture the conservation zone gets a higher intersect_length and pct, and it should be the protection zone with the higher numbers....

What I am struggling with

what is the optimal tolerance for this query? should I be intersecting linestrings and linestring or polygons and linestring?

GeoJSON from pic example

Here is the null_test table geojson

      {
"type": "FeatureCollection",
"name": "null_test",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3424" } },
"features": [
{ "type": "Feature", "properties": { "_uid_": 1.0, "nid": 4978, "null_shape": "SRID=3424;MULTIPOLYGON(((344454.802735861 687039.165515535,344412.472969595 687036.342903748,344418.018271845 687045.062163599,344464.241037186 687048.153686531,344465.205930281 687048.227505282,344472.698923375 687048.874546885,344476.812944379 687041.3", "null_lengt": 138.30200633145699, "final_zone": null }, "geometry": { "type": "LineString", "coordinates": [ [ 344454.802735861390829, 687039.165515534579754 ], [ 344412.472969595342875, 687036.34290374815464 ], [ 344418.018271844834089, 687045.062163598835468 ], [ 344464.241037186235189, 687048.153686530888081 ], [ 344465.205930281430483, 687048.227505281567574 ], [ 344472.698923375457525, 687048.874546885490417 ], [ 344476.812944378703833, 687041.303665973246098 ], [ 344455.570450853556395, 687039.228507541120052 ], [ 344454.802735861390829, 687039.165515534579754 ] ] } }
]
}

zones table

 {
"type": "FeatureCollection",
"name": "zones_test",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3424" } },
"features": [
{ "type": "Feature", "properties": { "_uid_": 1.0, "zid": 880, "zones_shap": "SRID=3424;POLYGON((345437.792011883 687350.612792753,345434.466492496 687343.090296954,345428.836168613 687335.076290473,345421.7516063 687328.313611709,345413.484639786 687323.061744273,345404.352472063 687319.52220197,345402.143786661 687318.906970404,", "final_zone": "Conservation" }, "geometry": { "type": "LineString", "coordinates": [ [ 344498.222791777690873, 687051.078868725686334 ], [ 344465.205737616866827, 687048.227750681340694 ], [ 344464.240802995860577, 687048.153826273977757 ], [ 344406.001681918860413, 687044.258611928322352 ] ] } },
{ "type": "Feature", "properties": { "_uid_": 2.0, "zid": 12544, "zones_shap": "SRID=3424;POLYGON((345651.125848219 688078.093516923,346707.275008645 686984.583528988,346675.620934352 686957.333075844,346668.770984292 686950.13285099,346663.480382003 686941.720087968,346659.95813819 686932.427141592,346658.343402643 686922.621139184", "final_zone": "Protection" }, "geometry": { "type": "LineString", "coordinates": [ [ 344400.592647661396768, 687035.550706182373688 ], [ 344454.802735831588507, 687039.165515519678593 ], [ 344455.570450838655233, 687039.22850751131773 ], [ 344502.59150800073985, 687043.821948615484871 ] ] } }
]
}

it seems PostGIS has a problem with topology, I cannot find a proper tolerance for this query. test out the query with above geojson

  • 1
    I don't suppose you can put up a couple of simple features (WKT) for each table to demonstrate what you are trying to achieve? – MickyT Dec 6 '17 at 18:32
  • thats a good idea, I will get on that – ziggy Dec 6 '17 at 18:36
  • @MickyT I provided the geojson for the pic example – ziggy Dec 6 '17 at 19:50
1

The issue is the order of your geometries in the ST_Snap. You are snapping the null_boundary to the zone, then intersecting the null boundary against it. This gives you the result of the whole null_boundary excluding the snapped portion of it.

Swapping them around and simplifying the query a little appears to give the right result.

WITH nullss AS (
    SELECT *
    FROM (VALUES
         (1.0, 4978, 'SRID=3424;MULTIPOLYGON(((344454.802735861 687039.165515535,344412.472969595 687036.342903748,344418.018271845 687045.062163599,344464.241037186 687048.153686531,344465.205930281 687048.227505282,344472.698923375 687048.874546885,344476.812944379 687041.3)))', 138.30200633145699, null, 'SRID=3424;LineString(344454.802735861390829 687039.165515534579754 , 344412.472969595342875 687036.34290374815464 , 344418.018271844834089 687045.062163598835468 , 344464.241037186235189 687048.153686530888081 , 344465.205930281430483 687048.227505281567574 , 344472.698923375457525 687048.874546885490417 , 344476.812944378703833 687041.303665973246098 , 344455.570450853556395 687039.228507541120052 , 344454.802735861390829 687039.165515534579754)'::Geometry)
        ) nulls(uid,nid,null_shape,null_length,final_zone,null_boundary)
)
,zones AS (
    SELECT *
    FROM (VALUES
         (1.0, 880,   'SRID=3424;POLYGON((345437.792011883 687350.612792753,345434.466492496 687343.090296954,345428.836168613 687335.076290473,345421.7516063 687328.313611709,345413.484639786 687323.061744273,345404.352472063 687319.52220197,345402.143786661 687318.906970404,','Conservation','SRID=3424;LineString(344498.222791777690873  687051.078868725686334 , 344465.205737616866827  687048.227750681340694 , 344464.240802995860577  687048.153826273977757 , 344406.001681918860413 687044.258611928322352)'::Geometry)
        ,(2.0, 12544, 'SRID=3424;POLYGON((345651.125848219 688078.093516923,346707.275008645 686984.583528988,346675.620934352 686957.333075844,346668.770984292 686950.13285099,346663.480382003 686941.720087968,346659.95813819 686932.427141592,346658.343402643 686922.621139184','Protection',  'SRID=3424;LineString(344400.592647661396768  687035.550706182373688 , 344454.802735831588507  687039.165515519678593 , 344455.570450838655233  687039.22850751131773  , 344502.59150800073985  687043.821948615484871)'::Geometry)
        ) zones(uid,zid,zones_shap, final_zone, zone_boundary)
)
select nid
     ,coalesce(st_length(st_intersection(null_boundary, st_snap(zone_boundary,null_boundary,0.001))) / st_length(null_boundary),0) pct
     ,st_length(st_intersection(null_boundary, st_snap(zone_boundary,null_boundary,0.001))) intersect_length
     ,st_astext(st_intersection(null_boundary, st_snap(zone_boundary,null_boundary,0.001))) line
     ,null_length
     ,zones.final_zone
     ,nullss.null_shape
from nullss 
    left join zones on st_dwithin(null_boundary,zone_boundary,0.001)

With feedback from @ziggy, a ST_Snap tolerance of 0.002 worked best. Remember to adjust the st_within tolerance as well

  • okay this worked for the above example, but why would switching the inputs on the st_snap accomplish this? – ziggy Dec 7 '17 at 19:15
  • running this now for all the features, ps* needed the case statement for because some features do not intersect and the nulls throw off a part of the query I did not include in the question – ziggy Dec 7 '17 at 19:15
  • oh I guess I could use the coalesce function instead... – ziggy Dec 7 '17 at 19:17
  • @ziggy I added a little explanation of what the problem was. Essentially you were intersecting the null_boundary with itself – MickyT Dec 7 '17 at 19:29
  • thanks for the answer, FYI I used a st_snap tolerance of 0.002 , that yielded the most accurate results for me – ziggy Dec 8 '17 at 13:54

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.