command-line geod fails when listing points along the equator

[proj4-bot]

Reported by mikebanahan on 3 Aug 2009 07:03 UTC Using geod to list points on a great circle route around the equator causes NaN values to be produced, e.g

geod +ellps=WGS84 +n_S=25 +lat_1=0dN +lon_1=-2dW +lat_2=0dN +lon_2=2dE -f %f

produces

0.000000 -2.000000

nan nan

nan nan

... etc

Looking at the code there are two clear problems. The first is not relevant to this but should be fixed. Line 18 in geod_for.c reads

if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {

whereas I believe it SHOULD be

if ((merid = (fabs(sina12 = sin(al12)) < MERI_TOL))) {

i.e instead of assigning the fabs value to merid, it should assign the truth value.

That then prompts what appears to be correct behaviour on north-south geodesics along the 0 longitude meridian

The cause of NaN is at lines 41-46:

41 if (merid) s1 = HALFPI - th1;

42 else {

43 s1 = (fabs(M) >= 1.) ? 0. : acos(M);

44 s1 = sinth1 / sin(s1);

45 s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);

46 }

Where on line 43 s1 becomes zero and line 44 immediately produces divide-by-zero

As a result I've had to switch to using Geographic Library for great-circle work, which is a shame. I don't understand the code well enough (too few comments here and what are to me, opaque variable names) to fix it. It doesn't help that I'm not a good enough geographer/mathematician to understand the algorithm of course.

Migrated-From: https://trac.osgeo.org/proj/ticket/48

[proj4-bot]

Comment by karney on 30 Nov 2012 19:34 UTC Note that the patch supplied with ticket #197 fixes this problem.

[proj4-bot]

Comment by warmerdam on 26 Jun 2013 06:12 UTC The code from #197 (reimplementing geod) is now checked in and the reported request gives the following which I suppose is appropriate:

src/geod +ellps=WGS84 +n_S=25 +lat_1=0dN +lon_1=-2dW +lat_2=0dN +lon_2=2dE -f %f
0.000000    -2.000000
0.000000    -1.840000
0.000000    -1.680000
0.000000    -1.520000
0.000000    -1.360000
0.000000    -1.200000
0.000000    -1.040000
0.000000    -0.880000
0.000000    -0.720000
0.000000    -0.560000
0.000000    -0.400000
0.000000    -0.240000
0.000000    -0.080000
0.000000    0.080000
0.000000    0.240000
0.000000    0.400000
0.000000    0.560000
0.000000    0.720000
0.000000    0.880000
0.000000    1.040000
0.000000    1.200000
0.000000    1.360000
0.000000    1.520000
0.000000    1.680000
0.000000    1.840000
0.000000    2.000000