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

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

我已经在 map 上绘制了这三个位置

我测量了两个距离 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

看来您可以使用 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

我没有检查结果（在那张图片中 perp_dist 大得多似乎有点令人惊讶），但也许我遗漏了什么。

如果你必须做一些特殊的事情（除了坐标系）才能得到 "haversine"（很抱歉我没有跟上这个速度），你可能需要通过projected_point 策略的第二个模板参数：“底层点对点距离策略”。

Answer 2

这个答案进行所有计算，没有提升。

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

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

点P也位于大圆内

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

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

。
平面 PL1 由其垂直矢量 pp1 定义。这个向量是通过向量 OA 和 OB.

的叉积获得的

因为平面PL2垂直于平面PL1，它必须包含向量pp1。所以平面PL2的垂直向量pp2可以通过OP和pp1.

的叉积得到

两个平面的交点OC中的向量ppi由pp1和pp2的叉积得到。

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

计算两点的大圆距离 Wikipedia link A, B，它依赖于直线 OA 和 OB 之间的角度 a。
为了在所有角度上获得最佳精度，我们使用 a = atan2(y, x)，其中使用半径 1，y= sin(a) 和 x= cos(a)。 sin(a)和cos(a)可以分别通过叉积(OA,OB)和点积(OA,OB)计算

将所有内容放在一起我们有这个 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

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

perpendicular geo distance from a point to a line using boost geometry

gis

geometry

boost

boost-geometry