search

OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.74k stars 788 forks source link

NaN result for inverse transform on ellipsoid with very close points #129

Closed proj4-bot closed 9 years ago

proj4-bot commented 9 years ago

Reported by nnseva on 7 Dec 2011 14:15 UTC I've passed very close coordinates with difference at fifth digit after dp in seconds to invgeod and received a NaN result:

echo 56d16'19.00001 43d58'35.00001 56d16'19.00002 43d58'35.00002 | invgeod +ellps=WGS84 -
-2147483648d-2147483648'-nan"   -2147483648d-2147483648'-nan"   -nan

What I am expecting: Basing on passed values it should be 0,0,0

Really the invgeod filters input (see # define DTOL at the top of the geod_inv.c file), but it is not enough as we can see.

The problem has been initially found and reported to pyproj project v.1.2.9/trunk. It also has been reproduced on invgeod v.4.7.1. As I can see in the geod_inv.c source, it should still remain in trunk version.

Last Modified by warmerdam on 26 Jun 2013 06:04 UTC

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

proj4-bot commented 9 years ago

Comment by nnseva on 8 Dec 2011 07:23 UTC I suppose that the DTOL value increased to 1e-8 works the problem around. The zero result then is returned till about 0.15 meters between points in my case, I think this is good enough (results below are from pyproj compiled with fixed DTOL):

>>> geod.inv(43.9763413300,56.2720659176,43.9763419000,56.272067000)
(0.0, 0.0, 0.0)
>>> geod.inv(43.9763413300,56.2720659176,43.976342000,56.272067000)
(19.00427040948362, -160.99572903328846, 0.13428452977133082)

BUT: The strong mathematical base should be used here to substantiate input filter algorithm instead of empiric like this.

proj4-bot commented 9 years ago

Comment by warmerdam on 3 Mar 2012 22:51 UTC Note that I'm not overly enthusiastic about working on the geodesic code in PROJ.4. There are better implementations elsewhere. See GeodesicCalculations

I'll leave this ticket open in case someone has a fix to offer, but I'm not planning to look into the issue.

proj4-bot commented 9 years ago

Modified by warmerdam on 26 Jun 2013 06:04 UTC

proj4-bot commented 9 years ago

Comment by warmerdam on 26 Jun 2013 06:07 UTC The geod algorithm has been reimplemented in PROJ.4 based on Charles Karney's code and now we get:

echo 56d16'19.00001 43d58'35.00001 56d16'19.00002 43d58'35.00002 | invgeod +ellps=WGS84 -
6d57'55.184"    -173d2'4.815"   0.019

I think this is a reasonable response. Closing.