java地理坐标系及投影间转换

一、坐标系简介

1.地心坐标

        地心坐标顾名思义所建立的地球椭球中心为地心，常用地心坐标系有WGS84与CGCS2000，主要参数有长半轴与扁率，地心坐标系椭球之间参数相近。

2.参心坐标

        参心坐标所建立的地球椭球是以某一片区域更为拟合，可以理解为一种局部最优的坐标系，在国内常用参心坐标系有西安80与北京54，各坐标系椭球间主要参数长半轴与扁率参数相差较大，且姿态也有所不同。

二、坐标系间转换

1.七参数转换

        七参数转换遵循布尔莎模型，根据至少三个已知在基础与目标坐标系中坐标的点可以计算得到，布尔莎模型易于查到不多赘述。

        上图中简单表示了不同坐标系间的姿态关系，以及七参数的实际意义。

        可以看到，坐标变换是针对平面坐标的，若获得的是经纬度，可根据后续提到的方法进行经纬度到平面坐标的投影转换。

public double[] sevenParamTrans(Datum7Paras datum7Paras,double X,double Y,double Z) {
        double[] point = new double[3];
        double dx = datum7Paras.getDx();
        double dy = datum7Paras.getDy();
        double dz = datum7Paras.getDz();
        double rx = datum7Paras.getRx() * 0.0000048481373323; // 1秒=0.0000048481373323 弧度
        double ry = datum7Paras.getRy() * 0.0000048481373323;
        double rz = datum7Paras.getRz() * 0.0000048481373323;
        double m = datum7Paras.getPpm();
        double k = m / 1000000;

        point[0] = (1 + k) * (X + this.rz * Y - this.ry * Z) + this.dx;
        point[1] = (1 + k) * (-this.rz * X + Y + this.rx * Z) + this.dy;
        point[2] = (1 + k) * (this.ry * X - this.rx * Y + Z) + this.dz;
}

class Datum7Paras{
    private double dx;
    private double dy;
    private double dz;
    private double rx;
    private double ry;
    private double rz;
    private double ppm;

    public Datum7Paras(double dx, double dy, double dz, double rx, double ry, double rz, double ppm) {
        this.dx = dx;
        this.dy = dy;
        this.dz = dz;
        this.rx = rx;
        this.ry = ry;
        this.rz = rz;
        this.ppm = ppm;
    }

    public double getDx() {
        return dx;
    }

    public void setDx(double dx) {
        this.dx = dx;
    }

    public double getDy() {
        return dy;
    }

    public void setDy(double dy) {
        this.dy = dy;
    }

    public double getDz() {
        return dz;
    }

    public void setDz(double dz) {
        this.dz = dz;
    }

    public double getRx() {
        return rx;
    }

    public void setRx(double rx) {
        this.rx = rx;
    }

    public double getRy() {
        return ry;
    }

    public void setRy(double ry) {
        this.ry = ry;
    }

    public double getRz() {
        return rz;
    }

    public void setRz(double rz) {
        this.rz = rz;
    }

    public double getPpm() {
        return ppm;
    }

    public void setPpm(double ppm) {
        this.ppm = ppm;
    }
}

2.四参数转换

        在已知两点一下的基础和目标坐标信息使，利用四参数转换也可以实现，但误差会更大，与七参数类似，但四参数转换不考虑高程，只针对平面的北坐标与东坐标。

public double[] transform4Para(Trans4Paras transPara,double X,double Y) {
        double X1 = transPara.getDx();
        double Y1 = transPara.getDy();

        double cosAngle = Math.cos(transPara.getA());
        double sinAngle = Math.sin(transPara.getA());

        double[] point = new double[2];
        point[0] = X1 + transPara.getK() * (cosAngle * X - sinAngle * Y);
        point[1] = Y1 + transPara.getK() * (sinAngle * X + cosAngle * Y);
    }

class Trans4Paras{
    private double _dx;
    private double _dy;
    private double _a;
    private double _k;
  
    public double getDx() {
        return _dx;
    }

    public void setDx(double dx) {
        _dx = dx;
    }

    public double getDy() {
        return _dy;
    }

    public void setDy(double dy) {
        _dy = dy;
    }

    public double getA() {
        return _a;
    }

    public void setA(double a) {
        _a = a;
    }

    public double getK() {
        return _k;
    }

    public void setK(double k) {
        _k = k;
    }

    public Trans4Paras(double dx, double dy, double a, double k) {
        _dx = dx;
        _dy = dy;
        _a = a;
        _k = k;
    }

}

三、投影转换

1.EPSG编码与依赖

        EPSG编码是一种用于定义不同坐标系的编码，在java中利用proj4j库配合EPSG编码即可实现不同投影转换。其实上面的不同坐标系下坐标变换也可以用此方法实现，但代码都较简单不多赘述。

        在Maven中添加下方依赖即可完成proj4j库的添加，注意版本号自行修改。

<dependency>
            <groupId>org.osgeo</groupId>
            <artifactId>proj4j</artifactId>
            <version>0.1.0</version>
</dependency>

2.投影方式

        UTM投影：

                特点：保长度不保方向，适用南纬80到北纬84度，通常应用于全球坐标，因此以六度分带，且分带从西经180度起始，对应国内带号加30。

                EPSG：32601-32660，按带号修改后两位。

        Web墨卡托投影：

                特点：保方向不保长度，适用于服务集成交互。

                EPSG：3857，只有一个是因为该投影方式认为地球为规则球体。

        高斯克吕格投影：

                特点：与UTM类似但包含所有纬度，多用于国内的坐标。

                EPSG：CGCS2000：4513-4554；北京54：2410-2442；西安80：2349-2390

        各椭球的EPSG：WGS84：4326；CGCS2000：4490；西安80：4610；北京54：4214

3.示例

        简单示例，列出UTM投影、WGS84_N坐标系时，满足转换的两坐标正反投影的结果。

import org.osgeo.proj4j.*;

public class PointTest {
    public static void main(String[] args) {
        ToCoordinate(579573.572, 4269411.894,"EPSG:32651","EPSG:4326",true);
        ToCoordinate(123.91343418163903, 38.569600138042816,"32651","4326",false);
    }

    public static void ToCoordinate(double argX,double argY,String sys1,String sys2,boolean xyTolb){
        // 创建CRS工厂
        CRSFactory crsFactory = new CRSFactory();

        // 创建UTM坐标参考系统（CRS），这里以UTM 33N为例
        CoordinateReferenceSystem utmCrs = crsFactory.createFromName("EPSG:" + sys1);

        // 创建WGS84坐标参考系统
        CoordinateReferenceSystem wgs84Crs = crsFactory.createFromName("EPSG:" + sys2);

        // 创建坐标转换工厂
        CoordinateTransformFactory transformFactory = new CoordinateTransformFactory();
        if (xyTolb){
            // 创建从UTM到WGS84的坐标转换器
            CoordinateTransform transformToWgs84 = transformFactory.createTransform(utmCrs, wgs84Crs);
            // 创建UTM坐标点
            ProjCoordinate utmPoint = new ProjCoordinate(argX, argY); // X, Y in meters

            // 转换到WGS84坐标点
            ProjCoordinate wgs84Point = new ProjCoordinate();
            transformToWgs84.transform(utmPoint, wgs84Point);

            // 打印WGS84坐标
            System.out.println("WGS84: " + wgs84Point.x + ", " + wgs84Point.y);
        }else {
            // 如果需要，可以创建从WGS84到UTM的坐标转换器
            CoordinateTransform transformToUtm = transformFactory.createTransform(wgs84Crs,utmCrs);
            // 创建UTM坐标点
            ProjCoordinate wgs84Point = new ProjCoordinate(argX, argY);
            // 转换回UTM坐标点
            ProjCoordinate utmPointBack = new ProjCoordinate();
            transformToUtm.transform(wgs84Point, utmPointBack);

            // 打印转换回的UTM坐标
            System.out.println("UTM: " + utmPointBack.x + ", " + utmPointBack.y);
        }

    }
}