Search code examples
postgisopenstreetmap

Retrieve length of "linestring" sub-segment within given area (using ST_DWithin)


I need to retrieve the length of the piece of a linestring geometry object inside a circle (point + radius).

I'm currently using ST_DWithin function but unfortunately it returns the full object, not the piece of object inside the area. Indeed, changing the radius (100 to 70) does not change the retrieved length.

osm_data=# select highway, st_length(way) from planet_osm_line where ST_DWITHIN(ST_Transform(way,4326),ST_GeomFromText('POINT(4.1884918 51.144580)',4326)::geography,100) and highway is not null;
   highway    |    st_length     
--------------+------------------
 motorway     | 5079.24632083105
 unclassified | 1110.56834915499
 motorway     | 1499.83459080537
(3 rows)

osm_data=# select highway, st_length(way) from planet_osm_line where ST_DWITHIN(ST_Transform(way,4326),ST_GeomFromText('POINT(4.1884918 51.144580)',4326)::geography,70) and highway is not null;
   highway    |    st_length     
--------------+------------------
 motorway     | 5079.24632083105
 unclassified | 1110.56834915499
 motorway     | 1499.83459080537
(3 rows)

Solution

  • ST_DWithin is used to select geometries that are at least partially within the said distance from the reference point.

    To compute the actual length (or area for polygons) within, you would need to compute the intersection.

    select highway, 
       st_length(way) full_length,   
       st_length(
         ST_INTERSECTION(
            ST_Transform(way,4326)::geography,
            ST_BUFFER(
              ST_GeomFromText('POINT(4.1884918 51.144580)',4326)::geography,
             100)
          )
       )
      from planet_osm_line 
      where
        ST_DWITHIN(
          ST_Transform(way,4326)::geography,
          ST_GeomFromText('POINT(4.1884918 51.144580)',4326)::geography,100) 
        and highway is not null;
    

    To optimize this query, you would use

    WITH buff as (
       SELECT ST_BUFFER(
                ST_GeomFromText('POINT(4.1884918 51.144580)',4326)::geography, 
                1000) as geog)  
    select highway, 
       st_length(way) full_length,   
       st_length(
         ST_INTERSECTION(
            ST_Transform(way,4326)::geography,
            buff.geog
          )
       )
    from osm.osm_line
        INNER JOIN buff
            ON ST_DWITHIN(
              ST_Transform(way,4326),
              buff.geog,
              1000) 
    where highway is not null;
    

    Please note that a buffer is only an approximation of a circle, so the result may be slightly off