第08章 - 坐标系管理与转换

第08章 - 坐标系管理与转换

8.1 坐标系基础知识

8.1.1 地理坐标系与投影坐标系

地理坐标系（Geographic Coordinate System, GCS）

• 使用经纬度表示位置
• 单位：度（degree）
• 例如：CGCS2000地理坐标系（WKID: 4490）

投影坐标系（Projected Coordinate System, PCS）

• 将球面坐标投影到平面
• 单位：米（meter）
• 例如：CGCS2000 / 3度分带高斯-克吕格投影（WKID: 4524-4554）

8.1.2 CGCS2000坐标系

中国采用的国家大地坐标系统（CGCS2000），OGU4J内置支持：

坐标系	WKID	说明
CGCS2000地理坐标系	4490	经纬度
CGCS2000 3度带 25N	4513	投影，中央经线75°
CGCS2000 3度带 26N	4514	投影，中央经线78°
...	...	...
CGCS2000 3度带 45N	4534	投影，中央经线135°

3度分带对照表（常用）

带号	中央经线	WKID	覆盖区域
35	105°	4524	云南、四川西部
36	108°	4525	陕西、重庆
37	111°	4526	山西、河南
38	114°	4527	河北、山东西部
39	117°	4528	北京、天津、山东
40	120°	4529	江苏、上海
41	123°	4530	辽宁、吉林

8.1.3 WKID（Well-Known ID）

WKID是坐标系的唯一标识符：

// 常用WKID
public class CommonWkids {
    // 地理坐标系
    public static final int WGS84 = 4326;           // GPS坐标系
    public static final int CGCS2000 = 4490;        // 中国国家坐标系
    
    // 投影坐标系
    public static final int WEB_MERCATOR = 3857;    // Web墨卡托
    public static final int CGCS2000_39N = 4528;    // CGCS2000 39带
    public static final int CGCS2000_40N = 4529;    // CGCS2000 40带
}

8.2 CrsUtil工具类

8.2.1 类概述

CrsUtil是OGU4J的坐标系管理核心类，位于com.znlgis.ogu4j.engine.util包：

/**
 * 坐标系工具类
 * 提供坐标转换、带号计算、CRS管理等功能
 */
public class CrsUtil {
    
    // 坐标转换
    public static String transform(String wkt, int sourceWkid, int targetWkid);
    public static Geometry transform(Geometry geom, int sourceWkid, int targetWkid);
    public static OguLayer reproject(OguLayer layer, int targetWkid);
    
    // 带号相关
    public static int getDh(Geometry geometry);
    public static int getDh(String wkt);
    public static int getDh(int projectedWkid);
    
    // WKID相关
    public static Integer getWkid(Geometry geometry);
    public static Integer getProjectedWkid(int dh);
    public static Integer getProjectedWkid(Geometry geometry);
    
    // CRS相关
    public static CoordinateReferenceSystem getCRS(int wkid);
    public static boolean isProjectedCRS(CoordinateReferenceSystem crs);
    public static double getTolerance(int wkid);
    public static Map<Integer, CoordinateReferenceSystem> getSupportedCRSList();
}

8.3 坐标转换

8.3.1 WKT字符串转换

import com.znlgis.ogu4j.engine.util.CrsUtil;

// 地理坐标 -> 投影坐标
String wktGeo = "POINT(116.4 39.9)";  // 北京，经纬度
int sourceWkid = 4490;  // CGCS2000地理坐标系
int targetWkid = 4528;  // CGCS2000 39带

String wktProj = CrsUtil.transform(wktGeo, sourceWkid, targetWkid);
System.out.println("投影坐标: " + wktProj);
// 输出类似: POINT(500000.0 4420000.0)

// 投影坐标 -> 地理坐标
String wktBack = CrsUtil.transform(wktProj, targetWkid, sourceWkid);
System.out.println("地理坐标: " + wktBack);
// 输出: POINT(116.4 39.9)

8.3.2 Geometry对象转换

import org.locationtech.jts.geom.Geometry;

// 创建Geometry
Geometry geomGeo = GeometryUtil.wkt2Geometry("POLYGON((116 39, 116 40, 117 40, 117 39, 116 39))");

// 转换
Geometry geomProj = CrsUtil.transform(geomGeo, 4490, 4528);
System.out.println("转换后: " + geomProj.toText());

// 面积对比
double areaGeo = GeometryUtil.area(geomGeo);    // 约1（平方度）
double areaProj = GeometryUtil.area(geomProj);  // 约xxx（平方米）
System.out.println("地理坐标面积（平方度）: " + areaGeo);
System.out.println("投影坐标面积（平方米）: " + areaProj);

8.3.3 图层整体转换

// 读取图层
OguLayer layer = OguLayerUtil.readLayer(
    DataFormatType.SHP,
    "D:/data/cities.shp",
    null, null, null,
    GisEngineType.GEOTOOLS
);

System.out.println("原坐标系: WKID " + layer.getWkid());

// 投影转换
OguLayer reprojected = CrsUtil.reproject(layer, 4528);

System.out.println("转换后: WKID " + reprojected.getWkid());

// 保存转换后的图层
OguLayerUtil.writeLayer(
    DataFormatType.SHP,
    reprojected,
    "D:/output/cities_projected.shp",
    null, null,
    GisEngineType.GEOTOOLS
);

8.4 带号计算

8.4.1 什么是带号

高斯-克吕格投影将地球按经度分成若干带，每带的中央经线不同：

• 3度分带：每带跨3个经度，中央经线为带号×3
• 6度分带：每带跨6个经度，中央经线为带号×6-3

带号计算公式（3度分带）：

带号 = (经度 + 1.5) / 3 的整数部分
中央经线 = 带号 × 3

8.4.2 从几何计算带号

// 从Geometry计算带号
Geometry geom = GeometryUtil.wkt2Geometry("POINT(116.4 39.9)");
int dh = CrsUtil.getDh(geom);
System.out.println("带号: " + dh);  // 39

// 从WKT计算带号
int dh2 = CrsUtil.getDh("POLYGON((116 39, 116 40, 117 40, 117 39, 116 39))");
System.out.println("带号: " + dh2);  // 39

// 从投影坐标系WKID计算带号
int dh3 = CrsUtil.getDh(4528);
System.out.println("带号: " + dh3);  // 39

8.4.3 带号与WKID转换

// 带号 -> 投影坐标系WKID
int wkid = CrsUtil.getProjectedWkid(39);
System.out.println("39带WKID: " + wkid);  // 4528

// 从几何获取对应的投影WKID
Geometry geom = GeometryUtil.wkt2Geometry("POINT(116.4 39.9)");
Integer projWkid = CrsUtil.getProjectedWkid(geom);
System.out.println("投影WKID: " + projWkid);  // 4528

8.5 WKID与CRS

8.5.1 获取CRS对象

import org.opengis.referencing.crs.CoordinateReferenceSystem;

// 获取CRS对象
CoordinateReferenceSystem crsGeo = CrsUtil.getCRS(4490);
CoordinateReferenceSystem crsProj = CrsUtil.getCRS(4528);

// 获取CRS信息
System.out.println("名称: " + crsGeo.getName());
System.out.println("WKT: " + crsGeo.toWKT());

8.5.2 判断坐标系类型

// 判断是否为投影坐标系
CoordinateReferenceSystem crsGeo = CrsUtil.getCRS(4490);
CoordinateReferenceSystem crsProj = CrsUtil.getCRS(4528);

boolean isProjectedGeo = CrsUtil.isProjectedCRS(crsGeo);   // false
boolean isProjectedProj = CrsUtil.isProjectedCRS(crsProj); // true

8.5.3 获取容差

// 获取坐标系对应的容差值
// 地理坐标系容差较小（度），投影坐标系容差较大（米）
double tolGeo = CrsUtil.getTolerance(4490);   // 例如: 0.00001
double tolProj = CrsUtil.getTolerance(4528);  // 例如: 0.001

8.5.4 获取支持的坐标系列表

// 获取OGU4J支持的所有坐标系
Map<Integer, CoordinateReferenceSystem> crsList = CrsUtil.getSupportedCRSList();

for (Map.Entry<Integer, CoordinateReferenceSystem> entry : crsList.entrySet()) {
    System.out.printf("WKID: %d, 名称: %s%n", 
        entry.getKey(), 
        entry.getValue().getName());
}

8.6 实战案例

8.6.1 智能投影转换

public class SmartReproject {
    
    /**
     * 智能投影转换
     * 自动计算目标投影坐标系
     */
    public static OguLayer smartReproject(OguLayer layer) {
        if (layer.getWkid() == null) {
            throw new IllegalArgumentException("图层没有坐标系信息");
        }
        
        // 获取CRS
        CoordinateReferenceSystem crs = CrsUtil.getCRS(layer.getWkid());
        
        // 如果已经是投影坐标系，直接返回
        if (CrsUtil.isProjectedCRS(crs)) {
            System.out.println("已是投影坐标系，无需转换");
            return layer;
        }
        
        // 计算合适的带号
        Geometry unionGeom = createUnionGeometry(layer);
        int dh = CrsUtil.getDh(unionGeom);
        
        // 获取目标WKID
        Integer targetWkid = CrsUtil.getProjectedWkid(dh);
        System.out.printf("自动选择投影: %d带, WKID %d%n", dh, targetWkid);
        
        // 执行转换
        return CrsUtil.reproject(layer, targetWkid);
    }
    
    /**
     * 创建图层的并集几何
     */
    private static Geometry createUnionGeometry(OguLayer layer) {
        Geometry union = null;
        for (OguFeature feature : layer.getFeatures()) {
            Geometry geom = GeometryUtil.wkt2Geometry(feature.getWkt());
            if (union == null) {
                union = geom;
            } else {
                union = GeometryUtil.union(union, geom);
            }
        }
        return union;
    }
}

8.6.2 坐标系验证

public class CrsValidator {
    
    /**
     * 验证图层坐标系
     */
    public static void validateCrs(OguLayer layer) {
        Integer wkid = layer.getWkid();
        
        if (wkid == null) {
            System.out.println("警告: 图层没有坐标系信息");
            return;
        }
        
        try {
            CoordinateReferenceSystem crs = CrsUtil.getCRS(wkid);
            System.out.println("坐标系有效");
            System.out.println("  WKID: " + wkid);
            System.out.println("  名称: " + crs.getName());
            System.out.println("  类型: " + (CrsUtil.isProjectedCRS(crs) ? "投影" : "地理"));
            
            // 检查坐标值范围
            validateCoordinateRange(layer, CrsUtil.isProjectedCRS(crs));
            
        } catch (Exception e) {
            System.out.println("错误: 无效的WKID " + wkid);
        }
    }
    
    /**
     * 验证坐标值范围
     */
    private static void validateCoordinateRange(OguLayer layer, boolean isProjected) {
        for (OguFeature feature : layer.getFeatures()) {
            Geometry geom = GeometryUtil.wkt2Geometry(feature.getWkt());
            Coordinate[] coords = geom.getCoordinates();
            
            for (Coordinate coord : coords) {
                if (isProjected) {
                    // 投影坐标范围检查
                    if (coord.x < 0 || coord.x > 1000000000 ||
                        coord.y < 0 || coord.y > 1000000000) {
                        System.out.printf("警告: FID %s 坐标可能异常: (%.2f, %.2f)%n",
                            feature.getFid(), coord.x, coord.y);
                    }
                } else {
                    // 地理坐标范围检查
                    if (coord.x < -180 || coord.x > 180 ||
                        coord.y < -90 || coord.y > 90) {
                        System.out.printf("警告: FID %s 坐标超出范围: (%.4f, %.4f)%n",
                            feature.getFid(), coord.x, coord.y);
                    }
                }
            }
        }
    }
}

8.6.3 批量坐标转换

public class BatchTransform {
    
    /**
     * 批量转换WKT坐标
     */
    public static List<String> batchTransform(List<String> wktList,
            int sourceWkid, int targetWkid) {
        
        List<String> result = new ArrayList<>();
        for (String wkt : wktList) {
            try {
                String transformed = CrsUtil.transform(wkt, sourceWkid, targetWkid);
                result.add(transformed);
            } catch (Exception e) {
                System.err.println("转换失败: " + wkt);
                result.add(null);
            }
        }
        return result;
    }
    
    /**
     * 并行批量转换（大数据量时使用）
     */
    public static List<String> parallelTransform(List<String> wktList,
            int sourceWkid, int targetWkid) {
        
        return wktList.parallelStream()
            .map(wkt -> {
                try {
                    return CrsUtil.transform(wkt, sourceWkid, targetWkid);
                } catch (Exception e) {
                    return null;
                }
            })
            .collect(Collectors.toList());
    }
}

8.6.4 面积计算（带坐标转换）

public class AreaCalculator {
    
    /**
     * 计算准确面积（平方米）
     * 自动处理坐标系转换
     */
    public static double calculateArea(String wkt, int wkid) {
        CoordinateReferenceSystem crs = CrsUtil.getCRS(wkid);
        
        if (CrsUtil.isProjectedCRS(crs)) {
            // 已经是投影坐标系，直接计算
            Geometry geom = GeometryUtil.wkt2Geometry(wkt);
            return GeometryUtil.area(geom);
        } else {
            // 地理坐标系，需要先转换到投影坐标系
            int dh = CrsUtil.getDh(wkt);
            Integer projWkid = CrsUtil.getProjectedWkid(dh);
            
            String projWkt = CrsUtil.transform(wkt, wkid, projWkid);
            Geometry geom = GeometryUtil.wkt2Geometry(projWkt);
            return GeometryUtil.area(geom);
        }
    }
    
    /**
     * 计算图层所有要素的面积
     */
    public static void calculateLayerAreas(OguLayer layer) {
        Integer wkid = layer.getWkid();
        if (wkid == null) {
            System.out.println("图层没有坐标系信息，无法计算准确面积");
            return;
        }
        
        for (OguFeature feature : layer.getFeatures()) {
            double area = calculateArea(feature.getWkt(), wkid);
            feature.setValue("AREA_M2", area);
            feature.setValue("AREA_MU", area / 666.67);  // 转换为亩
            feature.setValue("AREA_HA", area / 10000);   // 转换为公顷
        }
    }
}

8.6.5 WGS84与CGCS2000转换

public class Wgs84ToCgcs2000 {
    
    // CGCS2000与WGS84非常接近，在大多数应用中可以近似等价
    // 如需精确转换，需要使用七参数或四参数转换
    
    /**
     * WGS84 -> CGCS2000（简化，直接改WKID）
     * 适用于大多数应用场景，误差小于1米
     */
    public static OguLayer wgs84ToCgcs2000Simple(OguLayer layer) {
        if (layer.getWkid() != 4326) {
            throw new IllegalArgumentException("输入必须是WGS84坐标系");
        }
        
        // 简单处理：直接修改WKID
        layer.setWkid(4490);
        return layer;
    }
    
    /**
     * CGCS2000地理坐标 -> CGCS2000投影坐标
     */
    public static OguLayer cgcs2000ToProjected(OguLayer layer) {
        if (layer.getWkid() != 4490) {
            throw new IllegalArgumentException("输入必须是CGCS2000地理坐标系");
        }
        
        // 计算合适的带号
        OguFeature firstFeature = layer.getFeatures().get(0);
        int dh = CrsUtil.getDh(firstFeature.getWkt());
        Integer targetWkid = CrsUtil.getProjectedWkid(dh);
        
        return CrsUtil.reproject(layer, targetWkid);
    }
}

8.7 常见问题

8.7.1 坐标系不匹配

问题：叠加分析时发现数据对不上

解决：

// 确保所有图层使用相同坐标系
int targetWkid = 4490;  // 统一使用CGCS2000

OguLayer layer1 = OguLayerUtil.readLayer(...);
OguLayer layer2 = OguLayerUtil.readLayer(...);

// 统一坐标系
if (layer1.getWkid() != targetWkid) {
    layer1 = CrsUtil.reproject(layer1, targetWkid);
}
if (layer2.getWkid() != targetWkid) {
    layer2 = CrsUtil.reproject(layer2, targetWkid);
}

8.7.2 带号选择

问题：数据跨多个带

解决：

// 计算数据中心点的带号
Geometry union = createUnionGeometry(layer);
Geometry centroid = GeometryUtil.centroid(union);
int centerDh = CrsUtil.getDh(centroid);

// 或者使用6度分带（跨度更大）
// 需要自行计算WKID

8.7.3 面积计算误差

问题：使用地理坐标系计算面积，结果不准确

解决：

// 始终转换到投影坐标系后计算面积
if (!CrsUtil.isProjectedCRS(CrsUtil.getCRS(layer.getWkid()))) {
    layer = smartReproject(layer);
}

// 然后计算面积
double area = GeometryUtil.area(GeometryUtil.wkt2Geometry(feature.getWkt()));

8.7.4 坐标精度丢失

问题：转换后坐标精度下降

解决：

// 使用足够的小数位数
// OGU4J默认保留足够精度

// 如果需要控制输出精度，可以在最后格式化
String wkt = feature.getWkt();
// 例如保留6位小数...

8.8 性能优化

8.8.1 转换缓存

public class TransformCache {
    
    private static final Map<String, MathTransform> transformCache = 
        new ConcurrentHashMap<>();
    
    /**
     * 获取缓存的转换对象
     */
    public static MathTransform getTransform(int sourceWkid, int targetWkid) {
        String key = sourceWkid + "->" + targetWkid;
        return transformCache.computeIfAbsent(key, k -> {
            try {
                CoordinateReferenceSystem sourceCRS = CrsUtil.getCRS(sourceWkid);
                CoordinateReferenceSystem targetCRS = CrsUtil.getCRS(targetWkid);
                return CRS.findMathTransform(sourceCRS, targetCRS, true);
            } catch (Exception e) {
                throw new RuntimeException(e);
            }
        });
    }
}

8.8.2 批量转换优化

public class OptimizedTransform {
    
    /**
     * 批量转换（复用Transform对象）
     */
    public static List<Geometry> batchTransform(List<Geometry> geometries,
            int sourceWkid, int targetWkid) {
        
        if (geometries.isEmpty()) {
            return Collections.emptyList();
        }
        
        // 只创建一次Transform
        MathTransform transform = TransformCache.getTransform(sourceWkid, targetWkid);
        GeometryFactory factory = new GeometryFactory();
        
        List<Geometry> results = new ArrayList<>();
        for (Geometry geom : geometries) {
            try {
                Geometry transformed = JTS.transform(geom, transform);
                results.add(transformed);
            } catch (Exception e) {
                results.add(null);
            }
        }
        
        return results;
    }
}

---

← 上一章：几何处理与空间分析 | 下一章：异常处理体系 →