Proper proximity calculation using Azimuthal Equidistant projection

[MarByteBeep]

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)

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.

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

First let's focus on the Proximity Vector layer:

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 :)

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

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

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

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

Additional context

No response

[gioman]

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.

@MarByteBeep this how it works https://gdal.org/programs/gdal_proximity.html

You should extend this request also to GDAL, and see if there is any room for improvement there.

[MarByteBeep]

@gioman sorry, but how do I do that? I also noticed that I forgot the trivial case: distance is 0 meters. I'll fix the code shortly.

[gioman]

sorry, but how do I do that?

@MarByteBeep https://github.com/OSGeo/gdal/issues

[MarByteBeep]

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

[nyalldawson]

@MarByteBeep you should consider publishing your algorithm as a plugin, it's a great fit for a solution of this nature

[MarByteBeep]

@nyalldawson haha, I just downloaded QGIS two days ago, so I have NO clue on how to make one, yet :) Any help would be greatly appreciated.

[roya0045]

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.

CRS selection is always critical. That's pretty much always the case for everything, from simple length calculations to area. I'm not sure how much of it should be improved...