Tuesday 30 January 2018

postgresql - postgis st_intersection of linestrings geometric miscalculations


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



Answer



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...