<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">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>