GeoTools坐标转换(scala)

发表于 更新于 分类于 Waline： 阅读次数： 本文字数： 617 阅读时长 ≈ 2 分钟

今天在使用 GeoTools 做坐标转换的时候，发现转出来的坐标与预期不符。经过查询 GeoTools 文档，发现需要通过设置 Axis Order 来指定经纬度哪个在前哪个在后，在此记录一下。

参考: Axis Order — GeoTools 27-SNAPSHOT User Guide

问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import org.geotools.geometry.jts.JTSFactoryFinder
import org.geotools.geometry.jts.JTS
import org.geotools.referencing.CRS
import org.locationtech.jts.geom.Coordinate
import org.locationtech.jts.geom.GeometryFactory
import org.locationtech.jts.geom.Point

// 创建个点用于转换
val geometryFactory: GeometryFactory = JTSFactoryFinder.getGeometryFactory
val coord: Coordinate = new Coordinate(121, 32)
val point: Point = geometryFactory.createPoint(coord)

val mathTransform = CRS.findMathTransform(CRS.decode("EPSG:4490"), CRS.decode("EPSG:4528"))
val newPoint: Point = JTS.transform(point, mathTransform).asInstanceOf[Point]

这段代码乍一看没啥问题，就是转出来的 Point 坐标是错误的。

经过查询 GeoTools 文档，可以知道，由于历史原因，GeoTools 无法自动判断何时返回与 EPSG 库一致的 CRS，以及何时返回带有轴序配置的 CRS。所以 GeoTools 决定无论轴序是什么样的，都认为是和 EPSG 一致的 CRS，也就是 x 为纬度，y 为经度。如果需要调整数据的轴序，可以由用户自行规定 CRS。

解决方案

按照我们一般的习惯，想要经度在前，纬度在后。

为了帮助旧应用程序实现过渡，可以通过设置系统属性org.geotools.referencing.forceXY为true，此时，会将提示值FORCE_LONGITUDE_FIRST_AXIS_ORDER设置为true，将会强制认为轴序为经度在前，纬度在后。

1
System.setProperty("org.geotools.referencing.forceXY", "true")

也可以在代码里手动指定，代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import com.supermap.bdt.base.algorithm.CRSTransform
import org.geotools.factory.Hints
import org.geotools.geometry.jts.JTSFactoryFinder
import org.geotools.geometry.jts.JTS
import org.geotools.referencing.CRS
import org.geotools.referencing.ReferencingFactoryFinder
import org.locationtech.jts.geom.Coordinate
import org.locationtech.jts.geom.GeometryFactory
import org.locationtech.jts.geom.Point

// 创建个点用于转换
val geometryFactory: GeometryFactory = JTSFactoryFinder.getGeometryFactory
val coord: Coordinate = new Coordinate(121, 32)
val point: Point = geometryFactory.createPoint(coord)

val hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, true)
val factory = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints)
val mathTransform = CRS.findMathTransform(factory.createCoordinateReferenceSystem("EPSG:4490"), factory.createCoordinateReferenceSystem("EPSG:4528"))
val newPoint: Point = JTS.transform(point, mathTransform).asInstanceOf[Point]

通过手动指定了轴序，转出来的结果就是正确的了。

SuperMap iObject for Spark

我是在 Spark 环境下使用坐标转换，超图的 SuperMap iObject for Spark 中也引了一下这个方法。这里改成超图 iObject 里的方法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import com.supermap.bdt.base.algorithm.CRSTransform
import org.geotools.factory.Hints
import org.geotools.geometry.jts.JTSFactoryFinder
import org.geotools.referencing.CRS
import org.geotools.referencing.ReferencingFactoryFinder
import org.locationtech.jts.geom.Coordinate
import org.locationtech.jts.geom.GeometryFactory
import org.locationtech.jts.geom.Point

// 创建个点用于转换
val geometryFactory: GeometryFactory = JTSFactoryFinder.getGeometryFactory
val coord: Coordinate = new Coordinate(121, 32)
val point: Point = geometryFactory.createPoint(coord)

val hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, true)
val factory = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints)
val mathTransform = CRSTransform.findMathTransform(factory.createCoordinateReferenceSystem("EPSG:4490"), factory.createCoordinateReferenceSystem("EPSG:4528"))
val newPoint: Point = CRSTransform.transform(point, mathTransform).asInstanceOf[Point]

直接转换 FeatureRDD，这边不写 FeatureRDD 怎么构建了。

1
2
3
4
5
6
7
8
9
10
11
12
import com.supermap.bdt.base.algorithm.CRSTransform
import com.supermap.bdt.FeatureRDD
import org.geotools.factory.Hints
import org.geotools.geometry.jts.JTSFactoryFinder
import org.geotools.referencing.ReferencingFactoryFinder

// 此处仅为示例
val featureRDD: FeatureRDD = null

val hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, true)
val factory = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints)
val newFeatureRDD: FeatureRDD = CRSTransform.transform(featureRDD, factory.createCoordinateReferenceSystem("EPSG:4528"))