<div style="font-family: Arial, sans-serif; font-size: 14px;">Hi Nicolas,</div><div style="font-family: Arial, sans-serif; font-size: 14px;"><br></div><div style="font-family: Arial, sans-serif; font-size: 14px;">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.</div><div style="font-family: Arial, sans-serif; font-size: 14px;"><br></div><div style="font-family: Arial, sans-serif; font-size: 14px;">Cheers!<br></div><div style="font-family: Arial, sans-serif; font-size: 14px;"><br></div>
<div class="protonmail_signature_block" style="font-family: Arial, sans-serif; font-size: 14px;">
<div class="protonmail_signature_block-user">
<div><div>Eloi<br></div></div><div><br></div>
</div>
<div class="protonmail_signature_block-proton protonmail_signature_block-empty">
</div>
</div>
<div style="font-family: Arial, sans-serif; font-size: 14px;"><br></div><div class="protonmail_quote">
------- Original Message -------<br>
On Wednesday, October 25th, 2023 at 14:51, Nicolas Ribot via postgis-users <postgis-users@lists.osgeo.org> wrote:<br><br>
<blockquote class="protonmail_quote" type="cite">
<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 height="477" width="563" alt="Screenshot 2023-10-25 at 14.48.22.png" class="proton-embedded" src="cid:ii_lo5r4jei1"><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 class="gmail_attr" dir="ltr">On Wed, 25 Oct 2023 at 11:28, Nicolas Ribot <<a href="mailto:nicky666@gmail.com" rel="noreferrer nofollow noopener" target="_blank">nicky666@gmail.com</a>> wrote:<br></div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><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 class="gmail_attr" dir="ltr">On Tue, 24 Oct 2023 at 18:27, Eloi Ribeiro via postgis-users <<a target="_blank" href="mailto:postgis-users@lists.osgeo.org" rel="noreferrer nofollow noopener">postgis-users@lists.osgeo.org</a>> wrote:<br></div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div style="font-family:Arial,sans-serif;font-size:14px"><span>Dear PostGIS users,</span><div><br></div><div><span>I have a <a target="_blank" title="raster" href="https://data.jrc.ec.europa.eu/dataset/35fb1231-849b-4017-89d8-9cadbaf9d555" rel="noreferrer nofollow noopener">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" class="proton-embedded" 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 target="_blank" href="mailto:postgis-users@lists.osgeo.org" rel="noreferrer nofollow noopener">postgis-users@lists.osgeo.org</a><br>
<a target="_blank" rel="noreferrer nofollow noopener" href="https://lists.osgeo.org/mailman/listinfo/postgis-users">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
</blockquote></div>
</blockquote></div>
</blockquote><br>
</div>