1

I am working on a project where I have to transform coordinates expressed in LUREF (LUxemburg REference System/Frame) to WGS84 (Lat/long) by using the bundle PostGIS of Postgres database.

The Proj4 texts for the different CRS are the following :

Source CRS

EPSG:9895

+proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1 +x_0=80000 +y_0=100000 +ellps=intl +units=m +no_defs +type=crs

Target CRS

EPSG:4937

+proj=longlat +ellps=GRS80 +no_defs +type=crs

Transformation

+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1 +x_0=80000 +y_0=100000 +ellps=intl +step +proj=cart +ellps=intl +step +proj=helmert +x=-189.228 +y=12.0035 +z=-42.6303 +rx=0.48171 +ry=3.09948 +rz=-2.68639 +s=0.46346 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1

I did the transformation in python, here is the code :

from pyproj import CRS,Transformer
from pyproj.transformer import TransformerGroup

## Check transformation parameters between 4937 and 9895
trans_group = TransformerGroup(CRS("EPSG:9895").to_3d(),CRS("EPSG:4937").to_3d())
trans_group

print(trans_group.transformers[0].description)
print(trans_group.transformers[0].to_proj4(pretty=True))

## compute transformation
# LUREF2020 Bursa is computing 3D information (input N,E,U)
print(trans_group.transformers[0].transform(92294, 59432, 0))

What I need to do, is doing the same transformation but with PostGIS using the bundle functions. I did it using ST_Transform but I got a result not correct :

ST_AsText(ST_Transform(st_pointz(59432,92294,0,9895), 4937) )

But with this method, I didn't use the expression of the transformation, I only used the two proj4 texts that describe the two CRSs. What I want to do, is to do the transformation with PostGIS but using the transformation text I wrote above.

Thank you in advance!

mabkoz
  • 13
  • 3
  • Where are you expecting the point to be plotted? I managed to get it somewhere close to LZR at the "Rue de l'Hôpital" in Luxembourg. Is it wrong? – Jim Jones Jun 09 '23 at 17:52
  • `POINT Z (5.881197954291744 49.763701272304665 0)` or if you don't care about Z: `SELECT ST_AsText( ST_Transform( ST_SetSRID(ST_MakePoint(59432,92294),9895), 4937) );` – Jim Jones Jun 09 '23 at 18:07
  • Hi Jim, thanks for your reply. I got the same results too, but that's not correct.. It has to return the result calculated in python which is : 5.882638904443646, 49.764768308262816. That's why I asked if there is a way to make a transformation based on transformation proj4 text with PostGIS. – mabkoz Jun 12 '23 at 07:44

0 Answers0