使用boost几何从点到线的垂直地理距离

这个答案可以做所有的计算, 无助力 .

考虑半径r=1的球。

great circle . 这个大圆圈 gcAB 也穿过中心点 O , B , O PL1 .

P 也在一个大圆圈里。

距离的最小距离(沿大圆的弧测量,而不是沿三维直线) 到大圈子去 是弧的长度 .
PL2号机组 gcPC 垂直于平面 PL1号机组 .

C ,它位于 ,这是上述两个平面的交点。

.
PL1号机组 pp1 . 该向量由 向量的数量 OA 和 OB

因为飞机 PL2号机组 垂直于平面 ,它必须包含向量 第1页 . 所以垂直向量 pp2 到平面 OP 第1页

向量 ppi 排队 OC 两个平面的交集由 第1页 和 第2页

正常化 矢量 ppi公司 R 我们得到了点的坐标 .
叉积是不可交换的。这意味着如果我们交换点A,B,我们得到相反的点 C' 在球体中。我们可以测试距离 PC PC' 得到它们的最小值。

计算 大圆距离 Wikipedia link A. B a 字里行间 办公自动化 .
a = atan2(y, x) 其中,使用半径1, y= sin(a) 和 x= cos(a) . sin(a) cos(a) 点积

把所有的东西放在一起,我们有这个C++代码:

#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;
}

对于给定的三个点,AP=21024.4 BP=12952.1和PC=499.493。

运行代码 here

看来你可以用 project_point 策略:

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;
}

印刷品

12.9521
763.713

我没有检查结果(在那张照片中,犯罪嫌疑人的范围如此之大似乎有点令人惊讶),但也许我遗漏了什么。

projected_point 策略:“ ".