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.
-
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];}
-
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!