GeoPoint distanceTo

if getSRS() is projected, this method calculate the distance of (dx, dy, dz).

but if getSRS() is geographic, this method only calculate the distance of (dx, dy). dz is ignored.

is it right?

You are correct - that is what the code does. It should probably be ground distance and the projected case should probably be 2D only. I have opened a github ticket for it. Thanks for spotting.

https://github.com/gwaldron/osgearth/issues/1668

//GeoData
/**
* Calculates the distance in meters from this geopoint to another.
* dim = 3, 3d(spatial) distance; dim = 2, 2d(ground) distance
*/
double distanceTo(const GeoPoint& rhs, int dim = 3) const;

//GeoData.cpp
double
GeoPoint::distanceTo(const GeoPoint& rhs, int dim) const
{
    if ( getSRS()->isProjected() && rhs.getSRS()->isProjected() )
    {
        if ( getSRS()->isEquivalentTo(rhs.getSRS()) )
        {
            osg::Vec3d vec = vec3d() - rhs.vec3d();
            if (dim == 2)//ground distance
                vec.z() = 0;
            return vec.length();
        }
        else
        {
            GeoPoint rhsT = rhs.transform(getSRS());
            osg::Vec3d vec = vec3d() - rhsT.vec3d();
           (dim == 2)//ground distance
                vec.z() = 0;
            return vec.length();
        }
    }
    else
    {
        // https://en.wikipedia.org/wiki/Geographical_distance#Ellipsoidal-surface_formulae

        GeoPoint p1 = transform( getSRS()->getGeographicSRS() );
        GeoPoint p2 = rhs.transform( p1.getSRS() );

        double Re = getSRS()->getEllipsoid()->getRadiusEquator();
        double Rp = getSRS()->getEllipsoid()->getRadiusPolar();
        double F  = (Re-Rp)/Re; // flattening

        double
            lat1 = osg::DegreesToRadians(p1.y()),
            lon1 = osg::DegreesToRadians(p1.x()),
            lat2 = osg::DegreesToRadians(p2.y()),
            lon2 = osg::DegreesToRadians(p2.x());

        double B1 = atan( (1.0-F)*tan(lat1) );
        double B2 = atan( (1.0-F)*tan(lat2) );

        double P = (B1+B2)/2.0;
        double Q = (B2-B1)/2.0;

        double G = acos(sin(B1)*sin(B2) + cos(B1)*cos(B2)*cos(fabs(lon2-lon1)));

        double
            sinG = sin(G),
            sinP = sin(P), sinQ = sin(Q),
            cosP = cos(P), cosQ = cos(Q),
            sinG2 = sin(G/2.0), cosG2 = cos(G/2.0);

        double X = (G-sinG)*((sinP*sinP*cosQ*cosQ)/(cosG2*cosG2));
        double Y = (G+sinG)*((cosP*cosP*sinQ*sinQ)/(sinG2*sinG2));

        double dist = Re*(G-(F/2.0)*(X+Y));

        // NaN could mean start/end points are the same
        dist = osg::isNaN(dist)? 0.0 : dist;

        if (dim == 3)//spatial distance
        {
            //altitudeMode of p1 and p2 must be ALTMODE_ABSOLUTE
            dist = hypot(dist, p2.z() - p1.z());//refer to proj_lpz_dist in proj4
        }
        return dist;
    }
}

This is my change.

Let's continue this discussion under the GitHub ticket - thanks