[Proj] trying to understand results of reprojecting from epsg 4326 to epsg 31370

[Noel Zinn]

Consider changing the signs of the three rotations if you can do that 
manually, or the rotation sense of the 7-parameter transformation.
-Noel

----- Original Message ----- 
From: "Moritz Lennert" <mlennert at club.worldonline.be>
To: <Proj at lists.maptools.org>
Sent: Tuesday, July 07, 2009 9:14 AM
Subject: [Proj] trying to understand results of reprojecting from epsg 4326 
to epsg 31370


Hello,

I'm trying to reproject the following WGS84 Lat-Long (EPSG 4326)
coordinates to the Belgian 1972 projection EPSG 31370:

4.468493 50.856057

Using the official software from the Belgian National Geographic
Institute, I get the following coordinates (confirmed by a contact in
the Institute):

157022.870 171745.084

Using cs2cs +init=epsg:4326 +to +init=epsg:31370, I get

157201.51    171803.24

i.e. way off the official response.

31370 is defined as

+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
+towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m
+no_defs

so datum transformation is defined.

The transformation parameters are confirmed by the Belgian Institute
(http://www.ign.be/FR/FR2-1-2.shtm), which considers ETRS89 to be
equivalent to WGS84:

Paramètres de ETRS89 vers BD72
Tx 106.868628
Ty -52.297783
Tz 103.723893
Rx 0.336570
Ry -0.456955
Rz 1.842183
scale 1.0000012747

As these parameters are in the direction fromwgs84, they should be
negated to get the inverse transformation. So, IIUC, 31370 should
actually be defined as:

+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
+towgs84=-106.869,52.2978,-103.724,-0.33657,0.456955,-1.84218,1 +units=m
+no_defs

When I try this, I get:
157080.62 171775.09

Much close, but still no good. The Institute also use a correction grid,
but this should only change accuracy at the centimeter level.

Interestingly enough, when I use the 3-parameter transformation given by
EPSG 1610:

+proj=lcc +lat_2=51.16666723333333 +lat_1=49.8333339 +lat_0=90
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
+towgs84=-125.8,79.9,-100.5 +units=m +no_defs

I get perfect results:
157022.80    171745.11

Does anyone have an explanation for this behaviour ?

Moritz
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj