Specialist needed: Seven Parameter Helmert Transform

[Martin Isenburg]

Hello,

Is someone a specialist for the Seven Parameter Helmert Transform? I have LIDAR in EPSG 31259 (MGI / Austria GK M34) and try to get to EPSG 25833 (ETRS89 / UTM 33N) via this seven Seven Parameter Helmert Transform 

601.705 84.263 485.227 -4.7354 -1.3145 -5.3930 -2.388739

Here is the input in 31259 (also attached)

las2txt -i four_points.laz -stdout  

725363.66 200468.42 267.47
725374.33 200469.63 267.15
725360.41 200488.13 264.75
725373.60 200489.69 264.70

Let's just care about the x and y for now. First I set the z coordinate to 0

las2las -i four_points.laz -clamp_z 0 0 -odix _zero -olaz

las2txt -i four_points_zero.laz -stdout
725363.66 200468.42 0.00
725374.33 200469.63 0.00
725360.41 200488.13 0.00
725373.60 200489.69 0.00
 

Then I go to ECEF 

  
las2las -i four_points_zero.laz -epsg 31259 -target_ecef -odix _ecef -olaz

las2txt -i four_points_zero_ecef.laz -stdout

4192592.71 1202975.77 4636996.36
4192588.88 1202985.77 4636997.22
4192579.79 1202968.60 4637009.81
4192575.02 1202980.95 4637010.91

Next I apply the Seven Parameter Helmert Transform 

las2las -i four_points_zero_ecef.laz -remove_all_vlrs ^
-transform_helmert 601.705,84.263,485.227,-4.7354,-1.3145,-5.3930,-2.388739 ^
-odix _wgs84 -olaz
las2txt -i four_points_zero_ecef_wgs84.laz -stdout  

4193186.30 1203054.00 4637469.61
4193182.47 1203064.00 4637470.47
4193173.38 1203046.83 4637483.06
4193168.61 1203059.18 4637484.16

from this I go to this output in 25833

las2las -i four_points_zero_ecef_wgs84.laz ^
-ecef -wgs84 -target_epsg 25833 -odix _25833 -olaz
las2txt -i four_points_zero_ecef_wgs84_25833.laz -stdout
576756.89 5199316.82 45.66

576767.54 5199318.22 45.66
576753.31 5199336.47 45.66
576766.47 5199338.25 45.66

But according to the client the coordinates should be. Note again here I only care about the x and y coordinate for now.

576764.00 5199319.54 313.78
576774.64 5199320.93 313.46
576760.41 5199339.18 311.06
576773.57 5199340.97 311.01

Where is my mistake in the steps above?

Regards,

Martin

[jaime garbanzo]

Hi, Mating I was playing around with your transformation. 

I think the problem is the definition of the parameters. I think as a standard it is expected to have miliarc seconds in the rotation parameters but in this case it seems that they are seconds. 

601.705 84.263 485.227 -4.7354 -1.3145 -5.3930 -2.388739

This is my transformation. 

So I started from 

"las2las -i four_points_zero.laz -epsg 31259 -target_ecef -odix _ecef -olaz

las2txt -i four_points_zero_ecef.laz -stdout

4192592.71 1202975.77 4636996.36
4192588.88 1202985.77 4636997.22
4192579.79 1202968.60 4637009.81
4192575.02 1202980.95 4637010.91"

then I apply the Bursa-Wolf transformation (the one that does not use the barycenter)

I got

4193182.498 1203060.323 4637471.409
4193178.668 1203070.323 4637472.27
4193169.578 1203053.153 4637484.859
4193164.808 1203065.503 4637485.96

I transformed it to WGS84

lat,lon, z

46d 56m 34.6661s, 16d 0m  31.11549 s,  45.672

46d 56m, 34.70698s, 16d, 0m 31.61993s, 45.670

46d, 56m 35.30409s, 16d 0m, 30.95806s, 45.671

46d, 56m, 35.35632s, 16d, 0m, 31.58156s, 45.669

because I don't know the projection transformation I do it with the epsg.io tool (https://epsg.io/transform#s_srs=4326&t_srs=25833)

I get

x, y

576763.98,5199319.53

576774.63, 5199320.93

576760.40, 5199339.18

576773.56,  5199340.97

and these are your client coordinates. 

576764.00 5199319.54 313.78
576774.64 5199320.93 313.46
576760.41 5199339.18 311.06
576773.57 5199340.97 311.01

I hope that it helps.

My best.

Jaime

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe
---
You received this message because you are subscribed to the Google Groups "LAStools - efficient tools for LiDAR processing" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lastools+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/CABSWR-EGEiWMD3qoo4wbzKqAn%2BQEXvbugevkYwZHx-6e0CxmzA%40mail.gmail.com.

--

M. Sc. Jaime Garbanzo León

Instructor Invitado

Universidad de Costa Rica

[Martin Isenburg]

Hello Jaime,

I'm happy you are able to get to the right coordinates and thank you for looking into this for me, but I am not sure what you did, and how you performed a Bursa-Wolf transformation instead. The user wants to use the Seven Parameter Helmert Transform with these parameters. I must mess up something and it's not a milliarcseconds to arcseconds mixup.

The code is online and I think it assumes arcseconds. Or is there a bug in the rotational part of my code? I looked at it again and again but often it needs a fresh pair of eyes.

class LASoperationTransformHelmert : public LASoperation
	{
	public:
	inline const CHAR* name() const { return "transform_helmert"; };
	inline I32 get_command(CHAR* string) const { return sprintf(string, "-%s %lf,%lf,%lf,%lf,%lf,%lf,%lf ", name(), dx, dy, dz, rx, ry, rz, m); };
	inline U32 get_decompress_selective() const { return LASZIP_DECOMPRESS_SELECTIVE_CHANNEL_RETURNS_XY | LASZIP_DECOMPRESS_SELECTIVE_Z; };
	inline void transform(LASpoint* point) {
	F64 x = scale*( ( point->get_x())-(rz_rad*point->get_y())+(ry_rad*point->get_z())) + dx;
	F64 y = scale*( (rz_rad*point->get_x())+( point->get_y())-(rx_rad*point->get_z())) + dy;
	F64 z = scale*(-(ry_rad*point->get_x())+(rx_rad*point->get_y())+( point->get_z())) + dz;
	if (!point->set_x(x))
	{
	overflow++;
	}
	if (!point->set_y(y))
	{
	overflow++;
	}
	if (!point->set_z(z))
	{
	overflow++;
	}
	};
	LASoperationTransformHelmert(F64 dx, F64 dy, F64 dz, F64 rx, F64 ry, F64 rz, F64 m) { this->dx = dx; this->dy = dy; this->dz = dz; this->rx = rx; this->ry = ry; this->rz = rz; this->m = m; rx_rad = 4.84813681109536e-6*rx; ry_rad = 4.84813681109536e-6*ry; rz_rad = 4.84813681109536e-6*rz; scale = 1.0+(1.0e-6*m); };
	private:
	F64 dx, dy, dz, rx, ry, rz, m, rx_rad, ry_rad, rz_rad, scale;
	};

https://github.com/LAStools/LAStools/blob/master/LASlib/src/lastransform.cpp#L314

Regards,

Martin

To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/CAFdXu%3DbbmbwvNCx4vqzTibYtQG7B%2BqYuzjG5zQ-b4b4PCg-mmQ%40mail.gmail.com.

[Michael Perdue]

Hi Martin,

There is nothing wrong with your code. The problem is that the transformation parameters use a position vector rotation rather than the coordinate frame rotation convention. This is common for a lot of European countries. To get the right answers, just flip the sign of your rotations.

>las2las -epsg 31259 -i four_points.laz -o four_points_zero_ecef.laz -target_ecef

>las2las -i four_points_zero_ecef.laz -remove_all_vlrs -transform_helmert 601.705,84.263,485.227,4.7354,1.3145,5.3930,-2.388739 -odix _wgs84 -olaz

>las2las -i four_points_zero_ecef_wgs84.laz -ecef -wgs84 -target_epsg 25833 -odix _25833 -olaz

>las2txt -i four_points_zero_ecef_wgs84_25833.laz -stdout
576763.99 5199319.54 313.15
576774.64 5199320.93 312.83
576760.40 5199339.19 310.42
576773.56 5199340.97 310.38

Cheers,

Mike

To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/CABSWR-FWwRXVbuqsDyfcOEHX_QZJxW6Z4zo-LbCuF%2BHuT7dA%3DA%40mail.gmail.com.

[jaime garbanzo]

Hi Martin, 

yes Michael is right, your signs in the transformation are different from mine. As I know it, this that you're doing is a 7 parameter Helmert transformation and is known as Bursa-Wolf.  I send you a nice review on the topic just in case you need it.

My best, 

Jaime

	F64 x = scale*( ( point->get_x())+(rz_rad*point->get_y())-(ry_rad*point->get_z())) + dx;
	F64 y = scale*( -(rz_rad*point->get_x())+( point->get_y())+(rx_rad*point->get_z())) + dy;
	F64 z = scale*(+(ry_rad*point->get_x())-(rx_rad*point->get_y())+( point->get_z())) + dz;

To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/CAP-P_3JYsCeuHAK5Fc%2B%2BNTAPB02e9esvYqqeo7ENbYajUj-ytQ%40mail.gmail.com.

[Martin Isenburg]

Hello Jaime,

Thanks lot for your input. That fixed it. I owe you a beer but you have to come to Samara to get it ... (-:

In the next release of LAStools the new "lasdatum" tool will be able to do datum shifts using either an NTv2 grid or a seven parameter Helmert transform.

Regards,

Martin

To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/CAFdXu%3DYfjjZQHc8C99xoP16AcqOsLz3cTHR0QyP6mGOM9toPrQ%40mail.gmail.com.