[postgis-users] Distance from Rast and LineString overlay

Fri Oct 27 01:08:43 PDT 2023

Hi Nicolas,

Thanks for your reply. Your solution worked perfectly and faster than anything I was trying, overlaying with Polygons or Raster. Using Linestrings together with Subdivide function made it much faster. Also, thanks for your correction on how to build the line in the right direction.

Cheers!

Eloi

------- Original Message -------
On Wednesday, October 25th, 2023 at 14:51, Nicolas Ribot via postgis-users <postgis-users at lists.osgeo.org> wrote:

> Here are some tests made with the data you pointed out:
>
> gdal_rasterize was not that fast for this big dataset, but very reasonable IMOO (46min):
>
> time
>
> gdal_polygonize.py Forest_2020_Europe.tif -f PostgreSQL PG:"port=16432" forest val
>
> \
>
> -lco FID=forest_id -lco GEOMETRY_NAME=geom
>
> #
>
> 46
>
> :
>
> 19.61
>
> total
>
> It generated 5 343 630 polygons in database.
>
> Then some sql to create lines, intersect with forest boundaries and compute distances:
>
> -- gdal inserted a new srid instead of using 3035 (tif srid as I understood), there are equivalent
>
> -- extract forest polygon boundaries and subdivide them for performance
>
> drop table if exists
>
> forest_bound
>
> ;
>
> create table
>
> forest_bound
>
> as
>
> select
>
> ogc_fid
>
> as
>
> forest_id
>
> ,
>
> st_setSRID
>
> (
>
> st_subdivide
>
> (
>
> st_boundary
>
> (
>
> wkb_geometry
>
> ))
>
> ,
>
> 3035
>
> )::
>
> geometry
>
> (linestring
>
> ,
>
> 3035
>
> )
>
> as
>
> geom
>
> from
>
> forest
>
> ;
>
> -- [2023-10-25 12:38:36] 10,153,846 rows affected in 28 s 259 ms
>
> alter table
>
> forest_bound
>
> add column
>
> forest_bound_id
>
> int generated always as identity primary key;
>
> create index on
>
> forest_bound
>
> using
>
> gist (
>
> geom
>
> )
>
> ;
>
> -- [2023-10-25 12:40:03] completed in 10 s 21 ms
>
> vacuum analyze
>
> forest_bound
>
> ;
>
> -- random points created in QGIS inside forest polygons, stored into a table named "point"
>
> -- lines from points (st_project can be used for that, though I saw strange behaviour when using data's srid)
>
> drop table if exists
>
> line
>
> ;
>
> create table
>
> line
>
> as
>
> select
>
> id
>
> as
>
> id_point
>
> ,
>
> st_transform
>
> (
>
> st_makeLine
>
> (
>
> st_transform
>
> (
>
> geom
>
> ,
>
> 4326
>
> )::
>
> geometry,
>
> st_project
>
> (
>
> st_transform
>
> (
>
> geom
>
> ,
>
> 4326
>
> )::
>
> geography,
>
> 1000
>
> ,
>
> radians
>
> (
>
> 300.0
>
> ))::
>
> geometry
>
> )
>
> ,
>
> 3035
>
> )::
>
> geometry
>
> (linestring
>
> ,
>
> 3035
>
> )
>
> as
>
> geom
>
> --        st_makeLine(geom, st_project(geom, 1000, radians(300.0)))::geometry(linestring, 3035) as geom -- does not work
>
> from
>
> point
>
> ;
>
> create index on
>
> line
>
> using
>
> gist (
>
> geom
>
> )
>
> ;
>
> vacuum analyze
>
> line
>
> ;
>
> -- intersection between lines and forest boundaries: will produce n points
>
> -- (or linestring in some rare cases, can be filtered out with st_collectionExtract(geomcol, 1))
>
> -- lead() windows function to get original point and the next intersection points between line and forest boundaries, ordered by distance
>
> -- we also add (with union all) the original point to the collection to have all geometries we need to compute distances
>
> -- finally compute distance between consecutive points, starting with original points
>
> with
>
> tmp
>
> as
>
> (
>
> select
>
> l.
>
> id_point
>
> ,
>
> (
>
> st_dump
>
> (
>
> st_intersection
>
> (l.
>
> geom
>
> ,
>
> fb.
>
> geom
>
> ))).geom
>
> as
>
> geom
>
> from
>
> line l
>
> join
>
> forest_bound fb
>
> on
>
> st_intersects
>
> (l.
>
> geom
>
> ,
>
> fb.
>
> geom
>
> )
>
> join
>
> point pt
>
> on
>
> l.
>
> id_point
>
> = pt.
>
> id
>
> union all
>
> select
>
> id
>
> ,
>
> geom
>
> from
>
> point
>
> )
>
> ,
>
> tmp1
>
> as
>
> (
>
> select
>
> id_point
>
> ,
>
> row_number
>
> ()
>
> over w as
>
> dist_num
>
> ,
>
> tmp.geom
>
> ,
>
> lead
>
> (tmp.geom)
>
> over w as
>
> next_geom
>
> from
>
> tmp
>
> join
>
> point p
>
> on
>
> tmp.
>
> id_point
>
> = p.
>
> id
>
> window w as
>
> (
>
> partition by
>
> id_point
>
> order by
>
> st_distance
>
> (p.
>
> geom
>
> ,
>
> tmp.geom))
>
> )
>
> select
>
> id_point
>
> ,
>
> format
>
> (
>
> 'distance_%s'
>
> ,
>
> chr
>
> (dist_num::
>
> int
>
> +
>
> 64
>
> ))
>
> as
>
> dist_num
>
> ,
>
> st_distance
>
> (geom
>
> ,
>
> next_geom)
>
> as
>
> dist
>
> from
>
> tmp1
>
> where
>
> next_geom
>
> is not null;
>
> -- only lines intersecting with forest boundaries are returned, lines fully inside a forest are filtered out here
>
> -- +--------+----------+------------------+
>
> -- |id_point|dist_num  |dist              |
>
> -- +--------+----------+------------------+
>
> -- |1       |distance_A|743.6784090392044 |
>
> -- |3       |distance_A|746.1313796860095 |
>
> -- |3       |distance_B|129.52215146305574|
>
> -- |4       |distance_A|74.89447054570846 |
>
> -- |4       |distance_B|500.32736839280386|
>
> -- |9       |distance_A|262.47720947329753|
>
> -- |9       |distance_B|113.28881636877391|
>
> -- |10      |distance_A|54.10140760654185 |
>
> -- |10      |distance_B|113.35161581910967|
>
> -- |10      |distance_C|221.62871568295094|
>
> -- |10      |distance_D|118.42613177437802|
>
> -- |11      |distance_A|182.13323464032416|
>
> -- |12      |distance_A|389.19304731092586|
>
> -- |12      |distance_B|472.5385803172015 |
>
> -- |13      |distance_A|506.77444169025375|
>
> -- |13      |distance_B|80.19996080704423 |
>
> -- |14      |distance_A|61.194892187708376|
>
> -- |14      |distance_B|454.7118937668803 |
>
> -- |14      |distance_C|1.6389730162778315|
>
> -- |14      |distance_D|342.2631500871444 |
>
> -- |14      |distance_E|114.08771669601367|
>
> -- ...
>
> -- +--------+----------+------------------+
>
> A screenshot of some data:
> [Screenshot 2023-10-25 at 14.48.22.png]
>
> I'm sure there are some edge cases to deal with, but on the small point dataset I tested it seems to work.
>
> Nicolas
>
> On Wed, 25 Oct 2023 at 11:28, Nicolas Ribot <nicky666 at gmail.com> wrote:
>
>> Hello,
>>
>> I would polygonize the raster with gdal_polygonize in the database, as it is very fast. Then extract the boundaries (linestring), then subdivide them for better performance.
>> From there, I would intersect the 1km lines for each point with forest boundaries and sort intersection points by distance. First point would give distanceA, second point distance B, and so on.
>> By using polygon boundaries, I don't think you need st_touches or other predicates: intersection between boundaries and direction line gives you the result.
>>
>> Nicolas
>>
>> On Tue, 24 Oct 2023 at 18:27, Eloi Ribeiro via postgis-users <postgis-users at lists.osgeo.org> wrote:
>>
>>> Dear PostGIS users,
>>>
>>> I have a [raster](https://data.jrc.ec.europa.eu/dataset/35fb1231-849b-4017-89d8-9cadbaf9d555) (of 100m pixel) of forest (val=1) / no forest (val=NULL) and a set of points that always fall inside forest areas.
>>>
>>> I need to calculate the distance (up to 1km), in a given direction (let's say West, 300 radians), to the border of the forest patch where it falls in (distance A in fig.). Then, the distance from the previous forest border to the next forest border (distance B in fig.), following the same direction. Like so:
>>>
>>> [Screenshot from 2023-10-24 18-19-16.png]
>>>
>>> I was able to do this with the following steps:
>>> 1. Vectorize forest map. Note the ST_Union, to eliminate the borders of the raster tiles and avoid these being confused as forest border.
>>>
>>> SELECT (ST_Dump(ST_Union(geom))).geom::geometry(Polygon,3035)
>>> FROM (
>>> SELECT dp.*
>>> FROM storm.forest_2020_europe,
>>> LATERAL ST_DumpAsPolygons(rast) AS dp
>>> ) As foo
>>> GROUP BY valHAVING val = 1;
>>>
>>> 2. Create Linestring from point, as starting point, following 300 radians and 1km length.
>>>
>>> SELECT gid,
>>> '1km 300 rad',
>>> ST_MakeLine(geom, ST_Translate(geom, sin(300)*1000, cos(300)*1000)) AS geomFROM storm.point;
>>>
>>> 3. Overlay Linestring from step 2 with vectorize forest map and get the distances from Linestring segments. This part has too much SQL to put it here. I have used ST_Intersection and ST_Difference, and then ST_Touches to make sure the segment touches the first point, for the first distance and touches the EndPoint of the first segment to get the second distance.
>>>
>>> This at national level worked just fine. But now at continent level, the resulting vectorized polygons of the forest patches are just too big and takes forever to vectorize. So, I need to do this directly with the raster forest map.
>>>
>>> Any hints how could this be accomplished?
>>> Thanks!
>>>
>>> Cheers,
>>>
>>> Eloi
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/postgis-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20231027/8f1d01aa/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot from 2023-10-24 18-19-16.png
Type: image/png
Size: 55040 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20231027/8f1d01aa/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot 2023-10-25 at 14.48.22.png
Type: image/png
Size: 263279 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20231027/8f1d01aa/attachment-0001.png>