Message Boards Message Boards

1
|
7341 Views
|
3 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Use World AzimuthalEquidistant instead of Sphere AzimuthalEquidistant?

Posted 8 years ago

Hi. I am working with some old geographic data and it was all manipulated in ArcMap and saved with a World Azimuthal Equidistant projection. I added more data from a GPS device, and it is currently in Equirectangular coordinates, and tried to transform it using

GeoGridPosition[GeoPosition[#],{"AzimuthalEquidistant","Centering"->mycenter }]&/@moredata

but it seems that everything is off by a tiny bit, compared to the ArcMap coordinates. I put my new data in ArcMap and if I use the Spherical Azimuthal Equidistant projection then ArcMap and Mathematica give the exact same result. My question is: Is there a way to force Mathematica to use the "World" instead of the "Shperical" Azimuthal Equidistant projection? Thanks! Fernando

POSTED BY: Fernando Perez
3 Replies

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.

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.

Posted 8 years ago

Thank you José! I've worked out all the points in the dataset with your trick. This is fantastic. Now, is there an "inverse function"? I want to move the old data to spherical, so it would mean to transform everythin from ellipsoidal to shperical. If I understand correctly, you showed me cartesian->ellipsoidal. And GeoProjection does spherical->cartesian. The inverse function would start ellipsoidal and end in shperical (via cartesian if necessary). Also. There are a bunch of polygons in the old data set too. Is there a GeoProjection that takes the ellipsoidal (and not "world" as in ArcMap -sorry for that-) projection to the spherical projection and/or viceversa? I am asking this only because some times I do not need the exact point, but I need a picture, and some of these "polygons" are filled curves, lines, etc so changing points one by one would be a mess. Thanks a lot. Fernando

POSTED BY: Fernando Perez
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract