towgs84 parameter

[otofoto]

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

[dr-ni]

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

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

[otofoto]

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]

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]

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]

GRS80.txt LGS21.txt

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

[dr-ni]

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]

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
LKS2.txt
:

[zvezdochiot]

24.10877524 56.94747212

Its no X Y Z! Its B and L!

[otofoto]

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

[zvezdochiot]

Helmert3d requires Z

See #5 (comment)

[dr-ni]

@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]

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]

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]

@dr-ni say:

&name

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

[dr-ni]

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

[dr-ni]

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

[otofoto]

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

[zvezdochiot]

@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]

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]

@otofoto say:

Strange results:

See #5 (comment)

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

[otofoto]

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

[zvezdochiot]

@otofoto say:

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

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

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

[dr-ni]

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

[dr-ni]

@otofoto

can you give some details about LGS21?

[zvezdochiot]

@dr-ni say:

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

Ok. Find.

helmert3d/src/helmblhtoxyz.c

Line 94 in 7b52d75

fprintf(stdout,"X-Y-Z -> B-H-L...\n");

[zvezdochiot]

@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]

and the bin helmblhtoxyz shold also be blh2xyz

[zvezdochiot]

@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]

but then

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

is wrong

[zvezdochiot]

@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]

@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]

@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]

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

[dr-ni]

@zvezdochiot

we should then either use helmblh2xyz or helmblhtoxyz

[zvezdochiot]

@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]

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

[zvezdochiot]

@dr-ni say:

helmellipsoid

helmeltrc (Helmert Ellipsoid Transform Coordinates).

[dr-ni]

helmeltrans

[zvezdochiot]

@dr-ni say:

helmeltrans

Or helmtransel?

[dr-ni]

I like helmeltrans

[otofoto]

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]

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

[otofoto]

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]

@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]

yes

[dr-ni]

@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]

@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]

ah ok I was missing the ellipsoid BESSEL for LGS21

[dr-ni]

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

Well done!

[zvezdochiot]

@dr-ni say:

I was only checking rz and DS

So the first three direct-flow:

419.6398186970 183.7356435974 549.3785472015

[otofoto]

Don't compare to my previous parameters it may differ slightly as it was calculated by less points.
Got such string using Excel also only ry needs sign change not rz,ry

+towgs84=419.6398111864,183.7356423902,549.3785528344,1.89700504716617,3.30502843140266,4.55748689876759,-2.66940000004379

Looks good now it is around 1 meter shift.