如何在PostGIS中使用指定Proj4转换文本实现LUREF到WGS84转换

背景 #

需要将LUREF（卢森堡参考系统，EPSG:9895）坐标转换为WGS84经纬度（EPSG:4937），直接使用PostGIS的ST_Transform调用两个CRS得到的结果不正确，但已通过Python pyproj使用自定义转换管道实现了正确转换，需在PostGIS中复用该转换逻辑。

相关参数 #

源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

目标CRS（EPSG:4937） #

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

自定义转换管道 #

+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

PostGIS实现步骤 #

1. 注册自定义坐标转换 #

使用ST_CreateCoordinateOperation将自定义转换管道注册到PostGIS的空间参考系统中：

-- 创建名为lueref_to_wgs84_custom的自定义转换
SELECT ST_CreateCoordinateOperation(
  'EPSG:9895',
  'EPSG:4937',
  'lueref_to_wgs84_custom',
  '<http://www.opengis.net/def/coordinateOperation/OGC/1.0/pipeline>',
  'PROJ4:' || '+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'
);

2. 执行自定义转换 #

在ST_Transform中指定使用刚注册的自定义转换名称，完成坐标转换：

-- 转换示例点（E:59432, N:92294, U:0）
SELECT ST_AsText(
  ST_Transform(
    ST_PointZ(59432, 92294, 0, 9895),
    4937,
    'lueref_to_wgs84_custom'
  )
);

3. 结果验证 #

可对比PostGIS转换结果与pyproj的输出，确保一致性。pyproj参考代码：

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

# 查看转换参数详情
trans_group = TransformerGroup(CRS("EPSG:9895").to_3d(),CRS("EPSG:4937").to_3d())
print(trans_group.transformers[0].description)
print(trans_group.transformers[0].to_proj4(pretty=True))

# 执行转换（输入顺序为N,E,U）
print(trans_group.transformers[0].transform(92294, 59432, 0))

内容的提问来源于stack exchange，提问作者mabkoz