高斯坐标和经纬度之间的转化

最近项目需要做工程项目在地图的标注，项目的数据是高斯投影坐标，而地图的标注需要用到的是经纬度，因此需要做处理。在网上查找了各种帖子，掉了不少坑，最终完善出了一个工具类，亲测可用，有问题可以提出探讨。

首先是常用坐标系的参数：

（1）北京54坐标系

长半轴（米）：6378245.0； 短半轴（米）：6356863 ；扁率之倒数：298.3；中央子午线比例系数：1；北向偏离（米）：0；东向偏离（米）：500000

（2）国家80坐标系

长半轴（米）：6378140.0； 短半轴（米）：6356755；扁率之倒数：298.25722101；中央子午线比例系数：1；北向偏离（米）：0；东向偏离（米）：500000

（3）WGS84坐标系

长半轴（米）：6378137.0； 短半轴（米）：6356755；扁率之倒数：298.257223563；中央子午线比例系数：1；北向偏离（米）：0；东向偏离（米）：500000

（4）92坐标系（来源于北京54坐标系）

长半轴（米）：6378245.0； 短半轴（米）：6356863 ；扁率之倒数：298.3；中央子午线比例系数：1；北向偏离（米）：2700000；东向偏离（米）：400000

@Slf4j
public class GaussBLUtil {

    public static void main(String[] args) {
        // 以北京54坐标系为例，其中中央子午线（第3个参数）和偏离参数（第4、5个参数）可根据实际情况调整。第6个参数是长半轴，第7个参数是系数除于扁率
        GaussBLUtil.gaussToBL(2720076.481,  454802.918，118.5，0， 500000，6378245.0，1/ 298.3);
    }

    /**
     * @description 高斯投影坐标（平面坐标）转化成经纬度（大地坐标） tip:不带号，直接指定中央经度
     * @param X 平面坐标x
     * @param Y 平面坐标y
     * @param centerL   中央子午线（度）
     * @param X0    北向偏离【米】（X向加常数）
     * @param Y0    东向偏离【米】（Y向加常数）
     * @param factor    卫星椭球坐标投影到平面地图坐标系的投影因子
     * @param param 坐标系参数
     * @date 2023/4/12 15:25
     * @author lyb
     * @return double[]
     */
    public static double[] gaussToBL(double X, double Y, double centerL, double X0, double Y0, double factor, double param){
        int ProjNo;
        // 带宽
        int ZoneWide;
        double[] output = new double[2];
        // latitude0,
        double longitude1, latitude1, longitude0, yval, xval;
        double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
        //3.1415926535898/180.0;
        iPI = 0.0174532925199433;
        a = factor;// 卫星椭球坐标投影到平面地图坐标系的投影因子
        f = param;// 54年北京坐标系参数
        ZoneWide = 6;

        // 中央经线
        longitude0 = centerL * iPI;
        // 带内大地坐标
        log.info("带内大地坐标(偏移前) x={},y={}", X, Y);
        yval = Y-Y0;
        xval = X-X0;
        e2 = 2 * f - f * f;
        e1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2));
        ee = e2 / (1 - e2);
        M = xval;
        u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
        fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32)
                * Math.sin(4 * u) + (151 * e1 * e1 * e1 / 96) * Math.sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.sin(8 * u);
        C = ee * Math.cos(fai) * Math.cos(fai);
        T = Math.tan(fai) * Math.tan(fai);
        NN = a / Math.sqrt(1.0 - e2 * Math.sin(fai) * Math.sin(fai));
        R = a * (1 - e2) / Math.sqrt((1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin
                (fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)));
        D = yval / NN;
        // 计算经度(Longitude) 纬度(Latitude)
        longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
                * D * D * D * D / 120) / Math.cos(fai);
        latitude1 = fai - (NN * Math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
                + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
        // 转换为度 DD
        output[0] = longitude1 / iPI;
        output[1] = latitude1 / iPI;
        log.info("大地坐标L={}度, B={}度", output[0], output[1]);
        return output;
    }

    /**
     * @description 由经纬度（大地坐标）转化成高斯投影坐标（平面坐标） tip:不带号，直接指定中央经度
     * @param longitude 经度
     * @param latitude  纬度
     * @param centerL   中央子午线（度）
     * @param X0    北向偏离【米】（X向加常数）
     * @param Y0    东向偏离【米】（Y向加常数）
     * @param factor    卫星椭球坐标投影到平面地图坐标系的投影因子
     * @param param 坐标系参数
     * @date 2023/4/12 11:02
     * @author lyb
     * @return void
     */
    public static double[] bLToGauss(double longitude, double latitude, double centerL, double X0, double Y0, double factor, double param) {
        int ProjNo = 0;
        int ZoneWide; // //带宽
        double[] output = new double[2];
        double longitude1, latitude1, longitude0, xval, yval;
        double a, f, e2, ee, NN, T, C, A, M, iPI;
        iPI = 0.0174532925199433; // //3.1415926535898/180.0;
        ZoneWide = 6; // 6度带宽
        a = factor;
        f = param;
        longitude0 = centerL;
        longitude0 = longitude0 * iPI; // 中央经线
        longitude1 = longitude * iPI; // 经度转换为弧度
        latitude1 = latitude * iPI; // 纬度转换为弧度
        e2 = 2 * f - f * f;
        ee = e2 * (1.0 - e2);
        NN = a/Math.sqrt(1.0 - e2 * Math.sin(latitude1) * Math.sin(latitude1));
        T = Math.tan(latitude1) * Math.tan(latitude1);
        C = ee * Math.cos(latitude1) * Math.cos(latitude1);
        A = (longitude1 - longitude0) * Math.cos(latitude1);
        M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1
                - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
                / 1024) * Math.sin(2 * latitude1)
                + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
                * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
                * Math.sin(6 * latitude1));

        // 因为是以赤道为Y轴的，与我们南北为Y轴是相反的，所以xy与高斯投影的标准xy正好相反;
        yval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee)
                * A * A * A * A * A / 120);
        xval = M + NN
                * Math.tan(latitude1)
                * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61
                - 58 * T + T * T + 600 * C - 330 * ee)
                * A * A * A * A * A * A / 720);
        yval = yval+Y0; xval = xval+X0;
        output[0] = xval;
        output[1] = yval;
        log.info("平面坐标 x={},y={}", output[0], output[1]);
        return output;
    }

}