Convert Easting/Northing to Lat Lng

Question

I've been struggling to convert northing/eastings to lat/long without using a nuget package. I'm sure my UTM zone is 30U (Great Britain).

How would one go about doing this in C#?

When manually done using this site. It provides the correct values and locations. Easting: 426342 Northing: 505339 Lat/Long should be: 54.44277977022131, -1.5953328509040021

I had previously found some code on stackoverflow but it's giving me the wrong values.

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

Refer the below code and run 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

I ran into the same issue. Ended up using this library: https://github.com/Tronald/CoordinateSharp

This sample will return the exact position in latlon decimal.

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