dr-ni / helmert3d

The 3D Helmert transformation is a frequently used method in geodesy to produce transformations between different Cartesian coordinate systems
http://helmparms3d.sourceforge.net/
Other
3 stars 8 forks source link

towgs84 parameter #5

Open otofoto opened 2 years ago

otofoto commented 2 years ago

Option to calculate valid 7 parameter +towgs84 value that could be used with proj.4/WKT

dr-ni commented 2 years ago

These tools are designed only for 3d similarity transformation. Maybe gdaltransform is what you are looking for?

https://gdal.org/programs/gdaltransform.html

otofoto commented 2 years ago

No I am looking for tool to calculate position vector transformation parameters(towgs84) for setting coordinate system in QGIS as without these parameters there is a 100m shift on map transformation which is wrong.

dr-ni commented 2 years ago

can you give me a full example with all points, transformation parameters and your needed towgs84 values that I can understand what you want?

otofoto commented 2 years ago

Input: helmert3d-1.0.1-win32\helmparms3d.exe" LGS21.txt GRS80.txt out.txt

result: .... +towgs84=419.640956,183.735831,549.377809,1.89693297,3.30507103,4.55745697,-2.6694418

PROJ.4 uses the Position Vector" standard. Angles from matrix elements:

Rz = -r12 * 3600 *180 / PI
Ry = r13 * 3600 * 180 / PI
Rx = -r23 * 3600 * 180 / PI

The scale factor must be converted to scale difference in ppm:

DS = (s-1) * 1000000

No need to change translation vector components.

The seven parameter case uses delta_x, delta_y, delta_z, Rx - rotation X, Ry - rotation Y, Rz - rotation Z, M_BF - Scaling. The three translation parameters are in meters as in the three parameter case. The rotational parameters are in seconds of arc. The scaling is apparently the scale change in parts per million. A more complete discussion of the 3 and 7 parameter transformations can be found in the EPSG database (trf_method's 9603 and 9606). Within PROJ.4 the following calculations are used to apply the towgs84 transformation (going to WGS84). The x, y and z coordinates are in geocentric coordinates. Three parameter transformation (simple offsets):

  x[io] = x[io] + defn->datum_params[0];
  y[io] = y[io] + defn->datum_params[1];
  z[io] = z[io] + defn->datum_params[2];

Seven parameter transformation (translation, rotation and scaling):

  #define Dx_BF (defn->datum_params[0])
  #define Dy_BF (defn->datum_params[1])
  #define Dz_BF (defn->datum_params[2])
  #define Rx_BF (defn->datum_params[3])
  #define Ry_BF (defn->datum_params[4])
  #define Rz_BF (defn->datum_params[5])
  #define M_BF  (defn->datum_params[6])

  x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
  y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
  z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;

Note that EPSG method 9607 (coordinate frame rotation) coefficients can be converted to EPSG method 9606 (position vector 7-parameter) supported by PROJ.4 by reversing the sign of the rotation vectors. The methods are otherwise the same.

LGS21.txt GRS80.txt

zvezdochiot commented 2 years ago

GRS80.txt LGS21.txt

"This" by any chance latitude-longitude? Actually, need X-Y-Z!

X-Y-Z
NB

dr-ni commented 2 years ago

Why are you using a 3d transformation for a 2d usecase in your example? Do I understand this correct that you want to use the 3d helmert transformation to make an aproximated transformation between different coordinate systems?

otofoto commented 2 years ago

Yes need to find transformation parameters from Bessel(Cassini-soldner.txt) to WGS84(LKS2.txt)

  1. source +proj=cass +lat_0=56.94831083 +lon_0=24.10886056 +k=1 +x_0=0 +y_0=0 +ellps=bessel +units=m
  2. target EPSG:3059/EPSG3059

    
    gdaltransform.exe -s_srs EPSG:4661 -t_srs "+proj=cass +lat_0=56.94831083 +lon_0=24.10886056 +k=1 +x_0=0 +y_0=0  +ellps=bessel  +units=m"
    Enter X Y [Z [T]] values separated by space, and press Return.
    24.10877524 56.94747212
    result: **-5.19174698205156 -93.3869234910952 0**
    expected: **0.0 0.0 0**

gdaltransform -s_srs EPSG:3059 -t_srs "+proj=cass +lat_0=56.94831083 +lon_0=24.10886056 +k=1 +x_0=0 +y_0=0 +ellps=bessel +units=m" Enter X Y [Z [T]] values separated by space, and press Return. 311544.05 506617.17 result: -5.19203592124977 -93.388642267098 0 expected: 0.0 0.0 0

gdaltransform -s_srs EPSG:3059 -t_srs "+proj=cass +lat_0=56.94831083 +lon_0=24.10886056 +k=1 +x_0=0 +y_0=0 +towgs84=419.640956,183.735831,549.377809,1.89693297,3.30507103,4.55745697,-2.6694418 +ellps=bessel +units=m" Enter X Y [Z [T]] values separated by space, and press Return. 311544.05 506617.17 result almost as expected:0.99 -0.84 0.06 expected: 0.0 0.0 0


[Cassini-soldner.txt](https://github.com/dr-ni/helmert3d/files/7599022/Cassini-soldner.txt)
[LKS2.txt](https://github.com/dr-ni/helmert3d/files/7599023/LKS2.txt)
:
zvezdochiot commented 2 years ago

24.10877524 56.94747212

Its no X Y Z! Its B and L!

otofoto commented 2 years ago

I attached XYZ also but seems need to find different way as Helmert3d requires Z.

zvezdochiot commented 2 years ago

Helmert3d requires Z

See https://github.com/dr-ni/helmert3d/issues/5#issuecomment-978223652

dr-ni commented 2 years ago

@zvezdochiot wow super fast nice :-) @otofoto does this new tool now solve your problem?

GRS80.txt:

GRS80 6378137.0 6356752.3141
56.32258375 21.22652722 0
56.4935 21.783 0
56.61155556 22.93708333 0
56.48455556 23.55005556 0
56.48191667 24.31558333 0
56.59155556 25.15272222 0
56.25533333 25.54025 0
56.01013889 26.22444444 0
55.69546942 26.61753336 0
56.94747222 24.10877528 0
56.94747556 24.1087875 0

helmblhtoxyz GRS80.txt XYH.txt

*******************************
*      helmblhtoxyz v1.0.2      *
*   (c) U. Niethammer 2020    *
*******************************
Reading points...
Found 11 points
Starting calculate...
Ellipsoid GRS80, a = 6378137.000000, b = 6356752.314100
...done
Results written to XYH.txt
pi@rp4:~/Downloads $ cat XYH.txt
3304502.114760 1283492.332909 5284443.399514
3277159.672088 1309642.578748 5294972.981868
3240015.717909 1371108.106641 5302218.527351
3235978.291148 1410406.235492 5294423.110409
3217068.860883 1453615.777477 5294260.856490
3186268.667455 1496134.135442 5300992.627323
3204220.092028 1531098.564097 5280287.448084
3206044.070934 1579268.610385 5265073.416838
3221055.755008 1614217.624473 5245407.547179
3182737.947107 1424292.525098 5322712.005840
3182737.358780 1424293.076570 5322712.208698
zvezdochiot commented 2 years ago

Hi @dr-ni .

Do not rush. I am already redoing to the format (1.0.3):

helmblhtoxyz commands blh|xyz_src_infilename ellipsoid_src_infilename [blh|xyz_dest_infilename]

commands:
 xyz - convert B-L-H to X-Y-Z
 blh - convert X-Y-Z to B-L-H

It will be ready soon.

dr-ni commented 2 years ago
make 
gcc -Wall -O2 -o helmert3d src/helmert3d.c -lm 
gcc -Wall -O2 -o helmparms3d src/helmparms3d.c libsvdm.a -lm 
gcc -Wall -O2 -o helmdiff3d src/helmdiff3d.c -lm 
gcc -Wall -O2 -o helmblhtoxyz src/helmblhtoxyz.c -lm 
src/helmblhtoxyz.c: In function ‘main’:
src/helmblhtoxyz.c:102:36: warning: format ‘%s’ expects argument of type ‘char *’, but argument 3 has type ‘char (*)[128]’ [-Wformat=]
             stat = sscanf( ibuf, "%s %lf %lf", &name, &a, &b);
                                   ~^           ~~~~~

and some spaces wrong

*******************************
*      helmblhtoxyz v1.0.2      *
*   (c) U. Niethammer 2020    *
*******************************
*******************************
*       blh2xyz v1.0.2        *
*   (c) zvezdochiot 2021      *
*******************************
zvezdochiot commented 2 years ago

@dr-ni say:

&name

I know. Rush. Therefore, I want to quickly convert to name.

dr-ni commented 2 years ago

Starting calculate... => Starting calculation...

dr-ni commented 2 years ago

I think we should leave this issue open since it may be interesting for others, too

otofoto commented 2 years ago

So with helmblhtoxyz it will be possible to calculate geocentric coodinates and then find helmert parameters between them?

zvezdochiot commented 2 years ago

@otofoto say:

will be possible to calculate

It's already possible. It will be possible to return soon. See https://github.com/dr-ni/helmert3d/releases/tag/1.0.2

otofoto commented 2 years ago

Strange results:

helmblhtoxyz
19.710868 7.877001 0.000000
19.548217 8.272404 0.000000
19.458599 8.481061 0.000000
19.343550 8.740282 0.000000
19.213787 9.021964 0.000000
19.152327 9.151712 0.000000
19.041679 9.379761 0.000000
18.976880 9.510178 0.000000
19.374972 8.670405 0.000000
19.374970 8.670409 0.000000
zvezdochiot commented 2 years ago

@otofoto say:

Strange results:

See https://github.com/dr-ni/helmert3d/issues/5#issuecomment-979253526

"Watch yourself, be careful." (c) Viktor Tsoi.

otofoto commented 2 years ago

понял не дурак, дурак бы не понял :)

zvezdochiot commented 2 years ago

@otofoto say:

дурак бы не понял

See https://github.com/dr-ni/helmert3d/releases/tag/1.0.3

"И не будет после нас тьмы..." (c) @zvezdochiot .

dr-ni commented 2 years ago

@zvezdochiot X-Y-Z -> B-H-L... should be X-Y-Z -> B-L-H...

dr-ni commented 2 years ago

@otofoto

can you give some details about LGS21?

zvezdochiot commented 2 years ago

@dr-ni say:

X-Y-Z -> B-H-L...

Ok. Find. https://github.com/dr-ni/helmert3d/blob/7b52d75ff14096d05bdc06a1cc1de2ee6997da22/src/helmblhtoxyz.c#L94

zvezdochiot commented 2 years ago

@dr-ni say:

details about LGS21?

Try to apply GRS80. If the point in Helmert comply, it means LGS21 is a shifted GRS80. (Although there was some mention about Bessel!)

dr-ni commented 2 years ago

and the bin helmblhtoxyz shold also be blh2xyz

zvezdochiot commented 2 years ago

@dr-ni say:

and the bin helmblhtoxyz shold also be blh2xyz

Not! I will not do it exactly! You do, but I will be against.

dr-ni commented 2 years ago

but then

*******************************
*       blh2xyz v1.0.3        *
*   (c) U. Niethammer 2020    *
*******************************

is wrong

zvezdochiot commented 2 years ago

@dr-ni say:

is wrong

No. One thing is the inner name that is closest to the appointment. Another thing is to solve problems with the imposition of other packages, after which you remember what is her name. All utils specially have prefix helm!

otofoto commented 2 years ago

@otofoto

can you give some details about LGS21?

It is Cassini-Soldner projection on Bessel ellipsoid(Bessel-1841 either 6377397.155 6356078.963 or 6377397.155 6356078.96325 with azimuth to Jelgava 215°24'04.38[215.40121667]) lv: Page -135 AUGSTĀKĀ ĢEODĒZIJA ru:Page 33 - К- Ю. МЕНЗИН en: ASPRS PE&RS Article 10-20-GD-Latvia.pdf

I have known affine parameters but only for BL-BL, there are no direct LGS21 transformation to WGS84/EPSG:3059. For some reason and my limited knowledge just providing LGS21 as proj string ""+proj=cass +lat_0=56.94831083 +lon_0=24.10886056 +k=1 +x_0=0 +y_0=0 +ellps=bessel +units=m" resulting coordinate are shifted in any GIS/CAD or GDAL.

I really appreciate your help.

zvezdochiot commented 2 years ago

@otofoto say:

It is Cassini-Soldner projection on Bessel ellipsoid(Bessel-1841 either 6377397.155 6356078.963 or 6377397.155 6356078.96282 with azimuth to Jelgava 215°24'04.38[215.40121667])

Then a question arose. See PS: #7 .

otofoto commented 2 years ago

input format is ok as ellipsoid sizes can be easily found.

dr-ni commented 2 years ago

@zvezdochiot

we should then either use helmblh2xyz or helmblhtoxyz

zvezdochiot commented 2 years ago

@otofoto say:

input format is ok as ellipsoid sizes can be easily found.

Yes. And many programs contain exactly (a b). But the b is the derivative! Astro-geodetic methods determine (a 1/f). In normative documents, (a 1/f) are required. But in most programs, (a b).

@dr-ni say:

we should then either use helmblh2xyz or helmblhtoxyz

Ok. helmblh2xyz. But it is inconvenient to type.

dr-ni commented 2 years ago

ok better name for it is helmellipsoid then the options make also sense because helmblhtoxyz blh does not make sense

zvezdochiot commented 2 years ago

@dr-ni say:

helmellipsoid

helmeltrc (Helmert Ellipsoid Transform Coordinates).

dr-ni commented 2 years ago

helmeltrans

zvezdochiot commented 2 years ago

@dr-ni say:

helmeltrans

Or helmtransel?

dr-ni commented 2 years ago

I like helmeltrans

otofoto commented 2 years ago

IMHO "helm el" sounds like helmel trans and not related to helmert. Doesn't it translate to geocentric coordinates and is Helmert actually used there? BLH_XYZ or Gds_trf

dr-ni commented 2 years ago

zvezdochiot wants that all tools are beginnin with "helm" maybe we use helmprojection

otofoto commented 2 years ago

helmreprojection or helmtrf or helmtransform or maybe just instead of new module add blh and ellipsoid options in helmparms3d which sets input data type.

zvezdochiot commented 2 years ago

@dr-ni say:

I like helmeltrans

Die is cast. Discussion is closed!

*******************************
*     helmeltrans v1.0.4      *
*   (c) U. Niethammer 2020    *
*******************************

Should I put 2021 in copyright?

dr-ni commented 2 years ago

yes

dr-ni commented 2 years ago

@otofoto zvezdochiot has now prepared a new version :-) can you please try if you can now generate your needed towgs84 parameters please post all results after each step

zvezdochiot commented 2 years ago

@dr-ni say:

can you please try if you can now generate your needed towgs84 parameters

Why wait? You can do it yourself. The data is on hand:

0.9999999996 -0.0000220952 0.0000160236 
0.0000220953 0.9999999997 -0.0000091965 
-0.0000160234 0.0000091969 0.9999999998 
419.6398186970 183.7356435974 549.3785472015 
0.9999973306

LGS21-GRS80.zip

@otofoto say:

+towgs84=419.640956,183.735831,549.377809,1.89693297,3.30507103,4.55745697,-2.6694418

dr-ni commented 2 years ago

ah ok I was missing the ellipsoid BESSEL for LGS21

dr-ni commented 2 years ago

seems to be looking good (I was only checking rz=4,559773758 and DS=−2,6694)

Well done!

zvezdochiot commented 2 years ago

@dr-ni say:

I was only checking rz and DS

So the first three direct-flow:

419.6398186970 183.7356435974 549.3785472015