地球空间网格编码正反算实现

三维网格码计算思路分享。

似乎这两年低空经济成为了潮流方向，多了一些政策倾斜。前阵子刚有公司来交流相关事宜，16日华师大又要新成立低空经济研究中心，这两天组里为了这些事情鞍前马后，竟然有些仿佛过年的热闹。我负责的工作是探索其中一小个部分，也是一个非常小众的应用方向：地球空间网格码，该工作目前的理论基础为国标GBT 40087-2021 地球空间网格编码规则。
地球空间网格基于 GeoSOT(Geographical coordinate global Subdivision based on One-dimensioninteger and Two ton"power)地球剖分模型,将地球空间统一剖分成不同尺度的网格单元,并按统一编码规则进行标识和表达,构建了网格化的地球空间数据组织参考框架。该框架支持地球表面空间和地球立体空间与地理空间信息的聚合,可有效解决物联网、大数据、云计算中海量空间信息在标识和表达上的唯一性、可读性、尺度性、关联性的瓶颈,实现了多源,多尺度数据网格化高效组织、处理和应用,突破了地理空间信息跨行业应用的技术壁垒,推动地球系统科学的发展。
具体来说，网格剖分分为经纬度的水平面剖分和高程的垂直剖分。具体过程相对复杂，参见上述国家标准，本文仅分享一些实现步骤。

1. 经纬度网格剖分

根据国标，经纬度剖分的步骤如下所示：

1. 将该点经、纬度坐标表示成GB/T16831规定的形式,即A’B’C.D";
2. 将该点的坐标按度、分、秒、秒小数部分分别转换为二进制数。即将度|A|由十进制数转换成8 bit 定长二进制数(A)2, 将分B由一进制数转换成6bit定长二进制数(B)2, 将秒C由十进制数转换成6bit定长二进制数©2, 将秒以下数D由一进制数转换成11bit定长二进制数(D)2;
3. 分别将经、纬度坐标度、分、秒及秒以下二进制数(A)2、(B)2、©2、(D)2直接拼接成31 bit定长二进制数(E)2, 即(E)2=(A)2(B)2©2(D)2, 分别得到两个 31 bit定长数经度(EL)2, 和纬度(EB)2;
4. 将纬度(EB)2,前置、经度(EL)2后置, 采用莫顿交叉的方式生成62bit的混合代码(F)2, 例如若(EB)2为100111 ,(EL)2为011010, 则(EB)2在前, (EL)2在后的莫顿交叉运算结果为(F)2为 100101101110;
5. 将二进制混合代码(F)2转成四进制编码(F)4;
6. 根据待求网格的级别m, 将(F)4中后32-m位四进制的码元去掉得到(F’)4;
7. 根据经度和纬度,前加上G0、G1、G2或G3即可得到网格代码。

代码实现步骤如下：
定义一个类：

class WGS84coordinate:
    def __init__(self, latitude, longitude, elevation):
        self.latitude = latitude
        self.longitude = longitude
        self.elevation = elevation

实现步骤1：

# 1. 将该点经、纬度坐标表示成GB/T16831规定的形式,即A'B'C.D";
    def decimal2dms(self, decimal):
        degrees = int(decimal)
        minutes = int((abs(decimal) - abs(degrees)) * 60)
        seconds = (abs(decimal) - abs(degrees)) * 60
        seconds = (seconds - minutes) * 60
        seconds_fractional, seconds_integer = math.modf(seconds)
        return [int(degrees), int(minutes), int(seconds_integer), seconds_fractional]

实现步骤2,3：

# 2. 将该点的坐标按度、分、秒、秒小数部分分别转换为二进制数, 拼为31 bit定长数
    def dms2bit(self, dms):
        # dms = [degrees, seconds_fractional, seconds_integer]
        degrees, minutes, seconds_integer, seconds_fractional = dms
        # print(f"度：{degrees}, 分：{minutes}, 秒整数：{seconds_integer}，秒小数：{seconds_fractional}")
        degrees_bin, minutes_bin, seconds_integer_bin = f"{degrees:08b}", f"{minutes:06b}", f"{seconds_integer:06b}"
        seconds_fractional_scaled = int(seconds_fractional * 2048)
        seconds_fractional_bin = f"{seconds_fractional_scaled:011b}"
        # print(f"度：{degrees_bin}, 分：{minutes_bin}, 秒整数：{seconds_integer_bin}，秒小数：{seconds_fractional_bin}")
        return ''.join([degrees_bin, minutes_bin, seconds_integer_bin, seconds_fractional_bin])

实现步骤4,5：

# 3. 将经纬度的编码莫顿交叉,转为四进制编码
    def mortoncross(self):
        latbit = self.dms2bit(self.decimal2dms(self.latitude))
        lonbit = self.dms2bit(self.decimal2dms(self.longitude))
        
        assert len(latbit) == 31 and len(lonbit) == 31, "Both binary strings must be 31 bits long."
        interleaved_bits = []
        for i in range(31):
            interleaved_bits.append(latbit[i])
            interleaved_bits.append(lonbit[i])
        bit2 = ''.join(interleaved_bits)
        bit4 = ''.join(str(int(bit2[i:i+2], 2)) for i in range(0, len(bit2), 2))

        return bit4

实现步骤6,7：

# 4. 根据待求网格的级别m,将(F)中后32-m位四进制的码元去掉得到(F);
    def gridcode(self, m):
        bit4 = self.mortoncross()
        gridcode = bit4[:m-32]

        gridcode = "G0" + gridcode
        # 高度码实现下面讲。
        elevationcode = self.elevationcode(m)

        return [gridcode, elevationcode]

2. 高程网格剖分

$$n = \frac{\theta_0}{\theta} \log_{1 + \theta_0} \left( \frac{H + r_0}{r_0} \right) \tag{5}$$

其中：

• ( n ) —— 从地面向上（或向下）数第n层立体网格，n为整数，地面以上n大于等于0，地面以下 n < 0；
• $\theta_0$—— 初始划分范围定义的基础网格（1°网格）对应的经（纬）跨度差，单位为弧度（rad）
  $\theta_0 = \frac{\pi}{180}$
• $\theta$ —— 该网格对应的经（纬）跨度差，单位为弧度（rad）；
• $r_0$ —— 地球长半轴，单位为米（m），$r_0 = 6\,378\,137$；
  去掉层数$n$ 高位的 0 形成的编码即为高度域编码。

# 计算高度的网格码
    def elevationcode(self, m):
        theta_0 = math.pi / 180
        theta = GridSpanlist[m-1] * math.pi / 180
        H = self.elevation
        r_0 = 6378137
  
        # 计算公式
        n = (theta_0 / theta) * math.log((H + r_0) / r_0, (1 + theta_0))
        n = math.floor(n)
        if n < 0:
            n = 0
        n = str(bin(n)[2:])

        if len(n) < m:
            n = '0' * (m - len(n)) + n
        n = "H" + n

        return n

3. 网格码反算

这一部分的具体求解国标并没有提及，因此只能尝试自行求解，代码准确性为止，仅供参考：

class GridCode:
    def __init__(self, grid_code, elevation_code, m):
        self.grid_code = grid_code
        self.elevation_code = elevation_code
        self.m = m

    def __str__(self):
        print(f"格网编码：{self.grid_code, self.elevation_code}")

    def degridcode_bit2dms(self, binary_str):
        degrees_bin = binary_str[0:8]    # 前8位是度
        minutes_bin = binary_str[8:14]   # 接下来的6位是分
        seconds_integer_bin = binary_str[14:20]  # 接下来的6位是秒的整数部分
        seconds_fractional_bin = binary_str[20:31]  # 最后的11位是秒的小数部分

        # 将每个二进制子串转换回整数
        degrees = int(degrees_bin, 2)
        minutes = int(minutes_bin, 2)
        seconds_integer = int(seconds_integer_bin, 2)
        seconds_fractional_scaled = int(seconds_fractional_bin, 2)

        # 将秒的小数部分恢复为原来的小数值
        seconds_fractional = seconds_fractional_scaled / 2048
        total_seconds = seconds_integer + seconds_fractional

        # 将秒数转换为分的部分（秒/60）
        total_minutes = minutes + total_seconds / 60

        # 将分数转换为度的部分（分/60）
        decimal = abs(degrees) + total_minutes / 60

        return decimal

    # 编码反算
    def degridcode(self):
        gridcode = self.grid_code
        gridcode = gridcode[2:]
        gridcode = gridcode.ljust(31, '0')

        bit2 = ''.join(bin(int(c, 4))[2:].zfill(2) for c in gridcode)

        latbit = []
        lonbit = []
        for i in range(1, 32):
            latbit.append(bit2[i * 2 - 2])
            lonbit.append(bit2[i * 2 - 1])

        latcode = ''.join(latbit)
        loncode = ''.join(lonbit)

        lat = self.degridcode_bit2dms(latcode)
        lon = self.degridcode_bit2dms(loncode)
        
        return lat, lon

    # 高度码反算
    def deelevationcode(self):
        n_binary, m = self.elevation_code, self.m

        theta_0 = math.pi / 180
        theta = GridSpanlist[m-1] * math.pi / 180
        r_0 = 6378137

        # 高度码反算
        n_binary = n_binary[1:]

        # 将二进制的 n 转换回整数
        n = int(n_binary, 2)
        
        # 根据公式反向解码，恢复 H
        H = r_0 * (math.exp((n * theta) / theta_0 * math.log(1 + theta_0)) - 1)
        
        return H
        
    # 整体反算：
    def decode(self):
        lat, lon = self.degridcode()
        H = self.deelevationcode()
        return lat, lon, H

这样做的目的是为了得到网格码的代表性坐标，但是这个坐标到底具体代表了网格码六面体的哪个位置，不得而知。本文仅代表个人的一些探索，代码并不一定正确。