How can I convert the Gaussian coordinates( Gauß-Krüger) to the WGS84 coordinates

Question

I have a list of WGS84 and Gauss coordinates, but how to convert another WGS84 coordinates.
I want to use PROJ for conversion.

What is PROJ string like?

---

| Name | x         | y        | H       | Lat(N)      | Lon(E)      |  
|------|-----------|----------|---------|-------------|-------------|  
| X1   | -1676.071 | 4056.312 | 10.4167 | 31.2203175  | 121.5097212 |  
| X2   | -1180.277 | 6330.704 | 10.4985 | 31.22477859 | 121.533594  |  
| X3   | 472.067   | 6226.347 | 10.4846 | 31.23968209 | 121.5325094 |  
| X4   | -73.008   | 4895.714 | 10.4383 | 31.23477255 | 121.5185389 |  
| X5   | -646.523  | 5682.03  | 10.4711 | 31.22959614 | 121.5267891 |  
| X6   | -398.196  | 5181.716 | 10.4711 | ?           | ?           |  

Answer 1

First, let's see what the mapping between a Y-X plane and a Lon-Lat plane would look like:

Notes:

• Y-X-H coordinates are Cartesian ones.

• Since you named it Gaussian coordinates, it is suposed to Y and X are Transverse Mercator projections of the geodetic coordinates.

• I am assuming that Y coordinates are Eastings and X coordinates are Northings. Also, H coordinates seems to be orthometric heights. But in a projected system we can go from the ellipsoid to a 2D plane and vice versa.

• If you had also ellipsoidal heights we can try a transformation between Cartesian geocentric coordinates (X-Y-Z) and a X-Y-H Cartesian system. Or we can try the projection but using a vertical datum to transform altitudes in heights. For now, we can only discard H coordinates.

---

Then, we generate the text files to be processed by PROJ:

lonlat.txt:

121.5097212 31.2203175
121.533594 31.22477859
121.5325094 31.23968209
121.5185389 31.23477255
121.5267891 31.22959614

yx.txt:

4056.312 -1676.071
6330.704 -1180.277
6226.347 472.067
4895.714 -73.008
5682.03 -646.523
5181.716 -398.196

---

Let's try the most simple way to estimate the projection: Using one given point as the center of the projection plane. We will use X4 point because its relative spatial location:

C:\test>proj -f "%.3f" +proj=tmerc +lon_0=121.5185389 +lat_0=31.23477255 +x_0=4895.714 +y_0=-73.008 lonlat.txt
4055.528        -1675.655
6330.157        -1180.975
6226.607        471.414
4895.714        -73.008
5681.749        -646.905

We are defining a custom Transverse Mercator projection, centered on X4, with the X4 coordinates defining false Easting a Northing.

We can see indeterminations in the order of the meter. Analyzing the errors for each point, I think we could calculate the coordinates of X6 with an indeterminacy close to half a meter.

C:\test>proj -I -f "%.8f" +proj=tmerc +lon_0=121.5185389 +lat_0=31.23477255 +x_0=4895.714 +y_0=-73.008 yx.txt
121.50972942    31.22031375
121.53359975    31.22478488
121.53250667    31.23968798
121.51853890    31.23477255
121.52679205    31.22959959
121.52154084    31.23183955  

So, I think that X6 is almost a half meter away from LonLat = (121.52154084 31.23183955).

---

We can start playing with tests. For example, in a projection centered in X4, with X4 coordinates as false Easting and Northing, where is the Y-X = (0 0) point located?

C:\test>proj -I -f "%.8f" +proj=tmerc +lon_0=121.5185389 +lat_0=31.23477255 +x_0=4895.714 +y_0=-73.008
0 0
121.46715049    31.23542076

So, can we use LonLat = (121.46715049 31.23542076) as the center of the projection without false Easting and Northing?

C:\test>proj -f "%.3f" +proj=tmerc +lon_0=121.46715049 +lat_0=31.23542076 +x_0=0 +y_0=0 lonlat.txt
4056.307        -1673.769
6330.706        -1178.031
6226.388        474.311
4895.748        -70.731
5682.050        -644.263

We are seeing very good approximations in the Y coordinate, but very bad in the X.

This must be because when moving away from the central meridian, the latitudes deform. But I recognize that I was surprised by the magnitude of these deformation (also, by the accuracy in the Y outputs).

We could try assuming that the projection center is at Y = 5000? In a projection centered in X4, with X4 coordinates as false Easting and Northing, where is the Y-X = (5000 0) point located?

C:\test>proj -I -f "%.8f" +proj=tmerc +lon_0=121.5185389 +lat_0=31.23477255 +x_0=4895.714 +y_0=-73.008
5000 0
121.51963355    31.23543103

Let's see then project using that center with false Easting 5000:

C:\test>proj -f "%.3f" +proj=tmerc +lon_0=121.51963355 +lat_0=31.23543103 +x_0=5000 +y_0=0 lonlat.txt
4055.512        -1675.646
6330.145        -1180.989
6226.612        471.402
4895.713        -73.007
5681.743        -646.912

Ok, very similar than project centered in X4.

---

Conclusion:

• We can easily estimate a custom projection with an indetermination of the order of one meter, with and for the data given.

• Accuracy in Eastings and unaccuracy in Northings, when moving away the central meridian, attracts attention. An analysis is still pending in order to find a better solution to the problem.

• The recommended action is to consult the information provider about the type of transformation and the parameters used.

Answer 2

A working formula for this problem would be the following (in JavaScript Code), it was found and rewritten to JavaScript from here: https://community.qlik.com/t5/QlikView-Documents/Convert-Gauss-Kruger-GK-coordinates-to-Latitude-Longitude/ta-p/1487143

let XCOORD_Field // (Rechtswert)
let YCOORD_Field // (Hochwert)

let rho = 57.29577951; //=180/PI
let e2 = 0.006719219;
let b1 = YCOORD_Field/10000855.7646;
let b2 = Math.pow(b1,2);
let bf = 325632.08677*b1*((((((0.00000562025*b2-0.0000436398)*b2+0.00022976983)*b2-0.00113566119)*b2+0.00424914906)*b2-0.00831729565)*b2+1) / 3600/ rho;
let fa = ((XCOORD_Field-(Math.floor(XCOORD_Field/1000000)*1000000)-500000)/(6398786.849/Math.sqrt(1+(Math.pow(Math.cos(bf),2)*e2))));
let LAT_Formula = (bf-(Math.pow(fa,2)*(Math.sin(bf)/Math.cos(bf))*(1+(Math.pow(Math.cos(bf),2) * e2))/2)
    +(Math.pow(fa,4)*(Math.sin(bf)/Math.cos(bf))*(5+(3*Math.pow(Math.sin(bf)/Math.cos(bf),2))+(6*(Math.pow(Math.cos(bf),2)
    * e2))-(6*(Math.pow(Math.cos(bf),2) * e2)*Math.pow(Math.sin(bf)/Math.cos(bf),2)))/24)) * rho;
let LON_Formula = ((fa-(Math.pow(fa,3)*(1+(2*Math.pow(Math.sin(bf)/Math.cos(bf),2))+(Math.pow(Math.cos(bf),2)
    * e2))/6)+(Math.pow(fa,5)*(1+(28*Math.pow(Math.sin(bf)/Math.cos(bf),2))+(24*Math.pow(Math.sin(bf)/Math.cos(bf),4)))/120))
    * rho/Math.cos(bf)) + (Math.floor(XCOORD_Field/1000000)*3);

return [LAT_Formula, LON_Formula]