The inverse operation to GeoDisplacement is GeoDestination, so we can construct an inverse projection using it:

polar[{x_, y_}] := GeoDisplacement[{Sqrt[x^2 + y^2], ArcTan[y, x]/Degree}];
InverseEAE[coords_, c_] := GeoDestination[c, polar[coords]];

Now you can do:

In[]:= EAE[Here, GeoPosition[{0, 0}]]
Out[]= {-7.55175*10^6, 6.34236*10^6}

In[]:= InverseEAE[%, GeoPosition[{0, 0}]] == Here
Out[]= True

Jose.

Hi Fernando,

I think this is the difference between using a projection with a spherical reference model and using it with an ellipsoidal reference model of the Earth. The projection engine of the Wolfram Language contains formulas with spherical models for all projections, and with ellipsoidal models for some projections, like the Mercator or the Transverse Mercator projections.

The ellipsoidal azimuthal equidistant (EAE) projection can be easily obtained from GeoDisplacement, which uses the "ITRF00" model by default:

cartesian[GeoDisplacement[{r_, a_}]] := r {Sin[a Degree], Cos[a Degree]};
EAE[p_, c_] := cartesian[GeoDisplacement[c, p]];

where p is the GeoPosition of the point you want to compute the projection, and c is is the projection center. Then for example, using the default center at {0, 0}:

In[]:= EAE[Here, GeoPosition[{0, 0}]]
Out[]= {-7.55175*10^6, 6.34236*10^6}

with the result in meters. Compare those numbers to the result from the spherical projection with sphere of average radius:

In[]:= R = GeodesyData["ITRF00", "MeanRadius"]
Out[]= Quantity[6.37101*10^6, "Meters"]

In[]:= GeoGridPosition[Here, {"AzimuthalEquidistant", "ReferenceModel" -> R}][[1]]
Out[]= {-7.53793*10^6, 6.35277*10^6}

Differences are of the order of 10 kilometers, as usual when comparing spherical and ellipsoidal models of the Earth.

Jose.