[postgis-users] Coordinate transformation problem

[Arnaud Lesauvage]

Jan Hartmann a écrit :
> Hi Arnaud,
> 
> I don't know about the Belgian situation, but we had the same problem in 
> Holland. The problem is that the National Grid has a different Datum 
> than WGS84: it uses a slightly different ellipsoid to approximate the 
> surface of the earth. With the standard PROJ epsg-file, this Datum 
> difference will *not* be computed, and you get differences up to 100 
> meters. For a complete transformation to WGS84 you need an extra 
> "+towgs" parameter in the PROJ parameter string. How this looks for 
> Holland, you can see in
>  
> http://postgis.refractions.net/pipermail/postgis-users/2005-November/009894.html
> 
> The EPSG number for the Belgian national grid is 31370. In the standard 
> PROJ epsg-file, and in the "spatial_ref_sys" table in PostGIS, this 
> number has the following parameters:
> 
> +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
> +units=m +no_defs
> 
> As I said, you need extra numbers for the Datum conversion. The Belgian 
> Geodetic Service (http://www.ngi.be/FR/FR4-4.shtm) gives the following 
> values for transformations from the Belgian grid to ETRS89 / WGS84:
> 
> Tx 	106.868628
> Ty 	-52.297783
> Tz 	103.723893
> échelle 1+S 	1.0000012747
> 
> Rotations (secondes sexagésimales)
> Rx 	0.336570
> Ry 	-0.456955
> Rz 	1.842183
> 
> This would result in the following PROJ-string:
> 
> +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
> +towgs=106.868628,-52.297783,103.723893,0.336570,-0.456955,1.842183,1.0000012747
> +units=m +no_defs
> 
> To test this, you need to replace the entry for EPSG:31370 in the 
> PostGIS  "spatial_ref_sys" table with this string (make a backup!). It 
> is possible that you have to reverse the sign of the fourth, fifth, and 
> sixth number of the towgs parameter:
> 
> +towgs=106.868628,-52.297783,103.723893,-0.336570,0.456955,-1.842183,1.0000012747
> 
> I don't have Belgian data, so I cannot test this. Let me know if this 
> works (or not!)


Hi Jan !
Thansk a lot for your answer ! I am now quite confident that 
the problem comes from this missing 'towgs' parameter !

I still have the same problem though : I see absolutely no 
difference when I add this parameter !
I had inserted the projection from the .prj file into my 
spatial_ref_sys table (SRID 31376), so I updated it.

tmp=# SELECT srid, proj4text FROM spatial_ref_sys WHERE srid 
IN (31376,31370);
srid   |      proj4text
-------+-----------------------------------------------
  31370 | +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 +units=m +no_defs

  31376 | +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 
+towgs=-99.1,53.3,-112.5,0.419,-0.830,1.885,-1 +units=m +no_defs
(2 rows)

My geometry column was already created in this SRID, so I 
thought that just changing the proj4text would change the 
result of the projection, but I have :

tmp=# SELECT
tmp-# Equals(
tmp-# Transform(SetSRID(cabu_geometry,31376),32631),
tmp-# Transform(SetSRID(cabu_geometry,31370),32631))
tmp-# FROM cadastre.cabu LIMIT 10;
  equals
--------
  t
  t
  t
  t
  t
  t
  t
  t
  t
  t
(10 rows)

So it looks like the projection from both SRIDs to 32631 is 
the same, doesn't it ?

--
Arnaud