Converting Belge 1972 / Belgian Lambert 72 (EPSG:31370) coordinates to WGS 84 (EPSG:4326)

Question

I am trying to find a code that can convert Lambert72 coordinates to WGS84, so that I can get the same results as I get in this site, going to Menu APPs->Transform Coordinates.

https://mygeodata.cloud/cs2cs/

As an example I tried the following coordinates pair in this site, choosing on the left, the code 31370 (Lambert72) and on the right, WGS84.: 149334.41 167411

The result is: 4.35930680453; 50.817136997

For this result to work on Google Maps, the values must be switched, so the result I need is rather: 50.817136997; 4.35930680453

I tried the code from this answer in a similar post but I converted to C instead of C++/C#:

https://stackoverflow.com/a/40589076/1911497

So, my code, modified to C is this:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static void Lambert72toWGS84latlong(double X, double Y);

int main(){
    double x, y;
    
    x = 149334.41, y = 167411.0;
    
    Lambert72toWGS84latlong(x, y);

    return 0;
}


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

 double eLamb = sqrt(eCarre);
 double eSur2 = eLamb / 2.0;

 double Tan1 = (X - 150000.012) / (5400088.437 - Y);
 double Lambda = LongRef + (1.0 / nLamb) * (0.000142043 + atan(Tan1));
 double RLamb = sqrt(pow((X - 150000.012) , 2.0) + pow   ((5400088.437 - Y) ,2.0));

 double TanZDemi = pow((RLamb / KLamb),(1.0 / nLamb));
 double Lati1 = 2.0 * atan(TanZDemi);

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

 double lat, lng ;
 int i=0; 
 do{
   eSin = eLamb * sin(Lati1);
   Mult1 = 1.0 - eSin;
   Mult2 = 1.0 + eSin;
   Mult = pow((Mult1 / Mult2) , (eLamb / 2.0));
   LatiN = (M_PI / 2.0) - (2.0 * (atan(TanZDemi * Mult)));
   Diff = LatiN - Lati1;
   printf("Diff: %d\n", abs(Diff));
   //Lati1 = LatiN;
   printf("Iterations: %d\n", i++);
  }while (abs(Diff)> 0.0000000277777);


  lat=LatiN;
  lng=Lambda;

  double SinLat = sin(lat);
  double SinLng = sin(lng);
  double CoSinLat = cos(lat);
  double CoSinLng = 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 / sqrt(1.0 - LWe2 * SinLat * SinLat);    
  double Rm = LWa * (1 - LWe2) /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) / M_PI;
  double LngWGS84 = ((lng + DLng) * 180.0) / M_PI;

  printf("Latitude: %.8f\t Longitude: %.8f\n", LatWGS84, LngWGS84);
}

However I don't get the same result. The results I get are: Latitude: 50.78274764 Longitude: 4.35930689

What would I need to change to be able to get the same result as the one returned by the mentioned site above?

Thanks PsySc0rpi0n

Sorry, will edit the original post and will add the result from that function. thanks Psy — PsySc0rpi0n, Jul 21 '21 at 15:01
Why did you comment out `Lati1 = LatiN;` (in the `do` loop)? — Scott Hunter, Jul 21 '21 at 15:11
That was me just testing what would be the difference without that line, because if you don't comment that line, the while loop won't ever run more than once, if I'm not mistaken. — PsySc0rpi0n, Jul 21 '21 at 15:20
I did that, I replied in comments in the answer below. It's close but still not exactly the same. Tiny differences here makes quite some difference. This is the difference in the map: https://i.ibb.co/nsSMsL0/imagem.png — PsySc0rpi0n, Jul 21 '21 at 15:34
*play with* `FLT_EVAL_METHOD`? See https://stackoverflow.com/questions/28851350/same-flt-eval-method-different-results-in-gcc-clang ... but I'm inclined to say there's not enough error propagation to matter. — pmg, Jul 21 '21 at 15:51
Have you verified that the original C++ produces the correct result? If it doesn't, your translation doesn't have much of a chance. — Scott Hunter, Jul 21 '21 at 15:52
What's the point of `Rm + 0.0` an later `Rn + 0.0`? Are those typos or the original code required some type conversions? — Bob__, Jul 21 '21 at 15:54
@pmg I will see what that constant means and I'll try to play with it to see if result changes. — PsySc0rpi0n, Jul 22 '21 at 07:17
@ScottHunter I didn't try the original code because I'm not confortable with c++. That's why I did the transaltion to C. Anyway, the only difference from one code to the other is in 2 functions and a constant. Math.Abs to fabs(), Math.Pow to pow() and Math.PI to M_PI. These were the only things I changed. — PsySc0rpi0n, Jul 22 '21 at 07:17
@Bob__ That's in the original code. I'm not familiar with the steps to convert from one system of coordinates to the other, so I just used the original code as is, except from the transalation. — PsySc0rpi0n, Jul 22 '21 at 07:17
Also, maybe try making `long double` variables (and operations... remember `sqrtl()`, `fabsl()`, `powl()`, ...) in your translated code. I don't have a C# compiler/debugger... would be fun to statement-by-statement debug both versions side-by-side :-) — pmg, Jul 22 '21 at 07:26
On my computer (where `FLT_EVAL_METHOD` is `0`) your code gives the wanted result (in 3 iterations) with `//Lati1 = LatiN;` uncommented and using `fabs()`. — pmg, Jul 22 '21 at 07:44
@pmg I'll try to change the functions to those you mention. I was also using an online C compiler, so I'm not sure I can change the FLT_EVAL_METHOD const variable in the code or if it's an OS env variable which I would need (if i'm not mistaken) access to a shell to be able to change it. — PsySc0rpi0n, Jul 22 '21 at 11:15
@pmg I tried with compiler flag `-mfpmath=387`, which I think it is equivalent to `FLT_EVAL_METHOD 2` and with the sqrtl(), powl() and fabsl() functions The result is still not the one I need. I get: **Latitude: 50.81742339 Longitude: 4.35956292** — PsySc0rpi0n, Jul 22 '21 at 11:27
@pmg for the input coordintes x =149334.41, y = 167411.0 you get **Latitude: 50.81713876 Longitude: 4.35930782** but I need them to be **50.817136997, 4.35930680453**. Sorry if I already mess up with the input and output values. I think, in the meantime, I tried other input values and was replying above for those new values. — PsySc0rpi0n, Jul 22 '21 at 12:51

pmg · Answer 1 · 2021-07-21T18:29:19.707

1

do {
    // stuff
} while (abs(Diff)> 0.0000000277777);
//     **int abs(int);**

You probably want

 do {
     // stuff
 } while (fabs(Diff) > 0.0000000277777);
 //     **double fabs(double);**

edited Jul 21 '21 at 18:29

answered Jul 21 '21 at 15:13

pmg

106,608
13
126
198

That is a necessary, but not sufficient, fix. – Scott Hunter Jul 21 '21 at 15:17
The original code is Math.Abs(). I just used what was in the original code but converted to C. Will try though and will report back. – PsySc0rpi0n Jul 21 '21 at 15:22
1

[C#'s `Math.Abs()`](https://learn.microsoft.com/en-us/dotnet/api/system.math.abs?view=net-5.0) is different than [C's `abs()`](https://pubs.opengroup.org/onlinepubs/9699919799/functions/abs.html) – pmg Jul 21 '21 at 15:27
1

Ok, with fabs() it's closer but still not good enough to what I need. I get: Latitude: 50.81713876 Longitude: 4.35930782 @pmg ok I didn't know that. Thanks – PsySc0rpi0n Jul 21 '21 at 15:28

Converting Belge 1972 / Belgian Lambert 72 (EPSG:31370) coordinates to WGS 84 (EPSG:4326)

1 Answers1