第07章 - 坐标系与空间参考

第07章 - 坐标系与空间参考

7.1 空间参考概述

7.1.1 SpatialReference 类

SpatialReference 是 geometry-api-java 中管理坐标系的核心类：

public abstract class SpatialReference {
    // 通过 WKID 创建
    public static SpatialReference create(int wkid);
    
    // 通过 WKT 创建
    public static SpatialReference create(String wktext);
    
    // 从 JSON 解析
    public static SpatialReference fromJson(String json);
    
    // 获取 WKID
    public abstract int getID();
    
    // 获取 WKT
    public abstract String getText();
    
    // 获取容差
    public double getTolerance();
    public abstract double getTolerance(int semantics);
}

7.1.2 创建空间参考

// 使用 WKID（推荐）
SpatialReference wgs84 = SpatialReference.create(4326);      // WGS84
SpatialReference webMercator = SpatialReference.create(3857); // Web Mercator
SpatialReference cgcs2000 = SpatialReference.create(4490);    // CGCS2000

// 使用 WKT
String wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",...]]";
SpatialReference srFromWkt = SpatialReference.create(wkt);

// 从 JSON 解析
String json = "{\"wkid\": 4326}";
SpatialReference srFromJson = SpatialReference.fromJson(json);

7.1.3 常用 WKID

// 地理坐标系
int WGS84 = 4326;        // WGS 1984
int CGCS2000 = 4490;     // China Geodetic Coordinate System 2000
int GCJ02 = -1;          // 不是标准坐标系，需要特殊处理

// 投影坐标系
int WEB_MERCATOR = 3857;  // Web Mercator (Auxiliary Sphere)
int WEB_MERCATOR_102100 = 102100; // Web Mercator (旧版)

// UTM 分带
int UTM_50N = 32650;     // UTM Zone 50N
int UTM_51N = 32651;     // UTM Zone 51N

// 高斯-克吕格投影（中国）
// CGCS2000 3度带
int CGCS2000_3DEG_ZONE_37 = 4526;  // 中央经线 111°
int CGCS2000_3DEG_ZONE_38 = 4527;  // 中央经线 114°
int CGCS2000_3DEG_ZONE_39 = 4528;  // 中央经线 117°

7.2 容差与精度

7.2.1 容差的概念

容差是空间参考的关键属性，定义了两个点被认为相同的最小距离：

SpatialReference sr = SpatialReference.create(4326);
double tolerance = sr.getTolerance();
System.out.println("容差: " + tolerance);  // 通常是 1e-9

// 不同维度的容差
double xyTolerance = sr.getTolerance(VertexDescription.Semantics.POSITION);
double zTolerance = sr.getTolerance(VertexDescription.Semantics.Z);

7.2.2 容差的作用

容差在以下操作中起关键作用：

1. 拓扑操作：合并、相交、差集等
2. 关系判断：相等、相交等
3. 简化操作：顶点合并

// 两个非常接近的点
Point p1 = new Point(116.4, 39.9);
Point p2 = new Point(116.400000001, 39.9);

SpatialReference sr = SpatialReference.create(4326);

// 如果差值小于容差，被视为相等
boolean equals = GeometryEngine.equals(p1, p2, sr);

7.2.3 容差与坐标系的关系

// 地理坐标系（度）的容差约 1e-9 度
// 对应约 0.1 毫米

// 投影坐标系（米）的容差需要相应调整
// 0.001 米 = 1 毫米

// geometry-api-java 根据坐标系自动设置合适的容差

7.3 投影转换

7.3.1 OperatorProject

geometry-api-java 本身不提供完整的投影转换功能，OperatorProject 需要在子类中实现：

// 注意：OperatorProject 需要外部实现
// geometry-api-java 只提供接口

// 通常需要配合 proj4j 或其他投影库使用

7.3.2 使用 proj4j 实现投影

// 添加 proj4j 依赖
// <dependency>
//     <groupId>org.locationtech.proj4j</groupId>
//     <artifactId>proj4j</artifactId>
//     <version>1.1.5</version>
// </dependency>

import org.locationtech.proj4j.*;

public class ProjectionHelper {
    private static final CRSFactory crsFactory = new CRSFactory();
    private static final CoordinateTransformFactory ctFactory = 
        new CoordinateTransformFactory();
    
    /**
     * 将几何从一个坐标系转换到另一个
     */
    public static Geometry project(Geometry geometry, 
            String fromCrsCode, String toCrsCode) {
        
        CoordinateReferenceSystem fromCrs = crsFactory.createFromName(fromCrsCode);
        CoordinateReferenceSystem toCrs = crsFactory.createFromName(toCrsCode);
        
        CoordinateTransform transform = ctFactory.createTransform(fromCrs, toCrs);
        
        return projectGeometry(geometry, transform);
    }
    
    private static Geometry projectGeometry(Geometry geometry, 
            CoordinateTransform transform) {
        
        switch (geometry.getType()) {
            case Point:
                return projectPoint((Point) geometry, transform);
            case Polyline:
                return projectPolyline((Polyline) geometry, transform);
            case Polygon:
                return projectPolygon((Polygon) geometry, transform);
            default:
                throw new UnsupportedOperationException();
        }
    }
    
    private static Point projectPoint(Point point, CoordinateTransform transform) {
        ProjCoordinate src = new ProjCoordinate(point.getX(), point.getY());
        ProjCoordinate dst = new ProjCoordinate();
        transform.transform(src, dst);
        return new Point(dst.x, dst.y);
    }
    
    // ... 其他几何类型的投影实现
}

7.3.3 常见投影转换场景

// WGS84 → Web Mercator（地图显示）
Geometry webMercatorGeom = ProjectionHelper.project(
    wgs84Geometry, "EPSG:4326", "EPSG:3857");

// Web Mercator → WGS84（GPS 定位）
Geometry wgs84Geom = ProjectionHelper.project(
    webMercatorGeometry, "EPSG:3857", "EPSG:4326");

// WGS84 → UTM（距离计算）
// 根据经度选择合适的 UTM 分带
int utmZone = (int) ((longitude + 180) / 6) + 1;
String utmCode = "EPSG:" + (32600 + utmZone);  // 北半球
Geometry utmGeom = ProjectionHelper.project(geometry, "EPSG:4326", utmCode);

7.4 测地计算

7.4.1 测地距离

// 计算 WGS84 椭球体上两点间的测地距离
Point beijing = new Point(116.397428, 39.90923);
Point shanghai = new Point(121.469170, 31.224361);

double distanceMeters = GeometryEngine.geodesicDistanceOnWGS84(beijing, shanghai);
System.out.println("距离: " + (distanceMeters / 1000) + " 公里");
// 约 1068 公里

7.4.2 测地长度

// 计算线的测地长度
Polyline route = new Polyline();
route.startPath(116.397428, 39.90923);   // 北京
route.lineTo(117.200983, 39.084158);     // 天津
route.lineTo(121.469170, 31.224361);     // 上海

SpatialReference sr = SpatialReference.create(4326);

OperatorGeodeticLength op = (OperatorGeodeticLength) OperatorFactoryLocal
    .getInstance().getOperator(Operator.Type.GeodeticLength);

double lengthMeters = op.execute(route, sr, 
    GeodeticCurveType.Geodesic, null);

System.out.println("路线长度: " + (lengthMeters / 1000) + " 公里");

7.4.3 测地面积

// 计算多边形的测地面积
Polygon polygon = createPolygon();
SpatialReference sr = SpatialReference.create(4326);

OperatorGeodeticArea op = (OperatorGeodeticArea) OperatorFactoryLocal
    .getInstance().getOperator(Operator.Type.GeodeticArea);

double areaSquareMeters = op.execute(polygon, sr, 
    GeodeticCurveType.Geodesic, null);

System.out.println("面积: " + (areaSquareMeters / 1000000) + " 平方公里");

7.4.4 测地缓冲区

// 创建测地缓冲区（考虑地球曲率）
Point point = new Point(116.4, 39.9);
SpatialReference sr = SpatialReference.create(4326);

OperatorGeodesicBuffer op = (OperatorGeodesicBuffer) OperatorFactoryLocal
    .getInstance().getOperator(Operator.Type.GeodesicBuffer);

double[] distances = {1000.0};  // 1000 米
double maxDeviation = 0.0;  // 使用默认偏差

GeometryCursor result = op.execute(
    new SimpleGeometryCursor(point),
    sr,
    GeodeticCurveType.Geodesic,
    distances,
    maxDeviation,
    false,  // 是否合并
    true,   // 是否保持输入方向
    null
);

Polygon buffer = (Polygon) result.next();

7.4.5 测地加密

// 按测地距离加密线（考虑地球曲率）
Polyline sparseLine = createLongPolyline();
SpatialReference sr = SpatialReference.create(4326);

OperatorGeodeticDensifyByLength op = 
    (OperatorGeodeticDensifyByLength) OperatorFactoryLocal
    .getInstance().getOperator(Operator.Type.GeodeticDensifyByLength);

// 每 10 公里加密一个点
double maxSegmentLength = 10000.0;

GeometryCursor result = op.execute(
    new SimpleGeometryCursor(sparseLine),
    sr,
    maxSegmentLength,
    GeodeticCurveType.Geodesic,
    null
);

Polyline densified = (Polyline) result.next();

7.5 坐标变换

7.5.1 Transformation2D

// 2D 仿射变换
Transformation2D transform = new Transformation2D();

// 平移
transform.setShift(100, 200);

// 旋转（弧度）
double angle = Math.PI / 4;  // 45度
transform.setRotate(Math.cos(angle), Math.sin(angle));

// 缩放
transform.setScale(2.0, 2.0);

// 应用变换
Point point = new Point(10, 10);
point.applyTransformation(transform);

7.5.2 Transformation3D

// 3D 变换
Transformation3D transform = new Transformation3D();

// 设置变换矩阵
// [xx xy xz xd]
// [yx yy yz yd]
// [zx zy zz zd]

Point3D point = new Point3D();
point.setCoords(10, 10, 10);
Point3D transformed = transform.transform(point);

7.5.3 批量坐标变换

/**
 * 批量应用坐标变换
 */
public void transformGeometries(List<Geometry> geometries, 
        Transformation2D transform) {
    for (Geometry geom : geometries) {
        geom.applyTransformation(transform);
    }
}

7.6 MapGeometry 类

7.6.1 MapGeometry 概述

MapGeometry 封装了几何对象及其空间参考：

public class MapGeometry {
    private Geometry m_geometry;
    private SpatialReference m_spatialReference;
    
    public Geometry getGeometry();
    public SpatialReference getSpatialReference();
}

7.6.2 使用 MapGeometry

// 从 JSON 解析得到 MapGeometry
String json = "{\"x\":116.4,\"y\":39.9,\"spatialReference\":{\"wkid\":4326}}";
MapGeometry mapGeom = GeometryEngine.jsonToGeometry(json);

Geometry geometry = mapGeom.getGeometry();
SpatialReference sr = mapGeom.getSpatialReference();

// 创建 MapGeometry
Polygon polygon = createPolygon();
SpatialReference sr2 = SpatialReference.create(4326);
MapGeometry mapGeom2 = new MapGeometry(polygon, sr2);

7.7 坐标系处理实战

7.7.1 自动识别坐标系

/**
 * 根据坐标范围推断坐标系
 */
public SpatialReference guessCoordinateSystem(Geometry geometry) {
    Envelope2D env = new Envelope2D();
    geometry.queryEnvelope2D(env);
    
    // 检查是否在经纬度范围内
    if (env.xmin >= -180 && env.xmax <= 180 &&
        env.ymin >= -90 && env.ymax <= 90) {
        
        // 可能是 WGS84 或 CGCS2000
        if (isInChina(env)) {
            return SpatialReference.create(4490); // CGCS2000
        } else {
            return SpatialReference.create(4326); // WGS84
        }
    }
    
    // 检查是否在 Web Mercator 范围内
    if (Math.abs(env.xmax) > 180 || Math.abs(env.ymax) > 90) {
        return SpatialReference.create(3857);
    }
    
    // 无法确定
    return null;
}

private boolean isInChina(Envelope2D env) {
    return env.xmin >= 73 && env.xmax <= 136 &&
           env.ymin >= 3 && env.ymax <= 54;
}

7.7.2 计算 UTM 分带

/**
 * 根据经度计算 UTM 分带
 */
public int getUtmZone(double longitude) {
    return (int) ((longitude + 180) / 6) + 1;
}

/**
 * 获取 UTM 分带的 EPSG 代码
 */
public int getUtmEpsgCode(double longitude, double latitude) {
    int zone = getUtmZone(longitude);
    if (latitude >= 0) {
        return 32600 + zone;  // 北半球
    } else {
        return 32700 + zone;  // 南半球
    }
}

/**
 * 获取适合几何的 UTM 坐标系
 */
public SpatialReference getAppropriateUtm(Geometry geometry) {
    Envelope2D env = new Envelope2D();
    geometry.queryEnvelope2D(env);
    
    double centerLon = env.getCenterX();
    double centerLat = env.getCenterY();
    
    int epsg = getUtmEpsgCode(centerLon, centerLat);
    return SpatialReference.create(epsg);
}

7.7.3 度分秒转换

/**
 * 度分秒转十进制度
 */
public double dmsToDegrees(int degrees, int minutes, double seconds) {
    double sign = degrees >= 0 ? 1 : -1;
    return sign * (Math.abs(degrees) + minutes / 60.0 + seconds / 3600.0);
}

/**
 * 十进制度转度分秒
 */
public String degreesToDms(double degrees) {
    int sign = degrees >= 0 ? 1 : -1;
    degrees = Math.abs(degrees);
    
    int d = (int) degrees;
    double remaining = (degrees - d) * 60;
    int m = (int) remaining;
    double s = (remaining - m) * 60;
    
    String direction = sign >= 0 ? "" : "-";
    return String.format("%s%d°%d'%.2f\"", direction, d, m, s);
}

7.7.4 中国 CGCS2000 分带

/**
 * 获取 CGCS2000 高斯-克吕格 3 度带
 */
public SpatialReference getCgcs2000GK3DegZone(double longitude) {
    // 3度带号 = (经度 - 1.5) / 3 + 1
    int zone = (int) ((longitude - 1.5) / 3) + 1;
    
    // CGCS2000 3度带 EPSG: 4520 + 带号
    int epsg = 4520 + zone;
    return SpatialReference.create(epsg);
}

/**
 * 获取 CGCS2000 高斯-克吕格 6 度带
 */
public SpatialReference getCgcs2000GK6DegZone(double longitude) {
    // 6度带号 = (经度 + 6) / 6
    int zone = (int) ((longitude + 6) / 6);
    
    // CGCS2000 6度带 EPSG: 4502 + 带号
    int epsg = 4502 + zone;
    return SpatialReference.create(epsg);
}

7.8 本章小结

本章详细介绍了 geometry-api-java 的坐标系与空间参考：

1. SpatialReference：空间参考类的使用
2. 容差与精度：理解和应用容差
3. 投影转换：使用 proj4j 实现投影
4. 测地计算：距离、长度、面积、缓冲区
5. 坐标变换：2D/3D 仿射变换
6. MapGeometry：几何与空间参考的封装
7. 实战应用：UTM 分带、度分秒转换等

关键要点：

• geometry-api-java 提供空间参考管理，但不包含完整投影库
• 测地计算基于 WGS84 椭球体
• 容差影响拓扑和关系操作的结果
• 需要配合 proj4j 等库实现完整投影转换

---

← 上一章：数据格式转换 | 下一章：OGC 兼容开发 →