[postgis-users] Distance from Rast and LineString overlay
Nicolas Ribot
nicky666 at gmail.com
Wed Oct 25 05:51:36 PDT 2023
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:
[image: 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:
>>
>> [image: 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 val
>> HAVING 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 geom
>> FROM 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/20231025/3831823d/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/20231025/3831823d/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/20231025/3831823d/attachment-0001.png>
More information about the postgis-users
mailing list