FGFS-4.1/OpenSceneGraph
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

When testing XYZ/ECEF coordinate values that convert to LLA points with latitude = 90.0 or latitude = -90.0, the heights calculated by EllipsoidModel::convertXYZToLatLongHeight have substantial errors. With the provided change, points at the pole are handled more correctly. This has been tested against US NGA GoldData.

Browse Source
This commit is contained in:
Frank Bausch 2018-04-04 08:50:39 +01:00 committed by Robert Osfield
parent 10f000f043
commit 72ab22e539
1 changed files with 32 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
include/osg/CoordinateSystemNode
View File

@ -168,6 +168,38 @@ inline void EllipsoidModel::convertLatLongHeightToXYZ(double latitude, double lo
inline void EllipsoidModel::convertXYZToLatLongHeight(double X, double Y, double Z, inline void EllipsoidModel::convertXYZToLatLongHeight(double X, double Y, double Z,
double& latitude, double& longitude, double& height) const double& latitude, double& longitude, double& height) const
{ {
// handle polar and center-of-earth cases directly.
if (X != 0.0)
longitude = atan2(Y,X);
else
{
if (Y > 0.0)
longitude = PI_2;
else if (Y < 0.0)
longitude = -PI_2;
else
{
// at pole or at center of the earth
longitude = 0.0;
if (Z > 0.0)
{ // north pole.
latitude = PI_2;
height = Z - _radiusPolar;
}
else if (Z < 0.0)
{ // south pole.
latitude = -PI_2;
height = -Z - _radiusPolar;
}
else
{ // center of earth.
latitude = PI_2;
height = -_radiusPolar;
}
return;
}
}
// http://www.colorado.edu/geography/gcraft/notes/datum/gif/xyzllh.gif // http://www.colorado.edu/geography/gcraft/notes/datum/gif/xyzllh.gif
double p = sqrt(X*X + Y*Y); double p = sqrt(X*X + Y*Y);
double theta = atan2(Z*_radiusEquator , (p*_radiusPolar)); double theta = atan2(Z*_radiusEquator , (p*_radiusPolar));
@ -179,7 +211,6 @@ inline void EllipsoidModel::convertXYZToLatLongHeight(double X, double Y, double
latitude = atan( (Z + eDashSquared*_radiusPolar*sin_theta*sin_theta*sin_theta) / latitude = atan( (Z + eDashSquared*_radiusPolar*sin_theta*sin_theta*sin_theta) /
(p - _eccentricitySquared*_radiusEquator*cos_theta*cos_theta*cos_theta) ); (p - _eccentricitySquared*_radiusEquator*cos_theta*cos_theta*cos_theta) );
longitude = atan2(Y,X);
double sin_latitude = sin(latitude); double sin_latitude = sin(latitude);
double N = _radiusEquator / sqrt( 1.0 - _eccentricitySquared*sin_latitude*sin_latitude); double N = _radiusEquator / sqrt( 1.0 - _eccentricitySquared*sin_latitude*sin_latitude);