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

Question

我想获得从点(t)到线段的垂直距离(p, q)。垂线不能与线相交[p, q]。在这种情况下，我想(p, q)假设地延长线，然后绘制垂直线以获得距离。p、q、t 都是 gps 坐标。我正在使用增强几何。

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

我在地图上标出了这三个位置

我衡量两个距离qt从和距离t来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;

这两个距离是相同的。这意味着它不计算垂直距离。相反，它计算了shortest从一个点到一条线的距离bounds。

如何以这种方式计算垂直距离，使其无论边界如何都必须垂直？

cpp.sh 中的工作示例

Answer 1

这个答案做所有的计算，没有 boost。

考虑一个半径为 R = 1 的球体。

A、B 点在一个大圆上。这个大圆gcAB也穿过球体的中心点O（大圆需要）。点A、B、O定义了一个平面PL1。

P点也位于一个大圆内。

从P到大圆的最小距离（沿大圆的弧测量，而不是沿 3D 直线测量）gcAB是圆弧PC的长度。 大圆的
平面PL2gcPC垂直于平面PL1。

我们想要点C，它位于线OC 上，这是两个提到的平面的交点。

.
平面PL1由其垂直向量 定义pp1。该向量是通过向量和的叉积获得的。OAOB

因为平面PL2垂直于平面PL1，所以它必须包含向量pp1。因此，垂直矢量pp2于平面PL2可以通过的交叉乘积来获得OP和pp1。

一种载体，ppi在该行OC两平面的相交通过的横产物获得pp1和pp2。

如果我们归一化向量ppi并将其分量乘以R地球的半径，我们得到点C的坐标。
叉积不是可交换的。这意味着如果我们交换点 A、B，我们将得到球体中的相反点C'。我们可以测试距离PC和PC'并得到他们的最低水平。

要计算两点A、B的大圆距离 维基百科链接，它依赖于线和之间的角度。 为了在所有角度上获得最佳精度，我们使用where，使用半径 1和。并且可以分别通过叉积（OA，OB）和点积（OA，OB）计算。aOAOB
a = atan2(y, x)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;
}

对于给定的三个点，Wich 给出了距离 AP = 21024.4 BP = 12952.1 和 PC = 499.493。

在这里运行代码