CGCS2000地理坐标系转投影坐标系 坑

需求：计算jts的Geometry对象面积

遇到的坑：

1. 地理坐标系（EPSG:4490）得到的面积单位是度，不是常规意义上的面积，所以不能直接用地理坐标系的对象计算面积；
2. geotools库里用于坐标系转换的包没加，找不到EPSG系列坐标系标准；
3. 4490对应的地区投影坐标系，傻傻分不清楚。

解决办法：

有问题就一个个解决。

1. 计算面积前先做投影转换，把地理坐标系转换成投影坐标系再计算面积。

        创建工具类：

package com.zjplan.internal.workbench.common.util;

import lombok.extern.slf4j.Slf4j;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

/**
 * 坐标系转换工具类
 */
@Slf4j
public class ReprojectUtil {
    private CoordinateReferenceSystem sourceCRS;
    private CoordinateReferenceSystem targetCRS;

    public ReprojectUtil(String srcCrsCode, String tarCrsCode, boolean longFirst) {
        try{
            sourceCRS = CRS.decode(srcCrsCode, longFirst);
            targetCRS = CRS.decode(tarCrsCode, longFirst);
        }catch (FactoryException e){
            throw new RuntimeException(e);
        }
    }

    public Geometry reproject(Geometry geometry){
        Geometry reprojectGeometry = null;
        try {

            MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
            reprojectGeometry = JTS.transform(geometry, transform);
        }catch (Exception e){
            throw new RuntimeException(e);
        }
        return reprojectGeometry;
    }

    public Envelope reproject(Envelope envelope){
        Envelope reprojectEnvelope = null;
        try {
            MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
            reprojectEnvelope = JTS.transform(envelope, transform);
        }catch (Exception e){
            throw new RuntimeException(e);
        }
        return reprojectEnvelope;
    }
}

调用：

String SRC_CRS_CODE = "EPSG:4490";
String TAR_CRS_CODE = "EPSG:4549";
ReprojectUtil util = new ReprojectUtil(SRC_CRS_CODE, TAR_CRS_CODE, true);
Geometry geomTar = util.reproject(geom);
Double area = geomTar.getArea();

 这就是完整的解决方法，如果这么操作没问题的，可以不用往下看了。

2. 报错No code "EPSG:4326" from authority "EPSG" found for object of type "EngineeringCRS"

原因是没有添加坐标系转换对应的包。版本就和本系统添加的geotools的版本对齐。我这边是27.0

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-hsql</artifactId>
            <version>${geotools.version}</version>
        </dependency>

3. 报错 Bursa wolf parameters required./ No transformation available from system "CartesianCS[Geocentric]" to "CartesianCS[Earth centred, earth fixed, righthanded 3D coordinate system, consisting of 3 (...) positive Z-axis parallel to mean earth rotation axis and pointing towards North Pole...

原因是转换的目标坐标系和源坐标系基准椭球面不同或者维度不同（二维转三维），导致无法转换。我出现这些问题的原因是第一次转换坐标系用了3857，以为墨卡托治百病；第二次找了CGCS2000的单位为米的坐标系，结果发现是地心坐标系，一定需要有z轴坐标，无法转换。

确定目标：一定要找投影坐标系！那如何找到最符合自己所在地地理位置的坐标系呢？这样得出的面积更符合精确值。

这就需要用到GIS软件了。

最终的解决方法是：在ArcGIS-工具箱-投影里，点击小地球，选择“投影坐标系”-CGCS2000，里面根据自己所在地的经度确定投影坐标系即可。

我这边的经度在120°E，所以选择的是投影坐标系    CGCS2000 3 Degree GK CM 120E，点开详细信息，可以看到WKID是4549。所以在代码中，我的目标坐标系就是4549。