Bug with grdtrack and -Ejustification and -C??

[Esteban82]

I strangely got some NaN with the third command where I only revert the sense of the main profile (-E)

gmt grdmath -R0/4/1/2 -I1 X = grid.nc
gmt grdtrack -ECT/CB -Ggrid.nc -C4/1/1
gmt grdtrack -ECB/CT -Ggrid.nc -C4/1/1

Output:

> Cross profile number -L"0-0" at 2/2 az=090.0
0       2       -2      90      0
1       2       -1      90      1
2       2       0       90      2
3       2       1       90      3
4       2       2       90      4
> Cross profile number -L"0-1" at 2/1 az=090.0
0       1       -2      90      0
1       1       -1      90      1
2       1       0       90      2
3       1       1       90      3
4       1       2       90      4
grdtrack [WARNING]: Some points along your profiles were outside the grid domain(s).
> Cross profile number -L"0-0" at 2/1 az=090.0
4       1       -2      270     4
3       1       -1      270     3
2       1       0       270     2
1       1       1       270     NaN
0       1       2       270     NaN
> Cross profile number -L"0-1" at 2/2 az=090.0
4       2       -2      270     NaN
3       2       -1      270     3
2       2       0       270     2
1       2       1       270     1
0       2       2       270     0
grdtrack [WARNING]: Some input points were outside the grid domain(s).

[PaulWessel]

This might be a simpler example of what @seisman pointed out (those towers near TL BL corners for grids subjected to blending. Darn, can you remind moe of that issue, @seisman ?

[joa-quim]

GenericMappingTools/pygmt#2511 ?

[PaulWessel]

Yes. I will try to debug @Esteban82 example - it is annoying but I spent a lot of time on the grid version already...

[PaulWessel]

This is what I have learned. grdtrack -E creates the profiles, then the cross track machinery computes the angle of the cross profiles to be 180. Then we can sincosd (angle, &s, &c) and get cos = -1 (good) and sin = -1.2246467991473532E-16 (oops). Then, coordinates are computed and instead of getting y = 1 we get 0.99999999999999988 and then it is determined by BCR that the point is outside the grid and we get the default NaN answer.

I think for some reason sincosd fails since it actually calls sincos (angle*D2R) and hence roundoff is likely. So where do we intercept? Make the macro check for angle being exact multiple of 90 and set those returns manually? This is in many places... Suggestions - gotta take a break.

[joa-quim]

Can the BCR be made a bit more tolerant on what is outside? If “sufficiently close accept it”? Go make the break(s). If necessary, this can stay as a bug (again, our release process is so heavy that it prevents more often releases. I would drop most of the steps in the check list). Only really important bug to fix is that padding effect that makes grdimage crash. From: Paul Wessel ***@***.***> Sent: Sunday, December 17, 2023 5:46 PM To: GenericMappingTools/gmt ***@***.***> Cc: Joaquim Manuel Freire Luís ***@***.***>; Comment ***@***.***> Subject: Re: [GenericMappingTools/gmt] Bug with grdtrack and -Ejustification and -C?? (Issue #8194) This is what I have learned. grdtrack -E creates the profiles, then the cross track machinery computes the angle of the cross profiles to be 180. Then we can sincosd (angle, &s, &c) and get cos = -1 (good) and sin = -1.2246467991473532E-16 (oops). Then, coordinates are computed and instead of getting y = 1 we get 0.99999999999999988 and then it is determined by BCR that the point is outside the grid and we get the default NaN answer. I think for some reason sincosd fails since it actually calls sincos (angle*D2R) and hence roundoff is likely. So where do we intercept? Make the macro check for angle being exact multiple of 90 and set those returns manually? This is in many places... Suggestions - gotta take a break. — Reply to this email directly, view it on GitHub<#8194 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAEDF2N2HO4T2LOUO55KXQ3YJ4VU3AVCNFSM6AAAAABAWXO2GOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNJZGIZTIMZSHE>. You are receiving this because you commented.Message ID: ***@***.******@***.***>>

[joa-quim]

Ah, and I have another dinner tonight. From: Paul Wessel ***@***.***> Sent: Sunday, December 17, 2023 5:46 PM To: GenericMappingTools/gmt ***@***.***> Cc: Joaquim Manuel Freire Luís ***@***.***>; Comment ***@***.***> Subject: Re: [GenericMappingTools/gmt] Bug with grdtrack and -Ejustification and -C?? (Issue #8194) This is what I have learned. grdtrack -E creates the profiles, then the cross track machinery computes the angle of the cross profiles to be 180. Then we can sincosd (angle, &s, &c) and get cos = -1 (good) and sin = -1.2246467991473532E-16 (oops). Then, coordinates are computed and instead of getting y = 1 we get 0.99999999999999988 and then it is determined by BCR that the point is outside the grid and we get the default NaN answer. I think for some reason sincosd fails since it actually calls sincos (angle*D2R) and hence roundoff is likely. So where do we intercept? Make the macro check for angle being exact multiple of 90 and set those returns manually? This is in many places... Suggestions - gotta take a break. — Reply to this email directly, view it on GitHub<#8194 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAEDF2N2HO4T2LOUO55KXQ3YJ4VU3AVCNFSM6AAAAABAWXO2GOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNJZGIZTIMZSHE>. You are receiving this because you commented.Message ID: ***@***.******@***.***>>

[PaulWessel]

Lamely, this fixes the issue

angle = 90.0 - az_cross - deviation;	/* The direction */
sincosd (angle, &sa, &ca);	/* Trig on the direction */
i_angle = rint (angle);
if (i_angle == 0 || i_angle == 360) sa = 0.0, ca = 1.0;
else if (i_angle == 90) sa = 1.0, ca = 0.0;
else if (i_angle == 180 || i_angle == -180) sa = 0.0, ca = -1.0;
else if (i_angle == 270 || i_angle == -90) sa = -1.0, ca = 0.0;