[GIS] Counting self intersections of linestring using PostGIS

intersectionlinestringpostgispostgresql

Is there a way of using PostgreSQL + PostGIS for finding the number of self intersections in a linestring?

I have a linestring that goes forward but in circles, and I need to find the number of times it crosses itself.

Best Answer

I don't know for sure, but I think something like this might solve your problem (this wont work on multilines, but you could stdump them or squish them) :

SELECT  ID, 
        COUNT(*) AS theCOunt
FROM
(
  SELECT  id, 
          ST_Intersection(A.geom, B.geom) AS geom
  FROM
  (
    SELECT  xx.id,
            xx.id2,
            st_makeline(xx.l_start, xx.l_end) AS geom
     FROM 
     ( 
       SELECT  x.id,
               generate_series(1, st_npoints(x.geom) - 1) AS id2,
               st_pointn(x.geom, generate_series(1, st_npoints(x.geom) - 1)) AS l_start,
               st_pointn(x.geom, generate_series(2, st_npoints(x.geom))) AS l_end
       FROM 
       ( 
         SELECT  DATA_TABLE.id,
                 DATA_TABLE.geom 
         FROM DATA_TABLE
        ) AS X
     ) AS XX
 ) AS A
 INNER JOIN
 (
    SELECT  xx.id,
            xx.id2,
            st_makeline(xx.l_start, xx.l_end) AS geom
     FROM 
     ( 
       SELECT  x.id,
               generate_series(1, st_npoints(x.geom) - 1) AS id2,
               st_pointn(x.geom, generate_series(1, st_npoints(x.geom) - 1)) AS l_start,
               st_pointn(x.geom, generate_series(2, st_npoints(x.geom))) AS l_end
       FROM 
       ( 
         SELECT  DATA_TABLE.id,
                 DATA_TABLE.geom 
         FROM DATA_TABLE
        ) AS X
     ) AS XX
 ) AS B
ON A.ID = b.ID
   AND
   A.id2 <> B.id2
   AND
   A.id2 <> B.id2-1
   AND
   A.id2 < B.id2
   AND
   ST_INTERSECTS(A.geom, B.geom)

Related Solutions

PostGIS – Finding and Ordering Intersections of LineString and Polygon

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

[GIS] Strange behavior of of ST_Intersects in PostGIS

The question has indeed been asked before, but it is still a good question :) The manual already gives you a hint towards the problem.

The main problem is that you are trying to use geograpy types (that work with real world coordinates) while the function st_intersection is not completely ready for geography.

In the real (spherical) world, your geometries would overlap, since the top of the big square you give would follow the circle of the earth (and show as convex arc on your 2D map). ST_Intersects knows this for geography and hence gives 'true' ST_Intersection on the other hand is less smart and uses a 2D projection to calculate where your square will remain a square. It chooses the 2D projection itself but doesn't always choose right and defaults to mercator (especially not for larger areas). There is a nice explanation overhere

Your solution will be to first transform the geography to some appropriate geometry (don't know by heart) and then do the intersection.

Best Answer

Related Solutions

PostGIS – Finding and Ordering Intersections of LineString and Polygon

[GIS] Strange behavior of of ST_Intersects in PostGIS

Related Question