如何使用 GeoTools MathTransform？

如何解决如何使用 GeoTools MathTransform？

我需要使用 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 不是一个有效的纬度。

当我创建源对象时，我用经度初始化了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,8307
) from dual;

提供以下结果：

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

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

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

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

解决方法

对于轴顺序的第一个问题，您已经陷入了一个常见的初学者问题，即假设您知道投影的轴顺序是什么，并且在 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 模块后，我看到与 EPSG registry 中提供的最准确变换相同的 ToWGS84[601.705,84.263,485.227,4.7354,-1.3145,-5.393,-2.3887] 矩阵。如果有可用的 NTv2 转换，那么 adding that to your project 应该会提高精度。