mapproject with -Glon/lat/unit doing strange things

Added by John about 1 year ago

My favorite GMT command is mapproject. Something funny is going on with the -Glon/lat/unit option. MacOSX 10.11.6, GMT Version 6.0.0_r19216.

Here's an example:

pointfile=pts.xy

cat << EOF >$pointfile
-109.4 39.8
-109.0 40.3
-108.9 40.1
EOF

mapproject $pointfile -G-109.1/40/+k

This yields monotonically increasing distances:

-109.4 39.8 33.9312819658
-109 40.3 68.3150139724
-108.9 40.1 88.67545481

But that isn't what was asked. Monotonically increasing (accumulating) distances should only occur when lon/lat isn't specified (Use -G[<unit>]+a to get accumulated distances along track [Default].)

The distance to the first point is correct. Distances to second and third points should be:

echo -109.1 40 -109.0 40.3 | mapproject -G+k+
-109.1 40 -109 40.3 34.3837320067

echo -109.1 40 -108.9 40.1 | mapproject -G+k+
-109.1 40 -108.9 40.1 20.3604408375

Looks like something is erroneously causing computation of accumulated distances when -Glon/lat parameters have been specified.


Replies (6)

RE: mapproject with -Glon/lat/unit doing strange things - Added by Paul about 1 year ago

I will see why it gets confused but note that the correct syntax should be

gmt mapproject pts.xy -G-109.1/40+u+k

which indeed gives the correct ellipsoidal incremental output (after fixing a related parsing issue in r19235.)

So for the old syntax (I am looking at GMT4 here) to get incremental distances you would do -G-109.1/40/-k and to get ellipsoidal distances you would use capital K as unit. I am not sure that your /+k was ever meant to do what you wanted to do. Where did you see that syntax?

RE: mapproject with -Glon/lat/unit doing strange things - Added by John about 1 year ago

I saw it in the GMT mapproject 5.3.0_r15370 man page. Under units there was a line "prepend + to perform exact geodesic calculations." Actually, this line still exists in the latest man page (6.0.0_r18654).

It is evident that the latest mapproject man page has changed the -G option description quite a bit.

However, I'm not able to get this next construct to work properly (to compute the geodesic distance between two points):

echo -100 41 -101 40 | mapproject -G+u+k+v
-100 41 -101 40 0

hopefully, I'm using the options properly. The result returned is zero.

On the other hand, this construct yields a non-zero distance

echo -100 41 | mapproject -G-101/40+u+k
-100 41 139.699

Something is still not quite right; the first construct should yield the same result. Hope this helps!

RE: mapproject with -Glon/lat/unit doing strange things - Added by Paul about 1 year ago

Yep, the wrong column was output. Try r19243. Note that +a is default so if you have many pairs it will add up the distances between the pairs. To just get increments per line add +i.
John, if you have time it would be great if you could put together a test script that demonstrates the many uses you put mapproject to. That way we can add it to our tests and make sure things pass and not reopen old wounds...

RE: mapproject with -Glon/lat/unit doing strange things - Added by John about 1 year ago

In the next day or two, I'll try to put together a simple script to help with future testing of mapproject.

I updated to r19243, this morning, but still get a null (zero) result for the command:

echo -100 41 -101 40 | mapproject -G+u+k+v

Quite a few of my scripts use this kind of command. Thanks!

RE: mapproject with -Glon/lat/unit doing strange things - Added by Paul about 1 year ago

Aargh, forgot to commit at home. This was fixed. Try now (r19245).

RE: mapproject with -Glon/lat/unit doing strange things - Added by John about 1 year ago

Yea! Thanks Paul. Sorry to start your morning with such a minor thing. Looks and works great! I use that construct often - very much appreciate!

(1-6/6)