<div dir="ltr">Here are some tests made with the data you pointed out:<div><br></div><div>gdal_rasterize was not that fast for this big dataset, but very reasonable IMOO (46min):</div><div><br></div><div><div style="background-color:rgb(43,43,43);color:rgb(169,183,198)"><pre style="font-family:Consolas,monospace;font-size:10.5pt"><span style="color:rgb(204,120,50)">time </span>gdal_polygonize.py Forest_2020_Europe.tif -f PostgreSQL PG:"port=16432" forest val <span style="background-color:rgb(54,65,53)">\</span><br> -lco FID=forest_id -lco GEOMETRY_NAME=geom<br># <span style="color:rgb(104,151,187)">46</span>:<span style="color:rgb(104,151,187)">19.61 </span>total <br></pre></div></div><div>It generated 5 343 630 polygons in database.</div><div><br></div><div>Then some sql to create lines, intersect with forest boundaries and compute distances:</div><div><br></div><div><div style="background-color:rgb(43,43,43);color:rgb(169,183,198)"><pre style="font-family:Consolas,monospace;font-size:10.5pt"><span style="color:rgb(128,128,128)">-- gdal inserted a new srid instead of using 3035 (tif srid as I understood), there are equivalent<br></span><span style="color:rgb(128,128,128)">-- extract forest polygon boundaries and subdivide them for performance<br></span><span style="color:rgb(204,120,50)">drop table if exists </span>forest_bound<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)">create table </span>forest_bound <span style="color:rgb(204,120,50)">as<br></span><span style="color:rgb(204,120,50)"> select </span><span style="color:rgb(152,118,170)">ogc_fid </span><span style="color:rgb(204,120,50)">as </span>forest_id<span style="color:rgb(204,120,50)">, </span><span style="color:rgb(255,198,109);font-style:italic">st_setSRID</span>(<span style="color:rgb(255,198,109);font-style:italic">st_subdivide</span>(<span style="color:rgb(255,198,109);font-style:italic">st_boundary</span>(<span style="color:rgb(152,118,170)">wkb_geometry</span>))<span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">3035</span>)::<span style="color:rgb(204,120,50)">geometry</span>(linestring<span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">3035</span>) <span style="color:rgb(204,120,50)">as </span>geom<br> <span style="color:rgb(204,120,50)">from </span>forest<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(128,128,128)">-- [2023-10-25 12:38:36] 10,153,846 rows affected in 28 s 259 ms<br></span><span style="color:rgb(128,128,128)"><br></span><span style="color:rgb(204,120,50)">alter table </span>forest_bound<br> <span style="color:rgb(204,120,50)">add column </span><span style="color:rgb(152,118,170)">forest_bound_id </span><span style="color:rgb(204,120,50)">int generated always as identity primary key;<br></span><span style="color:rgb(204,120,50)">create index on </span>forest_bound <span style="color:rgb(204,120,50)">using </span>gist (<span style="color:rgb(152,118,170)">geom</span>)<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(128,128,128)">-- [2023-10-25 12:40:03] completed in 10 s 21 ms<br></span><span style="color:rgb(204,120,50)">vacuum analyze </span>forest_bound<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)"><br></span><span style="color:rgb(128,128,128)">-- random points created in QGIS inside forest polygons, stored into a table named "point"<br></span><span style="color:rgb(128,128,128)"><br></span><span style="color:rgb(128,128,128)">-- lines from points (st_project can be used for that, though I saw strange behaviour when using data's srid)<br></span><span style="color:rgb(204,120,50)">drop table if exists </span>line<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)">create table </span>line <span style="color:rgb(204,120,50)">as<br></span><span style="color:rgb(204,120,50)">select </span><span style="color:rgb(152,118,170)">id </span><span style="color:rgb(204,120,50)">as </span>id_point<span style="color:rgb(204,120,50)">,<br></span><span style="color:rgb(204,120,50)"> </span><span style="color:rgb(255,198,109);font-style:italic">st_transform</span>(<span style="color:rgb(255,198,109);font-style:italic">st_makeLine</span>(<br> <span style="color:rgb(255,198,109);font-style:italic">st_transform</span>(<span style="color:rgb(152,118,170)">geom</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">4326</span>)::<span style="color:rgb(204,120,50)">geometry,<br></span><span style="color:rgb(204,120,50)"> </span><span style="color:rgb(255,198,109);font-style:italic">st_project</span>(<span style="color:rgb(255,198,109);font-style:italic">st_transform</span>(<span style="color:rgb(152,118,170)">geom</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">4326</span>)::<span style="color:rgb(204,120,50)">geography, </span><span style="color:rgb(104,151,187)">1000</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(255,198,109);font-style:italic">radians</span>(<span style="color:rgb(104,151,187)">300.0</span>))::<span style="color:rgb(204,120,50)">geometry<br></span><span style="color:rgb(204,120,50)"> </span>)<span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">3035</span>)::<span style="color:rgb(204,120,50)">geometry</span>(linestring<span style="color:rgb(204,120,50)">, </span><span style="color:rgb(104,151,187)">3035</span>) <span style="color:rgb(204,120,50)">as </span>geom<br><span style="color:rgb(128,128,128)">-- st_makeLine(geom, st_project(geom, 1000, radians(300.0)))::geometry(linestring, 3035) as geom -- does not work<br></span><span style="color:rgb(204,120,50)">from </span>point<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(168,192,35);font-style:italic"><br></span><span style="color:rgb(204,120,50)">create index on </span>line <span style="color:rgb(204,120,50)">using </span>gist (<span style="color:rgb(152,118,170)">geom</span>)<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)">vacuum analyze </span>line<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)"><br></span><span style="color:rgb(128,128,128)">-- intersection between lines and forest boundaries: will produce n points<br></span><span style="color:rgb(128,128,128)">-- (or linestring in some rare cases, can be filtered out with st_collectionExtract(geomcol, 1))<br></span><span style="color:rgb(128,128,128)">-- lead() windows function to get original point and the next intersection points between line and forest boundaries, ordered by distance<br></span><span style="color:rgb(128,128,128)">-- we also add (with union all) the original point to the collection to have all geometries we need to compute distances<br></span><span style="color:rgb(128,128,128)">-- finally compute distance between consecutive points, starting with original points<br></span><span style="color:rgb(204,120,50)">with </span>tmp <span style="color:rgb(204,120,50)">as </span>(<br> <span style="color:rgb(204,120,50)">select </span>l.<span style="color:rgb(152,118,170)">id_point</span><span style="color:rgb(204,120,50)">, </span>(<span style="color:rgb(255,198,109);font-style:italic">st_dump</span>(<span style="color:rgb(255,198,109);font-style:italic">st_intersection</span>(l.<span style="color:rgb(152,118,170)">geom</span><span style="color:rgb(204,120,50)">, </span>fb.<span style="color:rgb(152,118,170)">geom</span>))).geom <span style="color:rgb(204,120,50)">as </span>geom<br> <span style="color:rgb(204,120,50)">from </span>line l<br> <span style="color:rgb(204,120,50)">join </span>forest_bound fb <span style="color:rgb(204,120,50)">on </span><span style="color:rgb(255,198,109);font-style:italic">st_intersects</span>(l.<span style="color:rgb(152,118,170)">geom</span><span style="color:rgb(204,120,50)">, </span>fb.<span style="color:rgb(152,118,170)">geom</span>)<br> <span style="color:rgb(204,120,50)">join </span>point pt <span style="color:rgb(204,120,50)">on </span>l.<span style="color:rgb(152,118,170)">id_point </span>= pt.<span style="color:rgb(152,118,170)">id<br></span><span style="color:rgb(152,118,170)"> </span><span style="color:rgb(204,120,50)">union all<br></span><span style="color:rgb(204,120,50)"> select </span><span style="color:rgb(152,118,170)">id</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(152,118,170)">geom<br></span><span style="color:rgb(152,118,170)"> </span><span style="color:rgb(204,120,50)">from </span>point<br>)<span style="color:rgb(204,120,50)">, </span>tmp1 <span style="color:rgb(204,120,50)">as </span>(<br> <span style="color:rgb(204,120,50)">select </span><span style="color:rgb(152,118,170)">id_point</span><span style="color:rgb(204,120,50)">,<br></span><span style="color:rgb(204,120,50)"> </span><span style="color:rgb(255,198,109);font-style:italic">row_number</span>() <span style="color:rgb(204,120,50)">over w as </span>dist_num<span style="color:rgb(204,120,50)">,<br></span><span style="color:rgb(204,120,50)"> </span>tmp.geom<span style="color:rgb(204,120,50)">,<br></span><span style="color:rgb(204,120,50)"> </span><span style="color:rgb(255,198,109);font-style:italic">lead</span>(tmp.geom) <span style="color:rgb(204,120,50)">over w as </span>next_geom<br> <span style="color:rgb(204,120,50)">from </span>tmp <span style="color:rgb(204,120,50)">join </span>point p <span style="color:rgb(204,120,50)">on </span>tmp.<span style="color:rgb(152,118,170)">id_point </span>= p.<span style="color:rgb(152,118,170)">id<br></span><span style="color:rgb(152,118,170)"> </span><span style="color:rgb(204,120,50)">window w as </span>(<span style="color:rgb(204,120,50)">partition by </span><span style="color:rgb(152,118,170)">id_point </span><span style="color:rgb(204,120,50)">order by </span><span style="color:rgb(255,198,109);font-style:italic">st_distance</span>(p.<span style="color:rgb(152,118,170)">geom</span><span style="color:rgb(204,120,50)">, </span>tmp.geom))<br>)<br><span style="color:rgb(204,120,50)">select </span><span style="color:rgb(152,118,170)">id_point</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(255,198,109);font-style:italic">format</span>(<span style="color:rgb(106,135,89)">'distance_%s'</span><span style="color:rgb(204,120,50)">, </span><span style="color:rgb(255,198,109);font-style:italic">chr</span>(dist_num::<span style="color:rgb(204,120,50)">int </span>+ <span style="color:rgb(104,151,187)">64</span>)) <span style="color:rgb(204,120,50)">as </span>dist_num<span style="color:rgb(204,120,50)">,<br></span><span style="color:rgb(204,120,50)"> </span><span style="color:rgb(255,198,109);font-style:italic">st_distance</span>(geom<span style="color:rgb(204,120,50)">, </span>next_geom) <span style="color:rgb(204,120,50)">as </span>dist<br><span style="color:rgb(204,120,50)">from </span>tmp1<br><span style="color:rgb(204,120,50)">where </span>next_geom <span style="color:rgb(204,120,50)">is not null;</span></pre><pre style="font-family:Consolas,monospace;font-size:10.5pt">-- only lines intersecting with forest boundaries are returned, lines fully inside a forest are filtered out here<span style="color:rgb(204,120,50)"><br></span><span style="color:rgb(128,128,128)">-- +--------+----------+------------------+<br></span><span style="color:rgb(128,128,128)">-- |id_point|dist_num |dist |<br></span><span style="color:rgb(128,128,128)">-- +--------+----------+------------------+<br></span><span style="color:rgb(128,128,128)">-- |1 |distance_A|743.6784090392044 |<br></span><span style="color:rgb(128,128,128)">-- |3 |distance_A|746.1313796860095 |<br></span><span style="color:rgb(128,128,128)">-- |3 |distance_B|129.52215146305574|<br></span><span style="color:rgb(128,128,128)">-- |4 |distance_A|74.89447054570846 |<br></span><span style="color:rgb(128,128,128)">-- |4 |distance_B|500.32736839280386|<br></span><span style="color:rgb(128,128,128)">-- |9 |distance_A|262.47720947329753|<br></span><span style="color:rgb(128,128,128)">-- |9 |distance_B|113.28881636877391|<br></span><span style="color:rgb(128,128,128)">-- |10 |distance_A|54.10140760654185 |<br></span><span style="color:rgb(128,128,128)">-- |10 |distance_B|113.35161581910967|<br></span><span style="color:rgb(128,128,128)">-- |10 |distance_C|221.62871568295094|<br></span><span style="color:rgb(128,128,128)">-- |10 |distance_D|118.42613177437802|<br></span><span style="color:rgb(128,128,128)">-- |11 |distance_A|182.13323464032416|<br></span><span style="color:rgb(128,128,128)">-- |12 |distance_A|389.19304731092586|<br></span><span style="color:rgb(128,128,128)">-- |12 |distance_B|472.5385803172015 |<br></span><span style="color:rgb(128,128,128)">-- |13 |distance_A|506.77444169025375|<br></span><span style="color:rgb(128,128,128)">-- |13 |distance_B|80.19996080704423 |<br></span><span style="color:rgb(128,128,128)">-- |14 |distance_A|61.194892187708376|<br></span><span style="color:rgb(128,128,128)">-- |14 |distance_B|454.7118937668803 |<br></span><span style="color:rgb(128,128,128)">-- |14 |distance_C|1.6389730162778315|<br></span><span style="color:rgb(128,128,128)">-- |14 |distance_D|342.2631500871444 |<br></span><span style="color:rgb(128,128,128)">-- |14 |distance_E|114.08771669601367|<br></span><span style="color:rgb(128,128,128)">-- ...<br></span><span style="color:rgb(128,128,128)">-- +--------+----------+------------------+<br></span></pre></div></div><div><br></div><div>A screenshot of some data:</div><div><br></div><img src="cid:ii_lo5r4jei1" alt="Screenshot 2023-10-25 at 14.48.22.png" width="563" height="477"><br><div><br></div><div>I'm sure there are some edge cases to deal with, but on the small point dataset I tested it seems to work.</div><div><br></div><div>Nicolas</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, 25 Oct 2023 at 11:28, Nicolas Ribot <<a href="mailto:nicky666@gmail.com">nicky666@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hello,<div><br></div><div>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.</div><div>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.</div><div>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.</div><div><br></div><div>Nicolas</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, 24 Oct 2023 at 18:27, Eloi Ribeiro via postgis-users <<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px"><span>Dear PostGIS users,</span><div><br></div><div><span>I have a <a href="https://data.jrc.ec.europa.eu/dataset/35fb1231-849b-4017-89d8-9cadbaf9d555" title="raster" target="_blank">raster</a> (of 100m pixel) of forest (val=1) / no forest (val=NULL) and a set of points that always fall inside forest areas.<br><br></span></div><div><span>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<span> (distance B in fig.)</span>, following the same direction. Like so:<br></span></div><div><br><img alt="Screenshot from 2023-10-24 18-19-16.png" src="cid:18b66224535f5ff3ce11"><br><br></div><div><span>I was able to do this with the <span>following</span> steps:</span></div><div><span>1. Vectorize forest map. Note the ST_Union, to eliminate the borders of the raster tiles and avoid these being confused as forest border.</span></div><div><span></span></div><div style="font-family:"Droid Sans Mono","monospace",monospace;font-weight:normal;line-height:19px;white-space:pre-wrap;color:rgb(204,204,204);background-color:rgb(31,31,31)"><span><span style="color:rgb(206,145,120)">SELECT (ST_Dump(ST_Union(geom))).geom::geometry(Polygon,3035)</span></span><div><span style="color:rgb(206,145,120)">FROM (</span></div><div><span style="color:rgb(206,145,120)"> SELECT dp.*</span></div><div><span style="color:rgb(206,145,120)"> FROM storm.forest_2020_europe, </span></div><div><span style="color:rgb(206,145,120)"> LATERAL ST_DumpAsPolygons(rast) AS dp</span></div><div><span style="color:rgb(206,145,120)"> ) As foo</span></div><div><span style="color:rgb(206,145,120)">GROUP BY val</span></div><span><span style="color:rgb(206,145,120)">HAVING val = 1;</span></span></div><div><span></span></div><div><br></div><div><span>2. Create Linestring from point, as starting point, following 300 radians and 1km length.</span></div><div><span></span></div><div style="font-family:"Droid Sans Mono","monospace",monospace;font-weight:normal;line-height:19px;white-space:pre-wrap;color:rgb(204,204,204);background-color:rgb(31,31,31)"><span><span style="color:rgb(206,145,120)">SELECT gid, </span></span><div><span style="color:rgb(206,145,120)"> '</span><span>1km </span><span style="color:rgb(181,206,168)">300</span><span> rad</span><span style="color:rgb(206,145,120)">', </span></div><div><span style="color:rgb(206,145,120)"> ST_MakeLine(geom, ST_Translate(geom, sin(300)*1000, cos(300)*1000)) AS geom</span></div><span><span style="color:rgb(206,145,120)">FROM storm.point;</span></span></div><div><span></span></div><div><br></div><div><span>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 <span>ST_Intersection and ST_Difference, and then ST_Touches to make sure the segment touches the first point, for the first distance and <span><span>touches</span></span> the EndPoint of the first segment to get the second distance.</span><br></span></div><div><br></div><div><span>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.</span></div><div><br></div><div><span>Any hints how could this be accomplished?<br></span></div><span></span><br>Thanks!<br><br>Cheers,<br></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div>
<div style="font-family:Arial,sans-serif;font-size:14px">
<div><div><div>Eloi<br></div></div><div><br></div></div>
<div>
</div>
</div>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
</blockquote></div>
</blockquote></div>