[PROJ] Geodesic buffers/offsets?

[jmfluis at gmail.com]

Hi,

I have taken a shot and implemented a Geodesic buffer in
Julia+Proj+GDAL+GMT. 
https://github.com/GenericMappingTools/GMT.jl/blob/master/src/proj_utils.jl#
L237

It works by interpolating the line segments at 1/4 circle radius (buffer
width) and calculating the union them. At the end there is some filtering
and Douglas-Peucker line simplification.
Example fig here
https://forum.generic-mapping-tools.org/t/how-to-use-gmt-spatial-sb-correctl
y/1725/13?u=joaquim

I have no tool to assess the accuracy of the final solution. GMT comes close
but it's point-to-line function works only in spherical Earth. Computing a
buffer in the spherical approximation shows that the average distance from
points to line when the asked buffer width is 111111 m turns out to be (for
the example used) ~111062 m. Not bad, I think.

See some example figures and discussion in the GMT issue.
https://github.com/GenericMappingTools/gmt/issues/5295

Joaquim

-----Original Message-----
From: PROJ <proj-bounces at lists.osgeo.org> On Behalf Of Dave Knopp
Sent: Wednesday, May 5, 2021 5:01 PM
To: Elmar Brockmann via PROJ <proj at lists.osgeo.org>; nyall.dawson at gmail.com
Subject: [PROJ] Geodesic buffers/offsets?


Nyall Dawson nyall.dawson at gmail.com
Tue May 4 15:23:09 PDT 2021
> Something which I get asked on a recurring basis by users is how to 
> create geometry buffers using a real-world distance on a geographic 
> dataset. E.g. creating a 100km buffer on a linestring in such a way 
> that the geodesic distances between the buffer and original line are 
> always 100km.


Assuming there is some desired 'target resolution' associated with the
application (e.g. display pixel sizes, vector polyline sampling interval,
etc.) then perhaps a simple dilation algorithm might work for you.


E.g.

Allocate a working grid associated entire range of lon/lat values at the
resolution you want. Us full ellipsoid range or restrict range of values if
that's appropriate for application cases. With each cell, associate data
structure comprising a scalar and 3D vector.

Initialize the vector portion of the grid cells with corresponding ECEF
locations for relevant surface - e.g. on the ellipsoid, on terrain, etc).
Fill the scalar portion with NaN or whatever.

During iterations, the scalar will be set to the distance metric from the
source points (e.g. distance along geodesic, distance over undulating
terrain surface, etc).

Use lon/lat of source points to find associated grid cell then fill the
scalar part of those grid cells (e.g. using line drawing or other
interpolation if needed). For these cells, set the scalar distance to 0.

Perform a modified 'dilation' algorithm that iterates over cells just set in
last iteration and determine/sets scalar value on unset adjacent cells.

Assign the new scalar values based on the magnitude of *vector* difference
between cell vector location values. For potential collisions, use 'shortest
wins' strategy.

Iterate on this until grid is full - or optionally stop if all new distances
in iteration exceed the target 'max dist' value.

For vector representation of the buffer boundary, run a contour following
algorithm on the grid scalar values.

That's all off-the-cuff, so maybe not complete nor error free.

Before thinking this is "too much computation" consider the underlying
complexity of other algorithms - which is also high (at comparable detail).
Also this has considerable additional flexibility which might boost reuse
potential.

Dave
_______________________________________________
PROJ mailing list
PROJ at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj