将东/北转换为纬度 Lng

Question

我一直在努力在不使用 nuget 包的情况下将北/东转换为纬度/经度。我确定我的 UTM 区域是 30U（英国）。

如何在 C# 中执行此操作？

使用this 站点手动完成时。它提供了正确的值和位置。 东行：426342 北行：505339 纬度/经度应为：54.44277977022131，-1.5953328509040021

我之前在 stackoverflow 上有 found 一些代码，但它给了我错误的值。

utmZone = "30U"

  public static LatLng ToLatLon(double utmX, double utmY, string utmZone)
    {
    bool isNorthHemisphere = utmZone.Last() >= 'N';

    var diflat = -0.00066286966871111111111111111111111111;
    var diflon = -0.0003868060578;

    var zone = int.Parse(utmZone.Remove(utmZone.Length - 1));
    var c_sa = 6378137.000000;
    var c_sb = 6356752.314245;
    var e2 = Math.Pow((Math.Pow(c_sa, 2) - Math.Pow(c_sb, 2)), 0.5) / c_sb;
    var e2cuadrada = Math.Pow(e2, 2);
    var c = Math.Pow(c_sa, 2) / c_sb;
    var x = utmX - 500000;
    var y = isNorthHemisphere ? utmY : utmY - 10000000;

    var s = ((zone * 6.0) - 183.0);
    var lat = y / (c_sa * 0.9996);
    var v = (c / Math.Pow(1 + (e2cuadrada * Math.Pow(Math.Cos(lat), 2)), 0.5)) * 0.9996;
    var a = x / v;
    var a1 = Math.Sin(2 * lat);
    var a2 = a1 * Math.Pow((Math.Cos(lat)), 2);
    var j2 = lat + (a1 / 2.0);
    var j4 = ((3 * j2) + a2) / 4.0;
    var j6 = ((5 * j4) + Math.Pow(a2 * (Math.Cos(lat)), 2)) / 3.0;
    var alfa = (3.0 / 4.0) * e2cuadrada;
    var beta = (5.0 / 3.0) * Math.Pow(alfa, 2);
    var gama = (35.0 / 27.0) * Math.Pow(alfa, 3);
    var bm = 0.9996 * c * (lat - alfa * j2 + beta * j4 - gama * j6);
    var b = (y - bm) / v;
    var epsi = ((e2cuadrada * Math.Pow(a, 2)) / 2.0) * Math.Pow((Math.Cos(lat)), 2);
    var eps = a * (1 - (epsi / 3.0));
    var nab = (b * (1 - epsi)) + lat;
    var senoheps = (Math.Exp(eps) - Math.Exp(-eps)) / 2.0;
    var delt = Math.Atan(senoheps / (Math.Cos(nab)));
    var tao = Math.Atan(Math.Cos(delt) * Math.Tan(nab));

    double longitude = ((delt * (180.0 / Math.PI)) + s) + diflon;
    double latitude = ((lat + (1 + e2cuadrada * Math.Pow(Math.Cos(lat), 2) - (3.0 / 2.0) * e2cuadrada * Math.Sin(lat) * Math.Cos(lat) * (tao - lat)) * (tao - lat)) * (180.0 / Math.PI)) + diflat;

    return new LatLng
    {
        Latitude = latitude,
        Longitude = longitude
    };
}

Answer 1

参考以下代码并运行LatLonConversions.ConvertOSToLatLon(426342, 505339)

public class LatLonConversions
{
	const double a = 6377563.396;
	const double b = 6356256.91;
	const double e2 = (a - b) / a;
	const double n0 = -100000;
	const double e0 = 400000;
	const double f0 = 0.999601272;
	const double phi0 = 0.855211333;
	const double lambda0 = -0.034906585;
	const double n = (a - b) / (a + b);

	static double lat, lng;

	private LatLonConversions() { }

	private static double Deg2Rad(double x)
	{
		return x * (Math.PI / 180);
	}

	private static double Rad2Deg(double x)
	{
		return x * (180 / Math.PI);
	}

	private static double SinSquared(double x)
	{
		return Math.Sin(x) * Math.Sin(x);
	}

	private static double TanSquared(double x)
	{
		return Math.Tan(x) * Math.Tan(x);
	}

	private static double Sec(double x)
	{
		return 1.0 / Math.Cos(x);
	}

	private static void OSGB36ToWGS84()
	{
		var airy1830 = new RefEll(6377563.396, 6356256.909);
		var a = airy1830.maj;
		var b = airy1830.min;
		var eSquared = airy1830.ecc;
		var phi = Deg2Rad(lat);
		var lambda = Deg2Rad(lng);
		var v = a / (Math.Sqrt(1 - eSquared * SinSquared(phi)));
		var H = 0; // height
		var x = (v + H) * Math.Cos(phi) * Math.Cos(lambda);
		var y = (v + H) * Math.Cos(phi) * Math.Sin(lambda);
		var z = ((1 - eSquared) * v + H) * Math.Sin(phi);

		var tx = 446.448;
		var ty = -124.157;
		var tz = 542.060;
		var s = -0.0000204894;
		var rx = Deg2Rad(0.00004172222);
		var ry = Deg2Rad(0.00006861111);
		var rz = Deg2Rad(0.00023391666);

		var xB = tx + (x * (1 + s)) + (-rx * y) + (ry * z);
		var yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
		var zB = tz + (-ry * x) + (rx * y) + (z * (1 + s));

		var wgs84 = new RefEll(6378137.000, 6356752.3141);
		a = wgs84.maj;
		b = wgs84.min;
		eSquared = wgs84.ecc;

		var lambdaB = Rad2Deg(Math.Atan(yB / xB));
		var p = Math.Sqrt((xB * xB) + (yB * yB));
		var phiN = Math.Atan(zB / (p * (1 - eSquared)));
		for (var i = 1; i < 10; i++)
		{
			v = a / (Math.Sqrt(1 - eSquared * SinSquared(phiN)));
			double phiN1 = Math.Atan((zB + (eSquared * v * Math.Sin(phiN))) / p);
			phiN = phiN1;
		}

		var phiB = Rad2Deg(phiN);

		lat = phiB;
		lng = lambdaB;
	}

	public static LatLon ConvertOSToLatLon(double easting, double northing)
	{
		RefEll airy1830 = new RefEll(6377563.396, 6356256.909);
		double OSGB_F0 = 0.9996012717;
		double N0 = -100000.0;
		double E0 = 400000.0;
		double phi0 = Deg2Rad(49.0);
		double lambda0 = Deg2Rad(-2.0);
		double a = airy1830.maj;
		double b = airy1830.min;
		double eSquared = airy1830.ecc;
		double phi = 0.0;
		double lambda = 0.0;
		double E = easting;
		double N = northing;
		double n = (a - b) / (a + b);
		double M = 0.0;
		double phiPrime = ((N - N0) / (a * OSGB_F0)) + phi0;
		do
		{
			M =
			  (b * OSGB_F0)
				* (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n))
				  * (phiPrime - phi0))
				  - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n))
					* Math.Sin(phiPrime - phi0)
					* Math.Cos(phiPrime + phi0))
				  + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n))
					* Math.Sin(2.0 * (phiPrime - phi0))
					* Math.Cos(2.0 * (phiPrime + phi0)))
				  - (((35.0 / 24.0) * n * n * n)
					* Math.Sin(3.0 * (phiPrime - phi0))
					* Math.Cos(3.0 * (phiPrime + phi0))));
			phiPrime += (N - N0 - M) / (a * OSGB_F0);
		} while ((N - N0 - M) >= 0.001);
		var v = a * OSGB_F0 * Math.Pow(1.0 - eSquared * SinSquared(phiPrime), -0.5);
		var rho =
		  a
			* OSGB_F0
			* (1.0 - eSquared)
			* Math.Pow(1.0 - eSquared * SinSquared(phiPrime), -1.5);
		var etaSquared = (v / rho) - 1.0;
		var VII = Math.Tan(phiPrime) / (2 * rho * v);
		var VIII =
		  (Math.Tan(phiPrime) / (24.0 * rho * Math.Pow(v, 3.0)))
			* (5.0
			  + (3.0 * TanSquared(phiPrime))
			  + etaSquared
			  - (9.0 * TanSquared(phiPrime) * etaSquared));
		var IX =
		  (Math.Tan(phiPrime) / (720.0 * rho * Math.Pow(v, 5.0)))
			* (61.0
			  + (90.0 * TanSquared(phiPrime))
			  + (45.0 * TanSquared(phiPrime) * TanSquared(phiPrime)));
		var X = Sec(phiPrime) / v;
		var XI =
		  (Sec(phiPrime) / (6.0 * v * v * v))
			* ((v / rho) + (2 * TanSquared(phiPrime)));
		var XII =
		  (Sec(phiPrime) / (120.0 * Math.Pow(v, 5.0)))
			* (5.0
			  + (28.0 * TanSquared(phiPrime))
			  + (24.0 * TanSquared(phiPrime) * TanSquared(phiPrime)));
		var XIIA =
		  (Sec(phiPrime) / (5040.0 * Math.Pow(v, 7.0)))
			* (61.0
			  + (662.0 * TanSquared(phiPrime))
			  + (1320.0 * TanSquared(phiPrime) * TanSquared(phiPrime))
			  + (720.0
				* TanSquared(phiPrime)
				* TanSquared(phiPrime)
				* TanSquared(phiPrime)));
		phi =
		  phiPrime
			- (VII * Math.Pow(E - E0, 2.0))
			+ (VIII * Math.Pow(E - E0, 4.0))
			- (IX * Math.Pow(E - E0, 6.0));
		lambda =
		  lambda0
			+ (X * (E - E0))
			- (XI * Math.Pow(E - E0, 3.0))
			+ (XII * Math.Pow(E - E0, 5.0))
			- (XIIA * Math.Pow(E - E0, 7.0));


		lat = Rad2Deg(phi);
		lng = Rad2Deg(lambda);
		// convert to WGS84
		OSGB36ToWGS84();

		return new LatLon(lat, lng);
	}
}

public class RefEll
{
	public double maj, min, ecc;
	public RefEll(double major, double minor)
	{
		maj = major;
		min = minor;
		ecc = ((major * major) - (minor * minor)) / (major * major);
	}
}

public class LatLon
{
	public double Latitude;
	public double Longitude;

	public LatLon()
	{
		Latitude = 0;
		Longitude = 0;
	}

	public LatLon(double lat, double lon)
	{
		Latitude = lat;
		Longitude = lon;
	}
}

Answer 2

我遇到了同样的问题。最终使用了这个库：https://github.com/Tronald/CoordinateSharp

此示例将返回精确的纬度十进制位置。

UniversalTransverseMercator utm = new UniversalTransverseMercator("T", 30, 581177.3879, 4794824.5279);
Coordinate c = UniversalTransverseMercator.ConvertUTMtoLatLong(utm);
c.FormatOptions.Format = CoordinateFormatType.Decimal;
c.FormatOptions.Round = 7;

Debug.Log($"({c.Latitude}, {c.Longitude})");