GeoPoint distanceTo

classic Classic list List threaded Threaded
4 messages Options
liugaijin liugaijin
Reply | Threaded
Open this post in threaded view
|

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?

gwaldron gwaldron
Reply | Threaded
Open this post in threaded view
|

Re: GeoPoint distanceTo

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
Glenn Waldron / Pelican Mapping
liugaijin liugaijin
Reply | Threaded
Open this post in threaded view
|

Re: GeoPoint distanceTo

In reply to this post by liugaijin
//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.
gwaldron gwaldron
Reply | Threaded
Open this post in threaded view
|

Re: GeoPoint distanceTo

Let's continue this discussion under the GitHub ticket - thanks
Glenn Waldron / Pelican Mapping