第10章 - 坐标参考系统与投影转换
第10章 - 坐标参考系统与投影转换
10.1 CRS 基础概念
10.1.1 坐标参考系统
坐标参考系统（Coordinate Reference System, CRS）定义了如何将地球上的位置映射到平面或三维坐标。
┌─────────────────────────────────────────────────────────────────────┐
│ 坐标参考系统分类 │
├─────────────────────────────────────────────────────────────────────┤
│ │
│ 地理坐标系（Geographic CRS） │
│ ├── 使用经纬度表示位置 │
│ ├── 单位：度、分、秒 │
│ └── 例：WGS84 (EPSG:4326)、CGCS2000 (EPSG:4490) │
│ │
│ 投影坐标系（Projected CRS） │
│ ├── 将球面投影到平面 │
│ ├── 单位：米、英尺等 │
│ └── 例：Web Mercator (EPSG:3857)、UTM │
│ │
│ 地心坐标系（Geocentric CRS） │
│ ├── 以地心为原点的三维坐标 │
│ └── 用于卫星定位等 │
│ │
└─────────────────────────────────────────────────────────────────────┘
10.1.2 常用 EPSG 代码
| EPSG 代码 | 名称 | 类型 | 说明 |
|---|---|---|---|
| 4326 | WGS 84 | 地理 | GPS 全球坐标系 |
| 3857 | Web Mercator | 投影 | Web 地图标准 |
| 4490 | CGCS2000 | 地理 | 中国大地坐标系 |
| 4479 | CGCS2000 3D | 地理 | 中国大地坐标系（3D） |
| 32650 | UTM Zone 50N | 投影 | UTM 50N（北京地区） |
| 4547 | CGCS2000 / 3-degree Gauss-Kruger zone 39 | 投影 | 中国高斯投影 |
10.2 CRS 操作
10.2.1 获取 CRS
import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
import org.geotools.referencing.CRS;
public class CRSOperations {
public static void getCRS() throws Exception {
// 1. 通过 EPSG 代码获取
CoordinateReferenceSystem wgs84 = CRS.decode("EPSG:4326");
System.out.println("WGS84: " + wgs84.getName());
// 2. 强制经度优先（lon/lat）
CoordinateReferenceSystem wgs84LonFirst = CRS.decode("EPSG:4326", true);
// 3. 通过 WKT 获取
String wkt = "GEOGCS[\"WGS 84\"," +
"DATUM[\"WGS_1984\"," +
"SPHEROID[\"WGS 84\",6378137,298.257223563]]," +
"PRIMEM[\"Greenwich\",0]," +
"UNIT[\"degree\",0.0174532925199433]]";
CoordinateReferenceSystem fromWKT = CRS.parseWKT(wkt);
// 4. 从 Geometry 获取
// Point point = ...; // 设置了 CRS 的几何
// CoordinateReferenceSystem geomCRS =
// CRS.getHorizontalCRS(JTS.getCRS(point));
// 5. 预定义 CRS
CoordinateReferenceSystem defaultCRS =
org.geotools.referencing.crs.DefaultGeographicCRS.WGS84;
}
public static void crsInfo(CoordinateReferenceSystem crs) {
System.out.println("=== CRS 信息 ===");
System.out.println("名称: " + crs.getName());
System.out.println("标识符: " + crs.getIdentifiers());
System.out.println("域: " + crs.getDomainOfValidity());
System.out.println("范围: " + crs.getScope());
// 获取 EPSG 代码
try {
Integer epsgCode = CRS.lookupEpsgCode(crs, true);
System.out.println("EPSG 代码: " + epsgCode);
} catch (Exception e) {
System.out.println("无法获取 EPSG 代码");
}
// WKT 表示
System.out.println("\nWKT:");
System.out.println(crs.toWKT());
}
}
10.2.2 CRS 比较
public class CRSComparison {
public static void compareCRS() throws Exception {
CoordinateReferenceSystem crs1 = CRS.decode("EPSG:4326");
CoordinateReferenceSystem crs2 = CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem crs3 = CRS.decode("EPSG:3857");
// 严格比较
boolean strictEqual = CRS.equalsIgnoreMetadata(crs1, crs2);
System.out.println("4326 vs 4326(lonFirst): " + strictEqual);
// 检查是否等效
boolean equivalent = CRS.findMathTransform(crs1, crs2, true).isIdentity();
System.out.println("转换是否为恒等: " + equivalent);
// 检查轴顺序
System.out.println("CRS1 轴顺序: " +
CRS.getAxisOrder(crs1)); // LAT_LON 或 LON_LAT
System.out.println("CRS2 轴顺序: " +
CRS.getAxisOrder(crs2));
}
}
10.3 坐标转换
10.3.1 基本转换
import org.geotools.api.referencing.operation.MathTransform;
import org.geotools.geometry.jts.JTS;
public class CoordinateTransform {
public static void basicTransform() throws Exception {
// 源和目标 CRS
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");
// 创建转换器
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
// 转换坐标数组
double[] srcPt = {116.4074, 39.9042}; // 北京 (lon, lat)
double[] dstPt = new double[2];
transform.transform(srcPt, 0, dstPt, 0, 1);
System.out.printf("WGS84: (%.4f, %.4f)%n", srcPt[0], srcPt[1]);
System.out.printf("Web Mercator: (%.2f, %.2f)%n", dstPt[0], dstPt[1]);
// 反向转换
MathTransform inverse = transform.inverse();
double[] backPt = new double[2];
inverse.transform(dstPt, 0, backPt, 0, 1);
System.out.printf("反向: (%.4f, %.4f)%n", backPt[0], backPt[1]);
}
public static void transformGeometry() throws Exception {
GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();
// 创建 WGS84 点
Point wgs84Point = gf.createPoint(new Coordinate(116.4074, 39.9042));
// 获取转换器
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
// 转换几何
Geometry webMercatorPoint = JTS.transform(wgs84Point, transform);
System.out.println("转换后: " + webMercatorPoint);
// 转换坐标
Coordinate srcCoord = new Coordinate(116.4074, 39.9042);
Coordinate dstCoord = JTS.transform(srcCoord, null, transform);
System.out.println("坐标转换: " + dstCoord);
}
}
10.3.2 转换复杂几何
public class GeometryReproject {
public static Geometry reproject(Geometry geometry,
CoordinateReferenceSystem sourceCRS,
CoordinateReferenceSystem targetCRS) throws Exception {
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
return JTS.transform(geometry, transform);
}
public static SimpleFeatureCollection reprojectCollection(
SimpleFeatureCollection fc,
CoordinateReferenceSystem targetCRS) throws Exception {
SimpleFeatureType sourceType = fc.getSchema();
CoordinateReferenceSystem sourceCRS = sourceType.getCoordinateReferenceSystem();
// 创建目标类型
SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
typeBuilder.init(sourceType);
typeBuilder.setCRS(targetCRS);
SimpleFeatureType targetType = typeBuilder.buildFeatureType();
// 获取转换器
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
// 转换要素
DefaultFeatureCollection result = new DefaultFeatureCollection();
SimpleFeatureBuilder builder = new SimpleFeatureBuilder(targetType);
try (SimpleFeatureIterator iter = fc.features()) {
while (iter.hasNext()) {
SimpleFeature feature = iter.next();
for (int i = 0; i < feature.getAttributeCount(); i++) {
Object attr = feature.getAttribute(i);
if (attr instanceof Geometry) {
attr = JTS.transform((Geometry) attr, transform);
}
builder.add(attr);
}
result.add(builder.buildFeature(feature.getID()));
}
}
return result;
}
}
10.3.3 测地线计算
import org.geotools.referencing.GeodeticCalculator;
public class GeodeticOperations {
public static void calculateDistance() throws Exception {
GeodeticCalculator calculator = new GeodeticCalculator(
CRS.decode("EPSG:4326", true)
);
// 设置起点（北京）
calculator.setStartingGeographicPoint(116.4074, 39.9042);
// 设置终点（上海）
calculator.setDestinationGeographicPoint(121.4737, 31.2304);
// 计算测地线距离
double distance = calculator.getOrthodromicDistance();
System.out.printf("北京到上海距离: %.2f 公里%n", distance / 1000);
// 获取方位角
double azimuth = calculator.getAzimuth();
System.out.printf("方位角: %.2f 度%n", azimuth);
}
public static void calculateDestination() throws Exception {
GeodeticCalculator calculator = new GeodeticCalculator(
CRS.decode("EPSG:4326", true)
);
// 起点
calculator.setStartingGeographicPoint(116.4074, 39.9042);
// 设置方向和距离
calculator.setDirection(90, 100000); // 向东 100 公里
// 获取目的地
java.awt.geom.Point2D dest = calculator.getDestinationGeographicPoint();
System.out.printf("目的地: (%.4f, %.4f)%n", dest.getX(), dest.getY());
}
public static double calculateArea(Polygon polygon,
CoordinateReferenceSystem crs) throws Exception {
// 对于地理坐标系，需要投影到适当的投影坐标系计算面积
if (crs.equals(CRS.decode("EPSG:4326"))) {
// 获取多边形中心
Point centroid = polygon.getCentroid();
// 找到适当的 UTM 带
int zone = (int) Math.floor((centroid.getX() + 180) / 6) + 1;
String utmCode = centroid.getY() >= 0 ?
"EPSG:326" + String.format("%02d", zone) : // 北半球
"EPSG:327" + String.format("%02d", zone); // 南半球
CoordinateReferenceSystem utmCRS = CRS.decode(utmCode);
MathTransform transform = CRS.findMathTransform(crs, utmCRS, true);
Geometry projected = JTS.transform(polygon, transform);
return projected.getArea(); // 平方米
}
return polygon.getArea();
}
}
10.4 自定义 CRS
10.4.1 从 WKT 创建
public class CustomCRS {
public static CoordinateReferenceSystem createFromWKT() throws Exception {
// 自定义横轴墨卡托投影
String wkt = "PROJCS[\"Custom_Transverse_Mercator\"," +
"GEOGCS[\"GCS_WGS_1984\"," +
"DATUM[\"D_WGS_1984\"," +
"SPHEROID[\"WGS_1984\",6378137.0,298.257223563]]," +
"PRIMEM[\"Greenwich\",0.0]," +
"UNIT[\"Degree\",0.0174532925199433]]," +
"PROJECTION[\"Transverse_Mercator\"]," +
"PARAMETER[\"False_Easting\",500000.0]," +
"PARAMETER[\"False_Northing\",0.0]," +
"PARAMETER[\"Central_Meridian\",117.0]," +
"PARAMETER[\"Scale_Factor\",1.0]," +
"PARAMETER[\"Latitude_Of_Origin\",0.0]," +
"UNIT[\"Meter\",1.0]]";
return CRS.parseWKT(wkt);
}
public static CoordinateReferenceSystem createGaussKruger(int zone) throws Exception {
// 中国高斯-克吕格投影
double centralMeridian = zone * 3; // 3度带中央经线
String wkt = String.format(
"PROJCS[\"CGCS2000 / 3-degree Gauss-Kruger zone %d\"," +
"GEOGCS[\"China Geodetic Coordinate System 2000\"," +
"DATUM[\"China_2000\"," +
"SPHEROID[\"CGCS2000\",6378137,298.257222101]]," +
"PRIMEM[\"Greenwich\",0]," +
"UNIT[\"degree\",0.0174532925199433]]," +
"PROJECTION[\"Transverse_Mercator\"]," +
"PARAMETER[\"latitude_of_origin\",0]," +
"PARAMETER[\"central_meridian\",%.1f]," +
"PARAMETER[\"scale_factor\",1]," +
"PARAMETER[\"false_easting\",%d]," +
"PARAMETER[\"false_northing\",0]," +
"UNIT[\"metre\",1]]",
zone, centralMeridian, zone * 1000000 + 500000
);
return CRS.parseWKT(wkt);
}
}
10.5 处理轴顺序
10.5.1 轴顺序问题
public class AxisOrderHandling {
public static void handleAxisOrder() throws Exception {
// EPSG:4326 标准顺序是 lat/lon，但很多软件使用 lon/lat
// 默认（lat/lon）
CoordinateReferenceSystem crs1 = CRS.decode("EPSG:4326");
System.out.println("默认轴顺序: " + CRS.getAxisOrder(crs1));
// 强制 lon/lat
CoordinateReferenceSystem crs2 = CRS.decode("EPSG:4326", true);
System.out.println("强制轴顺序: " + CRS.getAxisOrder(crs2));
// 全局设置（不推荐，会影响整个应用）
// System.setProperty("org.geotools.referencing.forceXY", "true");
// 使用 Hints
Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
CoordinateReferenceSystem crs3 = CRS.decode("EPSG:4326", hints);
}
public static void transformWithAxisOrder() throws Exception {
// 源数据是 lon/lat 顺序
double lon = 116.4074;
double lat = 39.9042;
// 确保使用正确的轴顺序
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
double[] srcPt = {lon, lat};
double[] dstPt = new double[2];
transform.transform(srcPt, 0, dstPt, 0, 1);
System.out.printf("转换结果: (%.2f, %.2f)%n", dstPt[0], dstPt[1]);
}
}
10.6 批量重投影
10.6.1 重投影 FeatureCollection
public class BatchReproject {
public static SimpleFeatureCollection reproject(
SimpleFeatureCollection source,
String targetEPSG) throws Exception {
SimpleFeatureType sourceType = source.getSchema();
CoordinateReferenceSystem sourceCRS = sourceType.getCoordinateReferenceSystem();
CoordinateReferenceSystem targetCRS = CRS.decode(targetEPSG, true);
if (CRS.equalsIgnoreMetadata(sourceCRS, targetCRS)) {
return source; // 无需转换
}
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
// 构建目标类型
SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
typeBuilder.init(sourceType);
typeBuilder.setCRS(targetCRS);
SimpleFeatureType targetType = typeBuilder.buildFeatureType();
DefaultFeatureCollection result = new DefaultFeatureCollection();
SimpleFeatureBuilder builder = new SimpleFeatureBuilder(targetType);
try (SimpleFeatureIterator iter = source.features()) {
while (iter.hasNext()) {
SimpleFeature feature = iter.next();
for (int i = 0; i < feature.getAttributeCount(); i++) {
Object attr = feature.getAttribute(i);
if (attr instanceof Geometry) {
attr = JTS.transform((Geometry) attr, transform);
}
builder.add(attr);
}
result.add(builder.buildFeature(feature.getID()));
}
}
return result;
}
}
10.7 本章小结
本章详细介绍了坐标参考系统和投影转换：
-
CRS 基础
- 地理坐标系
- 投影坐标系
- EPSG 代码
-
CRS 操作
- 获取 CRS
- CRS 比较
- 信息查询
-
坐标转换
- 基本转换
- 几何转换
- 测地线计算
-
高级特性
- 自定义 CRS
- 轴顺序处理
- 批量重投影
关键要点
- 注意轴顺序问题（lon/lat vs lat/lon）
- 使用
CRS.decode(code, true)强制经度优先 - 地理坐标系下的面积需要投影后计算
- 使用 GeodeticCalculator 进行精确的距离计算
浙公网安备 33010602011771号