Thursday 6 February 2020

postgis - Find all intersections of a LineString and a Polygon and the order in which it intersects it



I would like to find out all intersections of a LineString (River) through Polygons (Municipalities). I would like to list the order of municipalities which the river crosses, using linear referencing. I already have found which municipalities are crossed by the river using st_intersects, but this is will not provide me the order and how many times municipalities are crossed. I am not expecting a full solution, rather a modus operandi that I will implement by my self.


I use PostGIS 2 extension with PostgreSQL 9.3


MWE:


WITH 
A AS (
SELECT
curves.gid AS cgid,
shapes.gid AS sgid,
curves.geom AS cgeom,
shapes.geom AS sgeom,

st_intersection(curves.geom, shapes.geom) AS mgeom
FROM
curves, shapes),
B AS (
SELECT
cgid, sgid, cgeom, sgeom,
(st_dump(mgeom)).geom AS geom
FROM A),
C AS (
SELECT

row_number() OVER ()::INTEGER as gid,
cgid, sgid,
cgeom, sgeom, geom,
st_startPoint(geom) AS p0,
st_endPoint(geom) AS p1
FROM B),
D AS (
SELECT
C.*,
st_Line_Locate_Point(cgeom, p0) AS r0,

st_Line_Locate_Point(cgeom, p1) AS r1
FROM C
ORDER BY r0)

SELECT * FROM D

The query above find entry points and may be used in Linear Referencing in order to follow the original LineString through municipalities. I have to discriminate many use case and I am a little lost. This is my very first queries in GIS. Table shapes contains Polygons that have no intersections and curves contains LineString that pass through those Polygons. The expected output is, for each LineString, the list of Polygons that it crosses (order and redondance matters).



Answer




but this is will not provide me the order and how many times municipalities are crossed.




Assuming every row in Rivers table is a linestring with an entire river. Here's a query that will get you how many municipalities are crossed, the subquery t , will get you all the municipalities that each river crosses but yes they will not be in any guaranteed order.


SELECT t.river_gid as river_gid, count(*) as n_mun FROM (
SELECT rivers.gid as river_gid, municipalities.gid as mun_gid
FROM rivers, municipalities
WHERE ST_Crosses(rivers.geom, municipalities.geom)
) as t group by t.river_gid;

here's a query that will get you all the points from the linestring that intersected the polygon in the correct order.


SELECT municipalities.id as mun_id, dp.path, dp.river_gid FROM municipalities, (

SELECT (ST_DumpPoints(geom)).path[1] as path,
(ST_DumpPoints(geom)).geom as geom,
rivers.id as river_gid from rivers
) as dp
WHERE ST_INTERSECTS(dp.geom, municipalities.geom)
ORDER BY dp.river_gid, dp.path;

Unfortunately there is no way to know how many times a river crossed a municipality without some code, but with this query well get you very close. As result you will get something like this:


mun_id | path | river_gid
1 1 1

1 2 1
2 3 1
2 4 1
1 5 1
1 6 1

In this example the river crossed 2 times the same municipality, and the way to know it is the mun_id column pattern (It was on mun_id 1 and then changed to 2, and then back to 1).


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