Bug #1248

project -N yields unexpected results

Added by Dan 4 months ago. Updated 4 months ago.

Status:In ProgressStart date:2018-06-04
Priority:NormalDue date:
Assignee:Paul% Done:

50%

Category:-
Target version:-
Affected version:4.5.16 Platform:

Description

I am getting unexpected results with project -N and wonder if there is a bug in the code, in the
manual page on which my expectations are based, or in-between my ears. Namely, I get the results I expect
when I use the spherical trig math (no -N), but I have to subtract 90deg from the azimuth
to get the expected answer when using the Cartesian math (-N).

A simple example where I expect the output PQ to equal the input XY because I specify -A090 and -C0/0.

This works with project (no -N):

% echo 2 3 "<--inputXY -->outputPQRS" | $PROJECT -Fxyzpqrs -C0/0 -A090
2 3 2 3 2 1.22435767802e-16

This gives unexpected PQRS results with project -N:

% echo 2 3 "<--inputXY -->outputPQRS" | project -Fxyzpqrs -C0/0 -A090 -N
2 3 <--inputXY -->outputPQRS -3 2 -1.83691e-16 3

But if I change to -A000, it yields my expected answer.

% echo 2 3 "<--inputXY -->outputPQRS" | project -Fxyzpqrs -C0/0 -A000 -N
2 3 <--inputXY -->outputPQRS 2 3 2 0

I observe this behavior at a host of azimuths and replacing -Aaziumuth with
-Ebx/by. I run this on a CentOS box using:

project 4.5.18_r10406 [64-bit] - Project data onto line or great circle, generate track, or translate coordinates

I also observe this behavior using:

project [core] 6.0.0_r20206 [64-bit] - Project data onto lines or great circles, generate tracks, or translate coordinates

Also, in the GMT5 SVN version, I note that output of z columns occurs only for numeric values, not alphanumeric columns
as in GMT4. This can be seen in my example, above. I don't know if this is a bug or a feature, but I prefer the GMT4 behavior.

t.pdf - Cartesian rotations (9.11 KB) Paul, 2018-06-05 17:28

Associated revisions

Revision 20207
Added by Paul 4 months ago

Partial fix for issue #1248 regarding z string output

History

#1 Updated by Paul 4 months ago

  • Status changed from New to In Progress
  • Assignee set to Paul
  • % Done changed from 0 to 50

How text is handled in GMT 6 is different but not outputting it (like in your example) was a bug - now fixed in r20207. In GMT 6 all table records can be of the form

[<numerical columns>][<trailing text>]

whereas some GMT 4 and 5 modules might produce output that has numerical columns followed by text followed by other numerical columns. THis makes the output extremely non-portable in other context and we now require module output to follow the above format. This means all GMT module now has the ability to carry along trailing text, not just a few (e.g., project, mapproject, gmtselect). Those who have workflows that cannot change to handle the new and more predictable format will have to use GMT 4 or update their workflow. For instance of your "z" actually consists of a few numerical columns and a string then the order on the output would place those z-columns in according to the -F order but the z string would be trailing text.

I will look into the unexpected results next.

#2 Updated by Paul 4 months ago

I am looking at the Cartesian case. See PDF and view from bottom panel towards top panel. The black dot is your (x,y) point (3,2) in the original (black) coordinate system while the projected point (p,q) is the coordinates relative to the new axis (red). The red dot is the point (r,s) which is the position in the (x, y) coordinate system of the point on the profile (q = 0 path) closest to (x, y). Obviously, for no rotation (bottom panel) the two coordinate systems are the same and you get 3,2 back. As you rotate the coordinate system (i.e., 15, 30 etc going up the column) the values change. At the top (90 degree rotation, your rotated x-axis (p-axis) is pointing in the negative y-direction and hence your p-coordinate is the negative y-coordinate and the q coordinate is the positive z coordinate. Hence (3,2) rotates to give (2,3). This looks correct to me.

#3 Updated by Dan 4 months ago

Thanks for fixing and explaining the z-string behavior in GMT 6.

Your PDF of the Cartesian case is very helpful. The alpha annotated in your examples relates to the -Aazimuth
value, correct? From the following two parts of the project manual page:

When moving in the increasing (p) direction, (toward B or in the azimuth direction),...

and

Flat Earth (Cartesian) coordinate transformations can also be made. Set -N and remember that azimuth is clockwise from North (the y axis),...

my expectation was that for -Aazimuth=0, p coincides with the y-axis and q coincides with the negative x-axis.
In the bottom panel of your graphic where azimuth=0, you show p coinciding with the x-axis and q coinciding with the y-axis,
which seems to disagree with the excerpts from the man-page [increasing p is in the azimuth direction, and azimuth is
clockwise from North (the y axis)].

Thoughts? My expectation also works when you remove -N, which seems a desirable similarity to maintain when folks like
me are working both in geographic and projected coordinate systems.

#4 Updated by Paul 4 months ago

Let me have Walter comment on this - he wrote the earliest version that defined p,q. If GMT 4 and 6 gives the same then this goes back some version so would be good to fix and/or improve documentation.

Also available in: Atom PDF