Ticket #48 (closed defect: fixed)

Opened 5 years ago

Last modified 17 months ago

command-line geod fails when listing points along the equator

Reported by: mikebanahan Owned by: warmerdam
Priority: major Milestone:
Component: CommandlineUtilities Version: unspecified
Keywords: geod, geodesic Cc:

Description

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.

Change History

Changed 2 years ago by karney

Note that the patch supplied with ticket #197 fixes this problem.

Changed 17 months ago by warmerdam

  • status changed from new to closed
  • resolution set to fixed

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
Note: See TracTickets for help on using tickets.