OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.93k stars 2.56k forks source link

Feature Request: Proper proximity calculation using Azimuthal Equidistant projection #5006

Open MarByteBeep opened 2 years ago

MarByteBeep commented 2 years ago

Feature description

(taken from: https://github.com/qgis/QGIS/issues/46509)

Feature description

The result of the current Raster -> Analysis -> Proximity implementation heavily depends on the chosen CRS. It will only work properly if the measured Cartesian and ellipsoid distances on a map are roughly equal. This only occurs when you have a local projection. For instance when using EPSG:28992 (Local Dutch projection) when working on a map of the Netherlands. This is all fine when working on a small map and when you have that projection available, but it quickly becomes very inaccurate when your map is larger, or when your CRS isn't correct.

Based on work done by Spatial Thoughts I created a simple algorithm that iteratively created new features and stores the distance to the original feature in a BURN field. These distances are calculated using Azimuthal Equidistant projection. These new features can then be used to create a more accurate proximity map.

Let me provide an example. Consider the following features (military objects around some local town)

image

Next I run the proximity script, with the following settings

Buffer Distance specifies the distance the buffer is drawn around a feature. In this case I defined a distance of 5000 meters.

Segments specifies how well curves are approximated. Lower is more jaggy.

Iterations specifies how many iso-distances are generated.

Rasterize render resolution specifies the rasterizer resolution in degrees.

You can deselect Open output file after running algorithm of the Proximity Vector layer, as that layer is mostly an intermediate layer. But for the sake of explaining what's going on I'll leave it switched on.

image

When hitting run two layers, Proximity Vector layer and Rasterized, are generated.

First let's focus on the Proximity Vector layer:

image

When inspecting one of the structures, you can see that all rings contain different BURN values, each representing the distance to the original feature. Obviously I could also have named the field: DISTANCE. But whatever :)

image

The other layer Rasterize is the one we're after:

image

When sampling this layer, the value in the band will represent the distance in meters.

image

Code available below as a GitHub gist. Any optimization suggestions or other ideas are more than welcome!

https://gist.github.com/MarByteBeep/e688759bc051a59a5651e0ecd45240ff

MarByteBeep commented 2 years ago

I added automatic rasterization as well. Updated the code/text.