OpenLayers创建网格周长与QGIS结果存在方位偏差的原因排查

我在以Google地图为底图的OpenLayers中，基于两个间距500米的坐标和250米宽度创建网格周长，计算出另外两个顶点后，和QGIS生成的结果相比，ST2点存在1-2米的方位偏差。我尝试了两种计算方法，但结果都和QGIS不符：

尝试的计算方法 #

方法1：平面几何90度计算 #

const rectPoints = (point1, point2, width) => {
    const height = Math.sqrt(Math.pow(point1[0] - point2[0], 2) + 
        Math.pow(point1[1] - point2[1], 2));
    const d3y = (point1[0] - point2[0]) * width / height;
    const d3x = (point1[1] - point2[1]) * width / height;
    const point3 = [point2[0] - d3x, point2[1] + d3y];
    const d4x = d3x;
    const d4y = d3y;
    const point4 = [point1[0] - d4x, point1[1] + d4y]
    return [point1, point2, point3, point4];
}

方法2：Vincenty公式结合方位角计算 #

var bearing = calculateBearing(pt1[1], pt1[0], pt2[1], pt2[0]);
var perpendicular_bearing = adjustBearing(bearing, 90);
var result = destVincenty(pt1[1], pt1[0], perpendicular_bearing, 250);

function adjustBearing(bearing, degrees) {
    let result = (bearing + degrees) % 360;
    if (result < 0) {
        result += 360;
    }
    return result;    
}
                    
function destVincenty(lat1, lon1, brng, dist, callback) {
    var a = 6378137, b = 6356752.3142, f = 1 / 298.257223563; 
    var s = dist;
    var alpha1 = toRad(brng);
    var sinAlpha1 = Math.sin(alpha1);
    var cosAlpha1 = Math.cos(alpha1);
            
    var tanU1 = (1 - f) * Math.tan(toRad(lat1));
    var cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
    var sigma1 = Math.atan2(tanU1, cosAlpha1);
    var sinAlpha = cosU1 * sinAlpha1;
    var cosSqAlpha = 1 - sinAlpha * sinAlpha;
    var uSq = cosSqAlpha * (a * a - b * b) / (b * b);
    var A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    var B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
            
    var sigma = s / (b * A), sigmaP = 2 * Math.PI;
    while (Math.abs(sigma - sigmaP) > 1e-12) {
        var cos2SigmaM = Math.cos(2 * sigma1 + sigma);
        var sinSigma = Math.sin(sigma);
        var cosSigma = Math.cos(sigma);
        var deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
            B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
        sigmaP = sigma;
        sigma = s / (b * A) + deltaSigma;
    }
            
    var tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
    var lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
        (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
    var lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
    var C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
    var L = lambda - (1 - C) * f * sinAlpha *
        (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    var lon2 = (toRad(lon1) + L + 3 * Math.PI) % (2 * Math.PI) - Math.PI;  // normalise to -180...+180
            
    var revAz = Math.atan2(sinAlpha, -tmp);  // final bearing, if required
            
    var result = { lat: toDeg(lat2), lon: toDeg(lon2), finalBearing: toDeg(revAz) };
            
    if (callback !== undefined && callback instanceof Function) {
        if (callback.length === 3) {
            callback(result.lat, result.lon, result.finalBearing);
        }
        else {
            callback(result);
        }
    }
            
    return result;
}
    
function calculateBearing(lat1, lon1, lat2, lon2) {
   // Convert latitude and longitude from degrees to radians
   let lat1Rad = toRad(lat1);
   let lon1Rad = toRad(lon1);
   let lat2Rad = toRad(lat2);
   let lon2Rad = toRad(lon2);

    // Difference in the coordinates
    let dLon = lon2Rad - lon1Rad;

    // Calculate bearing
    let x = Math.sin(dLon) * Math.cos(lat2Rad);
    let y = Math.cos(lat1Rad) * Math.sin(lat2Rad) - Math.sin(lat1Rad) * Math.cos(lat2Rad) * Math.cos(dLon);
    let initialBearing = Math.atan2(x, y);

    // Convert bearing from radians to degrees
    initialBearing = initialBearing * 180 / Math.PI;

    // Normalize the bearing
    let bearing = (initialBearing + 360) % 360;

    return bearing;
}

偏差原因分析 #

方法1的核心问题：平面几何忽略投影变形 #

第一种方法采用平面直角坐标系的垂直计算逻辑，但Google地图的底图是Web Mercator（EPSG:3857）投影，这是一种球面投影，会在非赤道区域产生距离和角度的非线性变形。直接用平面几何计算垂直方向，完全没有考虑椭球面的特性，必然和QGIS基于椭球的精准计算出现偏差。

方法2的潜在问题：细节处理不一致 #

第二种方法用Vincenty公式（椭球精准计算），但可能存在以下细节问题：

• 方位角方向错误：Vincenty公式的方位角是正北顺时针计算，你给原方位角加90度得到的是右侧垂直方向，但QGIS可能生成的是左侧垂直方向，方向相反会导致点的位置偏差。
• 坐标顺序混淆：你的calculateBearing和destVincenty接收的参数是(lat, lon)，但OpenLayers的坐标通常是[lon, lat]（EPSG:4326）或Web Mercator的[x, y]，如果输入坐标顺序和QGIS不一致，会导致方位角计算错误。
• 投影基准不统一：QGIS默认用WGS84（EPSG:4326）椭球计算，如果你直接用OpenLayers中的Web Mercator坐标进行Vincenty计算，没有先转换为WGS84经纬度，会产生系统偏差。

修正方案 #

1. 统一坐标基准：将OpenLayers中的Web Mercator坐标转换为WGS84（EPSG:4326）经纬度后，再进行Vincenty计算，确保和QGIS的椭球基准一致。
2. 验证垂直方向：尝试将方位角调整为bearing - 90（或bearing + 270），对比QGIS结果确认正确的垂直方向。
3. 使用OpenLayers内置工具：OpenLayers的ol/sphere模块提供了封装好的球面偏移计算方法，比自己实现Vincenty更可靠：

import { offset } from 'ol/sphere';

// 假设pt1、pt2是EPSG:4326格式的经纬度坐标 [lon, lat]
const bearing = calculateBearing(pt1[1], pt1[0], pt2[1], pt2[0]);
// 计算垂直偏移点，方向根据实际情况用+90或-90调整
const pt4 = offset(pt1, 250, bearing - 90);
const pt3 = offset(pt2, 250, bearing - 90);
// 将结果转换回Web Mercator用于OpenLayers展示

4. 对齐QGIS参数：检查QGIS创建网格时的设置，确认是用“椭球距离”还是“平面距离”，以及垂直方向的定义，确保你的计算逻辑和QGIS完全匹配。

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