东移和北移到经纬度问题

Question

我在将以下代码转换为经纬度时遇到问题。

下面的代码总是偏离 50 米，我看不出有什么问题：

例如，如果您输入东距和北距 517835.99 和 428747.39，它应该转换为经纬度 53.7418 和 -0.21484，但它会给出 53.74159 和 -0.21312。

// Define the constants for the BNG coordinate system
var a = 6377563.396; // Semi-major axis of the Airy 1830 ellipsoid
var b = 6356256.909; // Semi-minor axis of the Airy 1830 ellipsoid
var e0 = 400000; // False easting
var n0 = -100000; // False northing
var F0 = 0.9996012717; // Central meridian scale factor
var phi0 = 49 * (Math.PI / 180); // Central meridian (latitude of true origin) in radians
var lambda0 = -2 * (Math.PI / 180); // Central meridian (longitude of true origin) in radians

// Define the easting and northing coordinates
var easting = 517835.99; // Replace with your easting value
var northing = 428747.39; // Replace with your northing value

// Calculate the transformation variables
var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2); // Eccentricity squared
var n = (a - b) / (a + b);
var n2 = Math.pow(n, 2);
var n3 = Math.pow(n, 3);

var phi = phi0;
var M = 0;

do {
  phi = (northing - n0 - M) / (a * F0) + phi;

  var Ma = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (phi - phi0);
  var Mb = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.sin(phi - phi0) * Math.cos(phi + phi0);
  var Mc = ((15 / 8) * n2 + (15 / 8) * n3) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0));
  var Md = (35 / 24) * n3 * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0));
  M = b * F0 * (Ma - Mb + Mc - Md);
} while (northing - n0 - M >= 0.00001);

var cosPhi = Math.cos(phi);
var sinPhi = Math.sin(phi);
var tanPhi = Math.tan(phi);
var secPhi = 1 / cosPhi;
var nu = a * F0 / Math.sqrt(1 - e2 * Math.pow(sinPhi, 2));
var rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * Math.pow(sinPhi, 2), 1.5);
var eta2 = nu / rho - 1;

var VII = tanPhi / (2 * rho * nu);
var VIII = tanPhi / (24 * rho * Math.pow(nu, 3)) * (5 + 3 * Math.pow(tanPhi, 2) + eta2 - 9 * Math.pow(tanPhi, 2) * eta2);
var IX = tanPhi / (720 * rho * Math.pow(nu, 5)) * (61 + 90 * Math.pow(tanPhi, 2) + 45 * Math.pow(tanPhi, 4));
var X = secPhi / nu;
var XI = secPhi / (6 * Math.pow(nu, 3)) * (nu / rho + 2 * Math.pow(tanPhi, 2));
var XII = secPhi / (120 * Math.pow(nu, 5)) * (5 + 28 * Math.pow(tanPhi, 2) + 24 * Math.pow(tanPhi, 4));
var XIIA = secPhi / (5040 * Math.pow(nu, 7)) * (61 + 662 * Math.pow(tanPhi, 2) + 1320 * Math.pow(tanPhi, 4) + 720 * Math.pow(tanPhi, 6));

var dE = easting - e0;

var lat = phi - VII * Math.pow(dE, 2) + VIII * Math.pow(dE, 4) - IX * Math.pow(dE, 6);
var lon = lambda0 + X * dE - XI * Math.pow(dE, 3) + XII * Math.pow(dE, 5) - XIIA * Math.pow(dE, 7);

// Convert latitude and longitude to degrees
var latitude = lat * (180 / Math.PI);
var longitude = lon * (180 / Math.PI);
// The single point is considered as the centroid
var centroid = [latitude, longitude];


// Output the results
console.log("Latitude: " + latitude);
console.log("Longitude: " + longitude);
// Output the centroid coordinates
console.log("Centroid: " + centroid);  

Answer 1

解决方案适用于赫尔河畔金斯顿：

我已经采用了上面的代码，因为我需要它的原始数学形式或jsmaths。我只是简单地推了等值。我做了很多研究，但没有结果，这给了我所需的值，与实际值相差最多 1.6 米。

public class Main {

    public static void main(String[] args) {
        // Define the constants for the BNG coordinate system
        double a = 6377563.396; // Semi-major axis of the Airy 1830 ellipsoid
        double b = 6356256.909; // Semi-minor axis of the Airy 1830 ellipsoid
        double e0 = 400000; // False easting
        double n0 = -100000; // False northing
        double F0 = 0.9996012717; // Central meridian scale factor
        double phi0 = 49 * (Math.PI / 180); // Central meridian (latitude of true origin) in radians
        double lambda0 = -2 * (Math.PI / 180); // Central meridian (longitude of true origin) in radians

        // Define the easting and northing coordinates
        double easting = 510124.08; // Replace with your easting value
        double northing = 436343.05; // Replace with your northing value

        // Calculate the transformation variables
        double e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2); // Eccentricity squared
        double n = (a - b) / (a + b);
        double n2 = Math.pow(n, 2);
        double n3 = Math.pow(n, 3);

        double phi = phi0;
        double M = 0;

        do {
            phi = (northing - n0 - M) / (a * F0) + phi;

            double Ma = (1 + n + (5 / 4.0) * n2 + (5 / 4.0) * n3) * (phi - phi0);
            double Mb = (3 * n + 3 * n * n + (21 / 8.0) * n3) * Math.sin(phi - phi0) * Math.cos(phi + phi0);
            double Mc = ((15 / 8.0) * n2 + (15 / 8.0) * n3) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0));
            double Md = (35 / 24.0) * n3 * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0));
            M = b * F0 * (Ma - Mb + Mc - Md);
        } while (Math.abs(northing - n0 - M) >= 0.00001); // Use absolute value for convergence check

        double cosPhi = Math.cos(phi);
        double sinPhi = Math.sin(phi);
        double tanPhi = Math.tan(phi);
        double secPhi = 1 / cosPhi;
        double nu = a * F0 / Math.sqrt(1 - e2 * Math.pow(sinPhi, 2));
        double rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * Math.pow(sinPhi, 2), 1.5);
        double eta2 = nu / rho - 1;

        double VII = tanPhi / (2 * rho * nu);
        double VIII = tanPhi / (24 * rho * Math.pow(nu, 3)) * (5 + 3 * Math.pow(tanPhi, 2) + eta2 - 9 * Math.pow(tanPhi, 2) * eta2);
        double IX = tanPhi / (720 * rho * Math.pow(nu, 5)) * (61 + 90 * Math.pow(tanPhi, 2) + 45 * Math.pow(tanPhi, 4));
        double X = secPhi / nu;
        double XI = secPhi / (6 * Math.pow(nu, 3)) * (nu / rho + 2 * Math.pow(tanPhi, 2));
        double XII = secPhi / (120 * Math.pow(nu, 5)) * (5 + 28 * Math.pow(tanPhi, 2) + 24 * Math.pow(tanPhi, 4));
        double XIIA = secPhi / (5040 * Math.pow(nu, 7)) * (61 + 662 * Math.pow(tanPhi, 2) + 1320 * Math.pow(tanPhi, 4) + 720 * Math.pow(tanPhi, 6));

        double dE = easting - e0;

        double lat = phi - VII * Math.pow(dE, 2) + VIII * Math.pow(dE, 4) - IX * Math.pow(dE, 6);
        double lon = lambda0 + X * dE - XI * Math.pow(dE, 3) + XII * Math.pow(dE, 5) - XIIA * Math.pow(dE, 7);

        // Convert latitude and longitude to degrees
        double latitude = lat * (180 / Math.PI);
        double longitude = lon * (180 / Math.PI);

        // Add the corrections 1.6 metres is down from 24 metres inaccuracy
        latitude += 0.000239504;
        longitude -= 0.001692888;

        // The single point is considered as the centroid
        double[] centroid = { latitude, longitude };

        // Output the results
        System.out.println("Latitude: " + latitude);
        System.out.println("Longitude: " + longitude);
        System.out.println("Centroid: " + java.util.Arrays.toString(centroid));
    }
}