如何在Matlab中检测两条WGS84坐标系线段的交点？

WGS84是椭球地理坐标系，直接基于经纬度计算线段交点会引入较大误差，建议先将地理坐标转换为平面投影坐标（如UTM投影），完成平面线段交点检测后，再转换回WGS84经纬度。以下是具体实现方案：

步骤1：将WGS84经纬度转换为平面投影坐标 #

借助Matlab Mapping Toolbox的坐标转换工具，先确定UTM投影带，再完成经纬度到平面坐标的转换：

% 定义两条线段的WGS84经纬度端点
lat0 = [latStart0, latEnd0];
lon0 = [lonStart0, lonEnd0];
lat1 = [latStart1, latEnd1];
lon1 = [lonStart1, lonEnd1];

% 确定UTM投影带（以第一条线段起点的带号为准）
zone = utmzone(latStart0, lonStart0);

% 定义投影参考系
utmProj = projcrs(zone, 'Type', 'UTM', 'Datum', 'WGS84');
wgs84Proj = projcrs(4326); % WGS84的EPSG标准代码

% 转换为UTM平面坐标（x为东向距离，y为北向距离）
[x0, y0] = projfwd(utmProj, lat0, lon0);
[x1, y1] = projfwd(utmProj, lat1, lon1);

步骤2：检测平面线段交点 #

方法一：使用Matlab内置函数（R2021b及以上版本） #

利用intersectLineSegments直接计算线段交点：

% 构造线段端点矩阵（每行对应一个端点的[x,y]）
seg0 = [x0(1) y0(1); x0(2) y0(2)];
seg1 = [x1(1) y1(1); x1(2) y1(2)];

% 计算交点及相交状态
[intersectPt, isIntersect] = intersectLineSegments(seg0, seg1);

if isIntersect
    % 将平面交点转换回WGS84经纬度
    [latIntersect, lonIntersect] = projinv(utmProj, intersectPt(1), intersectPt(2));
    disp(['交点经纬度：', num2str(latIntersect), ', ', num2str(lonIntersect)]);
else
    disp('两条线段无交点');
end

方法二：手动实现线段交点算法（兼容旧版本） #

通过向量叉积判断线段相交性并计算交点：

% 定义线段端点
A = [x0(1), y0(1)]; B = [x0(2), y0(2)];
C = [x1(1), y1(1)]; D = [x1(2), y1(2)];

% 计算向量
AB = B - A;
AC = C - A;
AD = D - A;
CD = D - C;
CA = A - C;
CB = B - C;

% 叉积判断相交性
cross1 = cross(AB, AC);
cross2 = cross(AB, AD);
cross3 = cross(CD, CA);
cross4 = cross(CD, CB);

if sign(cross1) ~= sign(cross2) && sign(cross3) ~= sign(cross4)
    % 计算交点坐标
    t = cross(AC, CD) / cross(AB, CD);
    intersectPt = A + t * AB;
    % 转换回WGS84经纬度
    [latIntersect, lonIntersect] = projinv(utmProj, intersectPt(1), intersectPt(2));
    disp(['交点经纬度：', num2str(latIntersect), ', ', num2str(lonIntersect)]);
else
    % 检查端点重合情况
    if isequal(A,C) || isequal(A,D) || isequal(B,C) || isequal(B,D)
        if isequal(A,C)
            disp(['交点为端点：', num2str(lat0(1)), ', ', num2str(lon0(1))]);
        elseif isequal(A,D)
            disp(['交点为端点：', num2str(lat0(1)), ', ', num2str(lon0(1))]);
        elseif isequal(B,C)
            disp(['交点为端点：', num2str(lat0(2)), ', ', num2str(lon0(2))]);
        else
            disp(['交点为端点：', num2str(lat0(2)), ', ', num2str(lon0(2))]);
        end
    else
        disp('两条线段无交点');
    end
end

注意事项 #

• 若两条线段跨UTM带，建议选择覆盖整个区域的自定义平面投影（如兰伯特投影），避免投影误差。
• 仅在小范围区域（几公里内）可近似将经纬度视为平面坐标直接计算，大范围场景必须使用投影转换保证精度。

内容的提问来源于stack exchange，提问作者Kitiara