Lambert72에서 LatLong으로 C#

Question

Lambet72 좌표를 위도/경도로 변환하는 알고리즘을 찾고 있습니다. 예를 들어, Lambert72 좌표 148990,169450은 50.83546802746704, 4.354415218851164로 변환해야하지만 알고리즘의 대부분에는 약간의 오프셋이 있습니다 (첨부 파일 참조).

Image 2

이 가까운 내가 찾은 알고리즘 중 하나입니다

Image 1

,하지만 여전히 오류가 있습니다.

누군가 더 좋은 알고리즘을 가지고 있습니까?

static void Lambert72ToLatLong(double x, double y, ref double longitude, ref double latitude) 
    { 
     const double n = 0.77164219; 
     const double F = 1.81329763; 
     const double thetaFudge = 0.00014204; 
     const double e = 0.08199189; 
     const double a = 6378388; 
     const double xDiff = 150000; 
     const double yDiff = 5400088.44; 
     const double theta0 = 0.07604294; 

     double xReal = xDiff-x; 
     double yReal = yDiff-y; 

     double rho = Math.Sqrt(xReal*xReal + yReal * yReal); 
     double theta = Math.Atan(xReal/-yReal); 

     longitude = (theta0 + (theta + thetaFudge)/n) * 180/Math.PI; 
     latitude = 0; 
     for(int i=0; i<10; ++i) 
     { 
      latitude = (2 * Math.Atan(Math.Pow(F * a/rho, 1/n) * Math.Pow((1 + e * Math.Sin(latitude))/(1 - e * Math.Sin(latitude)), e/2))) - Math.PI/2; 
     } 
     latitude *= 180/Math.PI; 
    } 

Answer 1

Lambert72 (Belgium Dattum 72 LCC 3p)를 WGS84 GPS 좌표로 변환하려고 시도한 것 같습니다. 따라서 구면 위도로 변환하고 경도 좌표로 변환 한 후에는 의 WGS84 좌표로 끝나기 위해 추가 단계를 수행해야합니다. 은 인용 된 이미지에 표시되는 좌표를 쉽게 확인할 수 있기 때문에 C#의 대체 코드는 Lambert72에서 WGS74까지의 전체 변환 을 사용하여 필요한 소수 정확성을 준수합니다.

static void Lambert72toWGS84latlong(double X, double Y) 
{ 
double LongRef = 0.076042943;  // 
double nLamb = 0.7716421928; 
double aCarre = Math.Pow(6378388.0,2.0); 
double bLamb = 6378388.0 * (1.0 - (1.0/297.0)); 
double eCarre = (aCarre - Math.Pow(bLamb, 2.0))/aCarre; 
double KLamb = 11565915.812935; 

double eLamb = Math.Sqrt(eCarre); 
double eSur2 = eLamb/2.0; 

double Tan1 = (X - 150000.012)/(5400088.437 - Y); 
double Lambda = LongRef + (1.0/nLamb) * (0.000142043 + Math.Atan(Tan1)); 
double RLamb = Math.Sqrt(Math.Pow((X - 150000.012) , 2.0) + Math.Pow ((5400088.437 - Y) ,2.0)); 

double TanZDemi = Math.Pow((RLamb/KLamb),(1.0/nLamb)); 
double Lati1 = 2.0 * Math.Atan(TanZDemi); 

double eSin; 
double Mult1, Mult2, Mult; 
double LatiN, Diff; 

double lat, lng ; 
int i=0; 
do 
    { 
    eSin = eLamb * Math.Sin(Lati1); 
    Mult1 = 1.0 - eSin; 
    Mult2 = 1.0 + eSin; 
    Mult = Math.Pow((Mult1/Mult2) , (eLamb/2.0)); 
    LatiN = (Math.PI/2.0) - (2.0 * (Math.Atan(TanZDemi * Mult))); 
    Diff = LatiN - Lati1; 
    Lati1 = LatiN; 
    i++; 
    } while (Math.Abs(Diff)> 0.0000000277777); 


    lat=LatiN; 
    lng=Lambda; 

    double SinLat = Math.Sin(lat); 
    double SinLng = Math.Sin(lng); 
    double CoSinLat = Math.Cos(lat); 
    double CoSinLng = Math.Cos(lng); 

    double dx = -125.8; 
    double dy = 79.9; 
    double dz = -100.5; 
    double da = -251.0; 
    double df = -0.000014192702; 

    double LWf = 1.0/297.0; 
    double LWa = 6378388.0; 
    double LWb = (1 - LWf) * LWa; 
    double LWe2 = (2.0 * LWf) - (LWf * LWf); 
    double Adb = 1.0/(1.0 - LWf); 

    double Rn = LWa/Math.Sqrt(1.0 - LWe2 * SinLat * SinLat);  
    double Rm = LWa * (1 - LWe2) /Math.Pow((1.0 - LWe2 * lat * lat) ,1.5); 
    double DLat = -dx * SinLat * CoSinLng - dy * SinLat * SinLng + dz * CoSinLat; 
    DLat = DLat + da * (Rn * LWe2 * SinLat * CoSinLat)/LWa; 
    DLat = DLat + df * (Rm * Adb + Rn/Adb) * SinLat * CoSinLat; 
    DLat = DLat/(Rm + 0.0); 

    double DLng = (-dx * SinLng + dy * CoSinLng)/((Rn + 0.0) * CoSinLat); 
    double Dh = dx * CoSinLat * CoSinLng + dy * CoSinLat * SinLng + dz * SinLat; 
    Dh = Dh - da * LWa/Rn + df * Rn * lat * lat/Adb; 

    double LatWGS84 = ((lat + DLat) * 180.0)/Math.PI; 
    double LngWGS84 = ((lng + DLng) * 180.0)/Math.PI; 

    MessageBox.Show("WGS84-Latitude=" + LatWGS84.ToString("###.######") + 
    "--WGS84 Longitude=" + LngWGS84.ToString("###.######")); 

    } 

희망 사항.

Answer 2

당신은 올바른 값은 당신이

를 해결하는 것이이

 xDiff=149910; 
     yDiff=5400150; 

해야

  xDiff 

및

 yDiff 

의 값을 해결하는 데

 i<10 

위도의 계산을위한 용 루프 내부

 i<5 

되려면. 그런 다음 원하는 결과를 얻습니다. 도움이 되길 바랍니다.