Mapproject -E and -I

Hate to ruin someone’s weekend! This is for Version 6.5.0_2da67da_2023.11.22.

It seems that mapproject coordinate transformation: lon,lat,height to X,Y,Z; then, running the reverse transformation : X,Y,Z to lon,lat,height yields something strange. Consider this example at the north pole.

Using these commands (presuming default WGS84 ellipsoid):

> echo "0 90 0" | gmt mapproject -E --FORMAT_FLOAT_OUT=%.5f
0.00000	0.00000	6356752.31425

That looks about right, given that the semi-minor axis of WGS-84 is something like: 6356752.314245 m

Now feed those (X,Y,Z) coordinates back into mapproject with -I:

> echo "0.00000 0.00000 6356752.31425" | gmt mapproject -E -I --FORMAT_FLOAT_OUT=%.5f
0.00000	-90.00000	-6399593.62576

Oops! Somehow, the result goes to the south pole latitude with a height that puts it close to the center of the Earth, 42.8 km above the geocenter.

Excellent. I think I have never used -E but pretty sure I was the one who implemented it back in the last century… Will try to look tomorrow.

Think we got it (now in master):

echo 0 -90 0 | gmt mapproject -E --FORMAT_FLOAT_OUT=%.5f | gmt mapproject -E --FORMAT_FLOAT_OUT=%.5f -I
0.00000 -90.00000 -0.00000

Thanks, Pål, hope it didn’t ruin your weekend. After updating to Master, I confirm that this works just fine!

Did not take long so no worries.

Had to read a bit about ECEF. As I understand, going from lon,lat,alt -> ecef is easy. The inverse is not, as you need (many) iterations and convergence to get precision. What does gmt do?

We do a ~2 step iteration. That should be good enough (now that special bug was fixed)

1 Like


Doing this in cs2cs gives a different result:

With cs2cs:

$ echo 10 60 | cs2cs EPSG:4326 EPSG:4978
3140936.41	5440261.45 1100248.55

With mapproject:

echo 10 60 0 | gmt mapproject -E
3148533.38439	555171.38533	5500477.13394

Should not these be the same?

and it is :slight_smile:

C:\v>echo 60 10 | cs2cs EPSG:4326 EPSG:4978
3148533.38      555171.39 5500477.13

Thanks @Joaquim. Order matters…