Saturday, 7 July 2018

intersection - Recover M-values from PostGIS `ST_Intersection`?


I am querying the intersection of a path (LinestringZM) through a region (MultiPolygon) but the ST_Intersection returns only the Lat/Long/Alt coordinates and I need to recover time, stored in the M values.


SELECT path.name AS path_name,
ST_Length(path.geom::geography) AS distance,
ST_Length(geography(ST_Intersection(path.geom,region.geom))) AS region_distance,

ST_M(ST_EndPoint(path.geom)) - ST_M(ST_StartPoint(path.geom)) AS time,
ST_M(ST_EndPoint(ST_Intersection(path.geom,region.geom))) - ST_M(ST_StartPoint(ST_Intersection(path.geom,region.geom))) AS region_time
region.name AS region_name
FROM path INNER JOIN region ON ST_Intersects(path.geom, region.geom)
ORDER BY vehicle_name, region_distance DESC;

How can I get the exact elapsed time within a region? (the above region_time is null because M coordinate is dropped. I was thinking of using LATERAL select to get the equivalent matching segment of the original path like so: ST_SharedPaths(path.geom,rel.intersection_path) ... LATERAL (SELECT ST_Intersection(path.geom,region.geom) AS intersection_path) AS rel


But that gives me ERROR: GEOSSharedPaths: TopologyException: found non-noded intersection between LINESTRING (...


I'm using PostGIS 2.4 and PostgreSQL 10.1 beta, in case it matters.



Answer




I use several functions to do the deed which you can grab here https://github.com/pauldzy/pg_dz_lrs/blob/master/Functions/LRS_INTERSECTION.sql


We could log an enhancement but I believe all the tepid M support derives from the underlying GEOS implementation.


Cheers, Paul


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...