使用 OpenLayers 创建网格周边

Question

给定 2 个坐标和一个宽度，我在 Openlayer 中创建网格的周长，并在后台使用 Google 地图。 2个坐标之间的距离是500米。 网格宽度必须为250米。 我能够计算出网格周长的另外 2 个点，但是如果我与在 QGIS 等其他软件上创建的周长进行比较，我会发现大约 1 或 2 米的误差，不是在 250 米的长度上，而是在方向。 （角度不正确，周长角应该在ST2点上）。

为了计算另外两个点，我尝试了两种方法。

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. 使用 Vicenty 公式和方位角 + 90° 计算点位置。

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

这两种方法给我的结果都与 QGIS 中的结果不同。 你能帮我理解为什么不同吗？ 谢谢！

Answer 1

在 EPSG:3857 投影中，在 500 米的距离上，即使使用简单的几何图形，任何误差也应远小于 1%。然而，笛卡尔距离仅等于赤道处的实际距离（您可以在 https://openlayers.org/en/v3.20.1/examples/measure.html 中看到差异）。

在方法 1 中，高度计算是笛卡尔距离，并假设宽度采用相同的单位。

您当前的代码将生成一个包含点 4 和 3 的矩形，其距离点 1 和 2 的“宽度”笛卡尔单位。在纬度 60 处，实际距离仅为该距离的一半。

要获得正确的笛卡尔宽度以用于实际宽度，您必须调整纬度：

使用 OpenLayers 点解析方法（使用中点以获得最佳精度）

width = width / getPointResolution('EPSG:3857', 1, [(point1[0] + point2[0]) / 2, (point1[1] + point2[1]) / 2]);

或者直接使用纬度的余弦值

width = width / Math.cos(toLonLat([(point1[0] + point2[0]) / 2, (point1[1] + point2[1]) / 2])[1] * Math.PI / 180);