Create grid perimeter with OpenLayers

Question

Given 2 coordinates and a width, I'm creating the perimeter of a grid in Openlayer with Google map in background. The distance between the 2 coordinates is 500 meters. And the width of the grid has to be 250 meters. I was able to compute the 2 other points of the grid perimeter, but if I compare with the perimeter created on other software like QGIS, I have an error of about 1 or 2 meters, not in the length of the 250 meters but in the orientation. (the angle is not correct, the perimeter angle should be on the point ST2).

To compute the two other points, I tried 2 ways.

1. Compute the point with simple geometric calculation of a 90° angle

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. Compute the point position with Vicenty's formulae and the bearing + 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;
}

Both ways give me a result different than the one in QGIS. Can you help me understand why it is different? Thanks!

Answer 1

Over a distance of 500 meters in EPSG:3857 projection any errors should be much less than 1 percent even with simple geometry. However cartesian distances will equal the actual distances only at the equator (you can see the difference in https://openlayers.org/en/v3.20.1/examples/measure.html).

In your method 1 your height calculation is cartesian distance and assumes width is in the same units.

Your current code will produce a rectangle with points 4 and 3 at "width" cartesian units from points 1 and 2. At latitude 60 the actual distance would be only half that.

To get the correct cartesian width to use for the actual width you must adjust for the latitude:

either using OpenLayers point resolution method (using a midpoint for best accuracy)

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

or using the cosine of the latitude directly

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