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.

This commit is contained in:
Frank Bausch 2018-04-04 08:50:39 +01:00 committed by Robert Osfield
parent 7541bc9d48
commit ceed09afea

View File

@ -168,6 +168,38 @@ inline void EllipsoidModel::convertLatLongHeightToXYZ(double latitude, double lo
inline void EllipsoidModel::convertXYZToLatLongHeight(double X, double Y, double Z,
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
double p = sqrt(X*X + Y*Y);
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) /
(p - _eccentricitySquared*_radiusEquator*cos_theta*cos_theta*cos_theta) );
longitude = atan2(Y,X);
double sin_latitude = sin(latitude);
double N = _radiusEquator / sqrt( 1.0 - _eccentricitySquared*sin_latitude*sin_latitude);