北斗三维网格位置码介绍与实现详解（非极地）

北斗卫星导航系统是中国自主研发的全球卫星导航系统，具备高精度、高可靠性的定位、导航和授时功能。在北斗系统中，三维网格位置码（Beidou Grid Code）是一种用于描述三维空间位置的编码方式，特别适用于需要高精度空间定位和管理的场景。

规范见百度百科内容范围章节图片：北斗网格位置码_百度百科

二维北斗网格位置码实现讲解见：北斗网格位置码详解：经纬度到二维网格码的转换（非极地）_北斗网格位置编码算法-CSDN博客

1. 北斗网格位置码介绍

北斗三维网格位置码是一种基于地球表面和空间划分的编码机制。它将地球表面划分为若干个二维网格单元，并结合高度维度，形成三维网格结构。每个三维网格单元具有唯一的编码标识，便于快速定位、检索和管理地理信息。

2. 北斗三维网格位置码编码规则

北斗三维网格位置码由二维编码 + 高度维（三维）编码交叉组成，共32位码元组成，北斗三维网格位置码（第三维度编码部分）由12位码元组成，其结构与码元取值如图所示。

​

第三维度编码如下:

1. 第一位码元，地上、地下标识用 0、1 进行表示，其中：地表以下，用1标识、地表以上，标识为0，对应初始网格；
2. 第二位、第三位码元，编码标识采用 00 ~ 63，对应第一级网格；
3. 第四位码元，编码标识采用 0 ~ 7，对应第二级网格；
4. 第五位码元，编码标识采用 0 ~ 1，对应第三级网格；
5. 第六位码元，编码标识采用 0 ~ 9、A ~ E，对应第四级网格；
6. 第七位码元，编码标识采用 0 ~ 9、A ~ E，对应第五级网格；
7. 第八位码元，编码标识采用 0 ~ 1，对应第六级网格；
8. 第九位码元，编码标识采用 0 ~ 7，对应第七级网格；
9. 第十位码元,编码标识采用 0 ~ 7，对应第八级网格；
10. 第十一位码元，编码标识采用 0 ~ 7，对应第九级网格；
11. 第十二位码元,编码标识采用 0 ~ 7，对应第十级网格。

 3. 规范解析

从上面规范截图可知，整个空域的高度划分范围从大地高 -6302.106722602182 km ~ 528680.1711252437 km 的范围。1 ~ 10层级高度维取值均从第 10 层级的高度维编码截取，故求出第 10 层级的高度维值即可得出其他高度维编码。第 10 层级的高度维编码根据 C.12 和 C.14 的公式得出。其中，C.14 的公式 n(a:b) 中，n 为 C.14 公式的 32 位二进制值；a，b的跨度是根据北斗三维网格位置码编码规则中对应等级 bit 位的值；n(a:b) 则为取 C.14 公式的 32 位二进制值的第 a 至 b位的值。(*)2，(*)8，(*)16，(*)64 则为 n(a:b) 转十进制后再转的进制。

4. 代码实现

1. 添加依赖：

<dependency>
    <groupId>ch.obermuhlner</groupId>
    <artifactId>big-math</artifactId>
    <version>2.0.0</version>
</dependency>

2. 创建异常类：

public class BeidouGridException extends RuntimeException{
    private static final long serialVersionUID = 1L;
 
    public BeidouGridException() {
    }
 
    public BeidouGridException(String message) {
        super(message);
    }
 
    public BeidouGridException(String message, Throwable cause) {
        super(message, cause);
    }
 
    public BeidouGridException(Throwable cause) {
        super(cause);
    }
 
    public BeidouGridException(String message, Throwable cause, boolean enableSuppression, boolean writableStackTrace) {
        super(message, cause, enableSuppression, writableStackTrace);
    }
}

3. 创建网格信息实体类：

import com.zjp.exceptions.BeidouGridException;
import lombok.AllArgsConstructor;
import lombok.Builder;
import lombok.Data;
import lombok.NoArgsConstructor;
import lombok.experimental.Accessors;
 
import java.math.BigDecimal;
import java.math.MathContext;
 
@Data
@Builder
@Accessors(chain = true)
@NoArgsConstructor
@AllArgsConstructor
public class GridInfo {
    // 常量定义（1度=3600秒）
    private static final BigDecimal SECONDS_PER_DEGREE = new BigDecimal("3600");
    private static final BigDecimal SECONDS_PER_MINUTE = new BigDecimal("60");
    // 纬度步长常量定义（单位：秒）
    private static final BigDecimal[] LON_STEPS = {
            // level 1: 6°=6*3600=21600秒
            new BigDecimal("6").multiply(SECONDS_PER_DEGREE),
            // level 2: 30'=30*60=1800秒
            new BigDecimal("30").multiply(SECONDS_PER_MINUTE),
            // level 3: 15'=15*60=900秒
            new BigDecimal("15").multiply(SECONDS_PER_MINUTE),
            // level 4: 1'=1*60=60秒
            SECONDS_PER_MINUTE,
            // level 5: 4秒
            new BigDecimal("4"),
            // level 6: 2秒
            new BigDecimal("2"),
            // level 7: 0.25秒
            new BigDecimal("0.25"),
            // level 8: 1/32秒=0.03125秒
            new BigDecimal("1").divide(new BigDecimal("32"), MathContext.DECIMAL128),
            // level 9: 1/256秒=0.00390625秒
            new BigDecimal("1").divide(new BigDecimal("256"), MathContext.DECIMAL128),
            // level 10: 1/2048秒=0.00048828125秒
            new BigDecimal("1").divide(new BigDecimal("2048"), MathContext.DECIMAL128)
    };
    // 经度步长常量定义（单位：秒）
    private static final BigDecimal[] LAT_STEPS = {
            // level 1: 4°=4*3600=14400秒
            new BigDecimal("4").multiply(SECONDS_PER_DEGREE),
            // level 2: 30'=30*60=1800秒
            new BigDecimal("30").multiply(SECONDS_PER_MINUTE),
            // level 3: 10'=10*60=600秒
            new BigDecimal("10").multiply(SECONDS_PER_MINUTE),
            // level 4: 1'=60秒
            SECONDS_PER_MINUTE,
            // level 5: 4秒
            new BigDecimal("4"),
            // level 6: 2秒
            new BigDecimal("2"),
            // level 7: 0.25秒
            new BigDecimal("0.25"),
            // level 8: 1/32秒=0.03125秒
            new BigDecimal("1").divide(new BigDecimal("32"), MathContext.DECIMAL128),
            // level 9: 1/256秒=0.00390625秒
            new BigDecimal("1").divide(new BigDecimal("256"), MathContext.DECIMAL128),
            // level 10: 1/2048秒=0.00048828125秒
            new BigDecimal("1").divide(new BigDecimal("2048"), MathContext.DECIMAL128)
    };

    /**
     * 初始经度
     */
    private BigDecimal baseLon;
    /**
     * 初始纬度
     */
    private BigDecimal baseLat;
    /**
     * 高程方向的整数编码二进制
     */
    private String heightCodeBinary;
    /**
     * 左下角经度
     */
    private BigDecimal currentLon;
    /**
     * 左下角纬度
     */
    private BigDecimal currentLat;
    /**
     * 当前列号
     */
    private int column;
    /**
     * 当前行号
     */
    private int row;
    /**
     * 网格粒度经度
     */
    private BigDecimal gridStepLon;
    /**
     * 网格粒度纬度
     */
    private BigDecimal gridStepLat;
    /**
     * 网格粒度
     */
    private int gridStep = 1;
    /**
     * 网格编码
     */
    private String gridCode;

    public BigDecimal getGridStepLon() {
        if (gridStep < 1 || gridStep > LON_STEPS.length) {
            throw new BeidouGridException("编码级别超出范围: " + gridStep);
        }
        return LON_STEPS[gridStep - 1];
    }

    public BigDecimal getGridStepLat() {
        if (gridStep < 1 || gridStep > LON_STEPS.length) {
            throw new BeidouGridException("编码级别超出范围: " + gridStep);
        }
        return LAT_STEPS[gridStep - 1];
    }
}

4. 北斗网格位置码工具类：

import ch.obermuhlner.math.big.BigDecimalMath;
import com.zjp.exceptions.BeidouGridException;
import com.zjp.pojo.GridInfo;
import lombok.experimental.UtilityClass;
import lombok.extern.slf4j.Slf4j;
 
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.function.Function;

@Slf4j
@UtilityClass
public class BeidouGridUtil {
    // 3D编码策略
    private static final GridEncoder<GridInfo> GRID_ENCODER_3D = new GridEncoder<>();

    static {
        // 注册3D编码器
        register3DEncoders();
    }

    /**
     * 注册所有3D编码器
     */
    private void register3DEncoders() {
        GRID_ENCODER_3D.register(1, BeidouGridUtil::calculate3DFirstLevelCode);
        GRID_ENCODER_3D.register(2, BeidouGridUtil::calculate3DSecondLevelCode);
        GRID_ENCODER_3D.register(3, BeidouGridUtil::calculate3DThirdLevelCode);
        GRID_ENCODER_3D.register(4, BeidouGridUtil::calculate3DForthLevelCode);
        GRID_ENCODER_3D.register(5, BeidouGridUtil::calculate3DFifthLevelCode);
        GRID_ENCODER_3D.register(6, BeidouGridUtil::calculate3DSixthLevelCode);
        GRID_ENCODER_3D.register(7, BeidouGridUtil::calculate3DSeventhLevelCode);
        GRID_ENCODER_3D.register(8, BeidouGridUtil::calculate3DEighthLevelCode);
        GRID_ENCODER_3D.register(9, BeidouGridUtil::calculate3DNinthLevelCode);
        GRID_ENCODER_3D.register(10, BeidouGridUtil::calculate3DTenthLevelCode);
    }

    /**
     * 生成三维北斗网格码
     *
     * @param longitude 经度(单位：秒)
     * @param latitude  纬度(单位：秒)
     * @param height    高程(单位：米)
     * @param level     编码级别(1-10)
     * @return 北斗网格编码
     */
    public String generate3DGridCode(BigDecimal longitude, BigDecimal latitude, BigDecimal height, int level) {
        validateCoordinates(longitude.doubleValue(), latitude.doubleValue());
        validateHeight(height);
        validateLevel(level);

        GridInfo gridInfo = createInitialGridInfo(longitude, latitude, height);
        for (int i = 1; i <= level; i++) {
            Function<GridInfo, GridInfo> encoder = GRID_ENCODER_3D.get(i);
            if (encoder == null) {
                throw new BeidouGridException("不支持的编码级别: " + i);
            }
            gridInfo = encoder.apply(gridInfo);
        }
        return gridInfo.getGridCode();
    }



    /**
     * 校验经纬度是否在有效范围内（经度：-180°~180°，纬度：-88°~88°），不包含极地
     *
     * @param longitude 经度(单位：秒)
     * @param latitude  纬度(单位：秒)
     */
    private void validateCoordinates(double longitude, double latitude) {
        if (longitude < -180 * 3600 || longitude > 180 * 3600 ||
                latitude < -88 * 3600 || latitude > 88 * 3600) {
            throw new BeidouGridException("经纬度超出范围");
        }
    }

    /**
     * 校验高度是否在有效范围内（-6302106.722602182m ~ 528680171.1252437m）
     *
     * @param height 高度(单位：米)
     */
    private void validateHeight(BigDecimal height) {
        if (height.compareTo(new BigDecimal("-6302106.722602182")) < 0 || height.compareTo(new BigDecimal("528680171.1252437")) > 0) {
            throw new BeidouGridException("高度超出范围");
        }
    }

    /**
     * 校验编码级别是否在有效范围内（1-10）
     *
     * @param level 编码级别
     */
    private void validateLevel(int level) {
        if (level < 1 || level > 10) {
            throw new BeidouGridException("编码级别超出范围: " + level);
        }
    }


    /**
     * 创建初始网格信息对象（用于三维编码）
     *
     * @param longitude 初始经度(单位：秒)
     * @param latitude  初始纬度(单位：秒)
     * @param height    初始高度(单位：米)
     * @return 初始化后的 GridInfo 对象
     */
    private GridInfo createInitialGridInfo(BigDecimal longitude, BigDecimal latitude, BigDecimal height) {
        // 获取地球半径（单位：米）
        BigDecimal earthRadius = BigDecimal.valueOf(6378137);
        // 获取基础纬度（单位：弧度）
        BigDecimal theta0 = BigDecimal.valueOf(Math.PI).divide(new BigDecimal("180"), MathContext.DECIMAL128);
        // 计算网格步长（单位：米）
        BigDecimal heightGridStep = BigDecimal.ONE.divide(new BigDecimal("2048")
                .multiply(new BigDecimal(3600)), MathContext.DECIMAL128);
        // 计算 (1 + theta_0)
        BigDecimal theta = BigDecimal.ONE.add(theta0);
        // 计算 n = theta_0 / (theta_0 * n) * ( log(rn/EARTH_RADIUS) /log(1 + theta_0))
        BigDecimal n = BigDecimal.ONE
                .divide(heightGridStep, MathContext.DECIMAL128)
                .multiply(BigDecimalMath.log(earthRadius.add(height)
                                .divide(earthRadius, MathContext.DECIMAL128), MathContext.DECIMAL128)
                        .divide(BigDecimalMath.log(theta, MathContext.DECIMAL128), MathContext.DECIMAL128));

        // 将 n 转换为32位二进制
        String heightCodeBinary = to32BitBinary(n.longValue());

        return GridInfo.builder()
                .baseLon(longitude)
                .baseLat(latitude)
                .heightCodeBinary(heightCodeBinary)
                .currentLon(longitude)
                .currentLat(latitude)
                .gridStep(1)
                .build();
    }

    /**
     * 将数字转换为32位二进制字符串
     *
     * @param number 数字
     * @return 32位二进制字符串
     */
    private String to32BitBinary(long number) {
        // 转换为二进制字符串
        String binary = Long.toBinaryString(number);

        // 取最后32位（自动兼容负数补码）
        if (binary.length() > 32) {
            binary = binary.substring(binary.length() - 32);
        }

        // 补前导0，确保长度为32
        return String.format("%32s", binary).replace(' ', '0');
    }


    /**
     * 计算第一级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度和高度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和网格编码
     */
    private GridInfo calculate3DFirstLevelCode(GridInfo gridInfo) {
        BigDecimal longitude = gridInfo.getBaseLon();
        BigDecimal latitude = gridInfo.getBaseLat();
        String heightCodeBinary = gridInfo.getHeightCodeBinary();

        // 北斗网格编码使用的32进制字符集
        char[] rowChars = "ABCDEFGHIJKLMNOPQRSTUV".toCharArray();

        // 判断半球并计算绝对纬度值
        char hemisphere = latitude.compareTo(BigDecimal.ZERO) >= 0 ? 'N' : 'S';

        // 计算经度方向编号（列号）
        int column = longitude
                .add(BigDecimal.valueOf(180 * 3600)) // 加上东经180°的偏移量
                .divide(gridInfo.getGridStepLon(), 0, RoundingMode.DOWN)
                .add(BigDecimal.ONE)
                .intValue();

        // 计算纬度方向编号（行号）
        int row = latitude.abs()
                .divide(gridInfo.getGridStepLat(), 0, RoundingMode.DOWN)
                .add(BigDecimal.ONE)
                .intValue();

        // 将纬度编号映射到对应的字母（1->A, 2->B,...）
        char rowChar;
        if (latitude.compareTo(BigDecimal.ZERO) < 0) {
            // 南半球使用倒序映射
            rowChar = rowChars[rowChars.length - row];
        } else {
            // 北半球使用正常映射
            rowChar = rowChars[row - 1];
        }

        // 计算高程编号
        String groundUndergroundFlag = heightCodeBinary.substring(0, 1);
        String layer = Long.toUnsignedString(Long.parseLong(heightCodeBinary.substring(1, 7), 2), 64);

        // 返回更新后的网格信息
        return gridInfo.setColumn(column)
                .setRow(row)
                .setGridCode(String.format("%c%s%02d%c%02d",
                        hemisphere,
                        groundUndergroundFlag,
                        column,
                        rowChar,
                        Integer.valueOf(layer)));
    }

    /**
     * 计算第二级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DSecondLevelCode(GridInfo gridInfo) {
        BigDecimal longitude = gridInfo.getBaseLon();
        BigDecimal latitude = gridInfo.getBaseLat();
        String heightCodeBinary = gridInfo.getHeightCodeBinary();

        // 计算基准经纬度 λ₁, φ₁
        BigDecimal baseLongitude = (longitude.compareTo(BigDecimal.ZERO) < 0
                ? new BigDecimal(gridInfo.getColumn())
                : new BigDecimal(gridInfo.getColumn() - 1))
                .multiply(gridInfo.getGridStepLon())
                .subtract(new BigDecimal("648000")); // 180°*3600秒

        BigDecimal baseLatitude = new BigDecimal(gridInfo.getRow() - 1)
                .multiply(gridInfo.getGridStepLat())
                .multiply(BigDecimal.valueOf(latitude.compareTo(BigDecimal.ZERO)));

        // 更新网格粒度级别为2
        gridInfo.setGridStep(2);

        // 计算子网格列号和行号
        int column = longitude
                .subtract(baseLongitude)
                .divide(gridInfo.getGridStepLon(), 0, RoundingMode.DOWN)
                .abs()
                .add(BigDecimal.ONE)
                .intValue();

        int row = latitude
                .subtract(baseLatitude)
                .divide(gridInfo.getGridStepLat(), 0, RoundingMode.DOWN)
                .abs()
                .add(BigDecimal.ONE)
                .intValue();

        // 计算高度方向编号
        String layer = Long.toOctalString(Long.parseLong(heightCodeBinary.substring(7, 10), 2));

        // 定义编码字符集
        char[] columnChars = "0123456789AB".toCharArray();
        char[] rowChars = "01234567".toCharArray();

        // 返回更新后的网格信息
        return gridInfo.setCurrentLon(baseLongitude)
                .setCurrentLat(baseLatitude)
                .setColumn(column)
                .setRow(row)
                .setGridCode(String.format("%s%c%c%s",
                        gridInfo.getGridCode(),
                        rowChars[column - 1],
                        columnChars[row - 1],
                        layer));
    }

    /**
     * 计算第三级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DThirdLevelCode(GridInfo gridInfo) {
        int[][] zOrderMatrix = {
                {0, 2, 4},
                {1, 3, 5}
        };
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = heightCodeBinary.substring(10, 11);
        return calculate3DLevelCode(gridInfo, 3, zOrderMatrix, layer);
    }

    /**
     * 计算第四级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DForthLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toHexString(Long.parseLong(heightCodeBinary.substring(11, 15), 2));
        return calculate3DLevelCode(gridInfo, 4, false, false, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789".toCharArray(),
                layer);
    }

    /**
     * 计算第五级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DFifthLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toHexString(Long.parseLong(heightCodeBinary.substring(15, 19), 2));
        return calculate3DLevelCode(gridInfo, 5, true, true, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789ABCDE".toCharArray(),
                layer);
    }

    /**
     * 计算第六级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DSixthLevelCode(GridInfo gridInfo) {
        int[][] zOrderMatrix = {
                {0, 2},
                {1, 3}
        };
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = heightCodeBinary.substring(19, 20);
        return calculate3DLevelCode(gridInfo, 6, zOrderMatrix, layer);
    }

    /**
     * 计算第七级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DSeventhLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toOctalString(Long.parseLong(heightCodeBinary.substring(20, 23), 2));
        return calculate3DLevelCode(gridInfo, 7, false, false, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789ABCDE".toCharArray(),
                layer);
    }

    /**
     * 计算第八级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DEighthLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toOctalString(Long.parseLong(heightCodeBinary.substring(23, 26), 2));
        return calculate3DLevelCode(gridInfo, 8, true, true, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789ABCDE".toCharArray(),
                layer);
    }

    /**
     * 计算第九级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DNinthLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toOctalString(Long.parseLong(heightCodeBinary.substring(26, 29), 2));
        return calculate3DLevelCode(gridInfo, 9, true, true, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789ABCDE".toCharArray(),
                layer);
    }

    /**
     * 计算第十级北斗三维网格编码
     *
     * @param gridInfo 网格信息对象，包含基础经纬度、高度及当前网格粒度信息
     * @return 包含计算结果的网格信息对象，包括列号、行号、高程号和更新后的网格编码
     */
    private GridInfo calculate3DTenthLevelCode(GridInfo gridInfo) {
        String heightCodeBinary = gridInfo.getHeightCodeBinary();
        String layer = Long.toOctalString(Long.parseLong(heightCodeBinary.substring(29, 32), 2));
        return calculate3DLevelCode(gridInfo, 10, true, true, true,
                "0123456789ABCDE".toCharArray(),
                "0123456789ABCDE".toCharArray(),
                layer);
    }

    /**
     * 计算三维网格编码信息
     * 该方法用于根据当前网格信息和指定的粒度级别，计算并更新网格的三维编码信息
     *
     * @param gridInfo     网格信息对象，包含网格的基本信息和当前状态
     * @param level        粒度级别，用于确定网格的细分程度
     * @param zOrderMatrix Z字形矩阵，用于确定网格在二维平面上的编码
     * @param layer        高程号
     * @return 返回更新后的网格信息对象
     */
    private GridInfo calculate3DLevelCode(GridInfo gridInfo, int level, int[][] zOrderMatrix, String layer) {
        // 获取网格基本信息
        BigDecimal longitude = gridInfo.getBaseLon();
        BigDecimal latitude = gridInfo.getBaseLat();

        // 计算当前粒度级别的基准经纬度 λ₁, φ₁
        BigDecimal baseLongitude = gridInfo.getCurrentLon()
                .add(new BigDecimal(gridInfo.getColumn() - 1)
                        .multiply(gridInfo.getGridStepLon())
                        .multiply(BigDecimal.valueOf(longitude.compareTo(BigDecimal.ZERO))));

        BigDecimal baseLatitude = gridInfo.getCurrentLat()
                .add(new BigDecimal(gridInfo.getRow() - 1)
                        .multiply(gridInfo.getGridStepLat())
                        .multiply(BigDecimal.valueOf(latitude.compareTo(BigDecimal.ZERO))));

        // 设置当前网格粒度级别
        gridInfo.setGridStep(level);

        // 计算子网格列号和行号
        int column = longitude
                .subtract(baseLongitude)
                .divide(gridInfo.getGridStepLon(), 0, RoundingMode.DOWN)
                .abs()
                .intValue();

        int row = latitude
                .subtract(baseLatitude)
                .divide(gridInfo.getGridStepLat(), 0, RoundingMode.DOWN)
                .abs()
                .intValue();

        // 更新网格信息并返回
        return gridInfo.setCurrentLon(baseLongitude)
                .setCurrentLat(baseLatitude)
                .setColumn(column)
                .setRow(row)
                .setGridCode(String.format("%s%d%s",
                        gridInfo.getGridCode(),
                        zOrderMatrix[column][row],
                        layer));
    }


    /**
     * 通用的三维北斗网格编码计算方法
     *
     * @param gridInfo             网格信息对象
     * @param level                当前网格粒度级别
     * @param useMinusOneForColumn 是否使用 column - 1 计算基准经度
     * @param useMinusOneForRow    是否使用 row - 1 计算基准纬度
     * @param addOneForSubGridCalc 是否在子网格列号和行号计算时加 1
     * @param columnChars          列号字符集
     * @param rowChars             行号字符集
     * @param layer                高程号
     * @return 更新后的 GridInfo 对象
     */
    private GridInfo calculate3DLevelCode(GridInfo gridInfo,
                                          int level,
                                          boolean useMinusOneForColumn,
                                          boolean useMinusOneForRow,
                                          boolean addOneForSubGridCalc,
                                          char[] columnChars,
                                          char[] rowChars,
                                          String layer) {
        BigDecimal longitude = gridInfo.getBaseLon();
        BigDecimal latitude = gridInfo.getBaseLat();

        // 计算当前粒度级别的基准经纬度 λ₁, φ₁
        BigDecimal baseLongitude = gridInfo.getCurrentLon()
                .add(new BigDecimal(useMinusOneForColumn ? gridInfo.getColumn() - 1 : gridInfo.getColumn())
                        .multiply(gridInfo.getGridStepLon())
                        .multiply(BigDecimal.valueOf(longitude.compareTo(BigDecimal.ZERO))));

        BigDecimal baseLatitude = gridInfo.getCurrentLat()
                .add(new BigDecimal(useMinusOneForRow ? gridInfo.getRow() - 1 : gridInfo.getRow())
                        .multiply(gridInfo.getGridStepLat())
                        .multiply(BigDecimal.valueOf(latitude.compareTo(BigDecimal.ZERO))));

        // 设置当前网格粒度级别
        gridInfo.setGridStep(level);

        // 计算子网格列号和行号
        int column = longitude
                .subtract(baseLongitude)
                .divide(gridInfo.getGridStepLon(), 0, RoundingMode.DOWN)
                .abs()
                .intValue();

        int row = latitude
                .subtract(baseLatitude)
                .divide(gridInfo.getGridStepLat(), 0, RoundingMode.DOWN)
                .abs()
                .intValue();

        // 是否加 1
        if (addOneForSubGridCalc) {
            column = column + 1;
            row = row + 1;
        }

        // 更新网格信息并返回
        return gridInfo.setCurrentLon(baseLongitude)
                .setCurrentLat(baseLatitude)
                .setColumn(column)
                .setRow(row)
                .setGridCode(String.format("%s%c%c%s",
                        gridInfo.getGridCode(),
                        columnChars[column - 1],
                        rowChars[row - 1],
                        layer));
    }
}