Greenspline for global data - spurious duplicate?

Q&A
1

Hi
I am trying to interpolate/grid irregularly spaced lon-lat-value data on the sphere in a way that has to be accurate at the poles (hence surface no good) and smooth (hence sphinterpolate no good - at least I can’t get it work properly). This is why I am trying greenspline which, from what I understand, should be able to do this, using something like

greenspline -Rg -Sp -Z3 -I1 -Gtmp.grd tmp.dat
(from the manual)

I have a bunch of globally distributed data, but things fail because the code seems to be complaining, to me, erroneously about duplicate location entries, even though they seem to be in different locations. I am probably missing something obvious, but for example if tmp.dat is

43.782 41.121 2.79
202.5 47.25 - 4.1

I get

image
image1127×218 14.4 KB

Thanks much for your help in advance!

Best

T

2

Here’s what I found. Your two points dist 89.999999918622919 degrees (computed geodetically by gmtmap_geodesic_dist_degree()), which later is taken the cos() and the result, 1.4202979405388712e-09, is smaller than the threshold (1e-8) used to determine when points are considered to overlap. If you change the the point coordinates just a tinny bit (e.g. 202.5 to 202.5001), the error goes away.

Don’t know what would be a good fix for this.

3

Thanks! So this is some sped up proximity check? I’ve also noticed strange behavior for sphinterpolate where numerous datapoints ar flagged as “identical” when they are numerically distinct.

The distance in radias

echo 43.782 41.121 202.5 47.25 | gawk -f ~/awk/dist_on_sphere.awk
1.56435210975521666e+00
is quite different from zero, but I presume it’s too close to pi/2 for something to be messed up?

Where is that check? I would assume the gmtmap_geodesic_dist_degree works fine.

4

gmt_map.c#L6517 calls
return (cosd (gmtmap_geodesic_dist_degree (GMT, lonS, latS, lonE, latE)));

note the cosd() call on result. That’s why the result is a number < 1e-8.

Me too (it’s likely used by several programs) but didn’t check.