perpendicular geo distance from a point to a line

问题:

I want to get perpendicular distance from a point (t) to a line segment (p, q). The perpendicular may not intersect the line [p, q]. In that case I want to extend the line (p, q) hypothetically and then draw the perpendicular to get the distance. p, q, t are all gps coordinates. I am using boost geometry.

typedef boost::geometry::model::point<
    double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>
> geo_point;
typedef boost::geometry::model::segment<geo_point> geo_segment;

geo_point p(88.41253929999999, 22.560206299999997);
geo_point q(88.36928063300775, 22.620867969497795);
geo_point t(88.29580956367181, 22.71558662052875);

I have plotted these three locations on map

I measure two distances qt and distance from t to pq

double dist_qt = boost::geometry::distance(q, t);
std::cout << dist_qt*earth_radius << std::endl;

geo_segment line(p, q);
double perp_dist = boost::geometry::distance(t, line);
std::cout << perp_dist*earth_radius << std::endl;

Both of these distances are same. This means it doesn't calculate the perpendicular distance. Rather it calculated shortest distance from a point to a line within bounds.

How can I calculate the perpendicular distance in such a way so that it is necessarily perpendicular irrespective of bounds ?

Working Example in cpp.sh

回答1:

This answer do all calculations, without boost.

---

Consider a sphere of radius R = 1.

Points A, B are on a great circle. This great circle gcAB goes also through the center point O of the sphere (required for great circles). Points A, B, O define a plane PL1.

Point P also lies in a great circle.

The minimum distance (measured along the arc of a great circle, not along a 3D straight line) from P to the great circle gcAB is the length of the arc PC.
Plane PL2 of great circle gcPC is perpendicular to the plane PL1.

We want point C, which lies in the line OC, which is the intersection of the two mentioned planes.

.
Plane PL1 is defined by its perpendicular vector pp1. This vector is obtained by the cross product of vectors OA and OB.

Because the plane PL2 is perpendicular to plane PL1, it must contain the vector pp1. So the perpendicular vector pp2 to plane PL2 can be obtained by the cross product of OP and pp1.

A vector ppi in the line OC intersection of both planes is obtained by the cross product of pp1 and pp2.

If we normalize vector ppi and multiply its components by the radius R of Earth, we get the coordinates of point C.
Cross product is not commutative. This means that if we interchange points A,B we get the opposite point C' in the sphere. We can test distances PC and PC' and get their minimum.

---

To calculate the great-circle distance Wikipedia link for two points A, B, it relies on the angle a between the lines OA and OB.
For best accuracy on all angles we use a = atan2(y, x) where, using radius 1, y= sin(a) and x= cos(a). sin(a) and cos(a) can be calculated by cross product (OA, OB) and dot product (OA, OB) respectively.

Putting all together we have this C++ code:

#include <iostream>
#include <cmath>

const double degToRad = std::acos(-1) / 180;

struct vec3
{
    double x, y, z;
    vec3(double xd, double yd, double zd) : x(xd), y(yd), z(zd) {}
    double length()
    {
        return std::sqrt(x*x + y*y + z*z);
    }
    void normalize()
    {
        double len = length();
        x = x / len;
        y = y / len;
        z = z / len;
    } 
};

vec3 cross(const vec3& v1, const vec3& v2)
{
    return vec3( v1.y * v2.z - v2.y * v1.z,
                 v1.z * v2.x - v2.z * v1.x,
                 v1.x * v2.y - v2.x * v1.y );
}

double dot(const vec3& v1, const vec3& v2)
{
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

double GCDistance(const vec3& v1, const vec3& v2, double R)
{
    //normalize, so we can pass any vectors
    vec3 v1n = v1;
    v1n.normalize();
    vec3 v2n = v2;
    v2n.normalize();
    vec3 tmp = cross(v1n, v2n);
    //minimum distance may be in one direction or the other
    double d1 = std::abs(R * std::atan2(tmp.length() , dot(v1n, v2n)));
    double d2 = std::abs(R * std::atan2(tmp.length() , -dot(v1n, v2n)));

    return std::min(std::abs(d1), std::abs(d2));
}                  

int main()
{
    //Points A, B, and P
    double lon1 = 88.41253929999999  * degToRad;
    double lat1 = 22.560206299999997 * degToRad;
    double lon2 = 88.36928063300775  * degToRad;
    double lat2 = 22.620867969497795 * degToRad;
    double lon3 = 88.29580956367181  * degToRad;
    double lat3 = 22.71558662052875  * degToRad;

    //Let's work with a sphere of R = 1
    vec3 OA(std::cos(lat1) * std::cos(lon1), std::cos(lat1) * std::sin(lon1), std::sin(lat1));
    vec3 OB(std::cos(lat2) * std::cos(lon2), std::cos(lat2) * std::sin(lon2), std::sin(lat2));
    vec3 OP(std::cos(lat3) * std::cos(lon3), std::cos(lat3) * std::sin(lon3), std::sin(lat3));
    //plane OAB, defined by its perpendicular vector pp1
    vec3 pp1 = cross(OA, OB);
    //plane OPC
    vec3 pp2 = cross(pp1, OP);
    //planes intersection, defined by a line whose vector is ppi
    vec3 ppi = cross(pp1, pp2);
    ppi.normalize(); //unitary vector

    //Radious or Earth
    double R = 6371000; //mean value. For more precision, data from a reference ellipsoid is required

    std::cout << "Distance AP = " << GCDistance(OA, OP, R) << std::endl;
    std::cout << "Distance BP = " << GCDistance(OB, OP, R) << std::endl;
    std::cout << "Perpendicular distance (on arc) = " << GCDistance(OP, ppi, R) << std::endl;
}

Wich gives the distances AP = 21024.4 BP = 12952.1 and PC= 499.493 for the given three points.

Running code here

回答2:

Seems like you could use the project_point strategy:

Live On Coliru

#include <string>
#include <iostream>
#include <boost/geometry.hpp>

namespace bg = boost::geometry;

int main(){
    double const earth_radius = 6371.0; // Km

    typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree>> geo_point;
    typedef bg::model::segment<geo_point> geo_segment;

    geo_point p(88.41253929999999, 22.560206299999997);
    geo_point q(88.36928063300775, 22.620867969497795);
    geo_point t(88.29580956367181, 22.71558662052875);

    double dist_qt = bg::distance(q, t);
    std::cout << dist_qt*earth_radius << std::endl;

    geo_segment line(p, q);
    double perp_dist = distance(t, line, bg::strategy::distance::projected_point<>{});
    std::cout << perp_dist*earth_radius << std::endl;
}

Prints

12.9521
763.713

I didn't check the results (in that picture it seems slightly surprising that the perp_dist be so much larger), but perhaps I'm missing something.

In case you have to do special things (apart from the coordinate system) to get that "haversine" (I'm sorry I'm not up to speed with this), you might need to pass the second template argument to the projected_point strategy: "underlying point-point distance strategy".