如何使用 GeoTools MathTransform？

Question

我需要使用 GeoTools Java 库将点从 EPSG:4312 转换为 WGS84。 但我不确定我是否正确使用它，我也不确定它是否提供了正确的结果。

请参阅以下代码示例：

    @Test
    public void testSingle4312ToWGS84() throws FactoryException, TransformException {
        double latitude = 29.0;
        double longitude = -99.0;

        CoordinateReferenceSystem sourceCrs = CRS.decode("EPSG:4312");
        CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");
        MathTransform engine = CRS.findMathTransform(sourceCrs, targetCrs, false);

        DirectPosition source = new DirectPosition2D(sourceCrs, longitude, latitude);
        DirectPosition target = new DirectPosition2D(targetCrs);
        engine.transform(source, target);

        System.out.println("x,y=" + target.getCoordinate()[0] + "," + target.getCoordinate()[1]);
    }

运行上述代码时，我得到以下结果：

x,y=-81.00483829765083,-150.99434404015307

这是错误的，因为 -150 不是有效的纬度。

当我创建 source 对象时，我用经度初始化 X，用纬度初始化 Y，遵循 X 是经度，Y 是纬度的基本原理（根据地理社区）。 但是，如果我反转它们：

DirectPosition source = new DirectPosition2D(sourceCrs, latitude, longitude);

我现在得到的结果看起来更真实：

x,y=29.003866857343915,-98.99222530180596

当然，我现在必须将 X 解释为纬度，将 Y 解释为经度。

那么，问题 #1，我们是否应该采用逆理性，将 X 视为纬度，将 Y 视为经度？

接下来，我想验证结果。首先，尝试使用 Oracle SQL SDO 函数转换同一点：

select sdo_cs.transform (
MDSYS.SDO_GEOMETRY(2001, 4312, MDSYS.SDO_POINT_TYPE(-99, 29, NULL), NULL, NULL),
8307
) from dual;

提供以下结果：

MDSYS.SDO_GEOMETRY(2001, 8307, MDSYS.SDO_POINT_TYPE(-98.9924748682774, 29.0036029261273, NULL), NULL, NULL)

转换后的点与 GeoTools 生成的点相似，但相差约 40 米。

接下来，我尝试了website 提供的转换，它提供了与 Oracle 相同的结果。 我还使用了 Luciad AIXM Viewer，这是一个 GML 的图形表示工具，它也会在 Oracle 计算的同一位置显示我的点。

所以，我有 3 个独立的工具提供相同的结果点，该结果点与 GeoTools 计算的结果点相距 40 米。 问题 #2，为什么 Geotools 会失败？ Geotools 对于 CRS 转换是否可靠？

Answer 1

对于轴顺序的第一个问题，您遇到了一个常见的初学者问题，即假设您知道投影的轴顺序是什么，并且在 ESPG:4326 的情况下它是固定的。

例如下面的代码：

CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");

System.out.println("EPSG:4326 - " + CRS.getAxisOrder(targetCrs));
System.out.println("WGS84 - " + CRS.getAxisOrder(DefaultGeographicCRS.WGS84));

给出这个输出：

EPSG:4326 - NORTH_EAST
WGS84 - EAST_NORTH

但是如果你添加Hints.putSystemDefault(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);这一行，输出变成：

EPSG:4326 - EAST_NORTH
WGS84 - EAST_NORTH

因此，一般来说，您永远不应该在代码中依赖“已知”的轴顺序，而是使用类似的东西：

double latitude = 29.0;
double longitude = -99.0;

CoordinateReferenceSystem sourceCrs = CRS.decode("EPSG:4312");
CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");

MathTransform engine = CRS.findMathTransform(sourceCrs, targetCrs, false);

DirectPosition2D source;
if (CRS.getAxisOrder(sourceCrs).equals(org.geotools.referencing.CRS.AxisOrder.EAST_NORTH)) {
  source = new DirectPosition2D(sourceCrs, longitude, latitude);
} else {
  source = new DirectPosition2D(sourceCrs, latitude, longitude);
}
DirectPosition target = new DirectPosition2D(targetCrs);
engine.transform(source, target);

if (CRS.getAxisOrder(targetCrs).equals(org.geotools.referencing.CRS.AxisOrder.EAST_NORTH)) {
  System.out.println("lon,lat=" + target.getCoordinate()[0] + "," + target.getCoordinate()[1]);
} else {
  System.out.println("lon,lat=" + target.getCoordinate()[1] + "," + target.getCoordinate()[0]);
}

至于转换的准确性取决于您导入的引用模块以及 Oracle 使用的 CRS 定义。加载gt-epsg-hsql 模块后，我看到ToWGS84[601.705, 84.263, 485.227, 4.7354, -1.3145, -5.393, -2.3887] 矩阵与EPSG registry 中提供的最准确的变换相同。如果有可用的 NTv2 转换，那么adding that to your project 应该会提高精度。