How to properly generate distance map?

Question

This question is a related to Proximity (raster distance) using ellipsoidal distance

What I'm trying to do is to generate a map of distances to roads. At any given point on a map, I want to be able to instantly say: "you are X kilometres away from a main road".

So my initial thoughts were, get road information from OSM, or draw them, and generate a distance map using Raster -> Analysis -> Proximity (Raster Distance) (after rasterizing that road vector layer).

The problem however seems to be that the Proximity algorithm isn't aware of any CRS. It simply calculates the distances naively without any transformation. Consider the following trivial example

That's one long road all the way from Africa to Russia in one "straight" line (at least on the map). I rasterized this road and then used the Proximity tool with a cut-off of 500 km, and georeferenced coordinates, and this is what I got:

So clearly this tool does not take into account any stretching of the map. Which results in a straight line with a uniform distance gradient. This means that depending on the position on the map, the sampled distance on the proximity layer is reasonably correct near the equator and completely off the more north or south you go. Which obviously makes sense if the generated distance map is based on Cartesian calculations and not ellipsoid calculations.

Is there any way for me to warp, transform or modify this distance map so that it more closely resembles the CRS I'm using (EPSG:3857)? Accuracy isn't a big problem; it's okay to be a few hundred meters off. As long as I'm able to process the road network of the large portion of the Earth in one go, I'm very happy.

Answer 1

For those battling the same issue in the future: I came up with a new solution, that I posted on GitHub.

TL;DR: This is the code

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

The long read:

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

After running 'the long road' example posted in the original question with the new algorithm, I ended up with something more intuitively correct: