TileID编码与TileID解码莫顿码NDS与WGS84坐标转换（Python、C++双实现）

NDS与WGS84坐标转换、Tile计算

本文出现的知识点

以上知识点点击链接 Partitioning of Geographic
Data（NDS，导航数据标准中的地理数据分区）
本文章是对以上的一些细节补充，以及代码实现。

莫顿码（Morton Code）

最容易疑惑的就是莫顿码。

莫顿码是将多维数据转化为一维数据的编码。
莫顿编码定义了一条 Z 形的空间填充曲线，因此莫顿编码通常也称Z阶曲线(Z-order curve)。

Morton Code 二维数据（x,y）的编码方式：

在 N 维空间中对于彼此接近的坐标具有彼此接近的莫顿码, 可以应用于为一个整数对产生一个唯一索引。例如，对于坐标系中的坐标点使用莫顿编码生成的莫顿码，可以唯一索引对应的点。
这些索引为“Z”形排序 。
如下图以Z形(左下->右下->左上->右上)分别代表 2×2、4×4、… n² 平方单位：

Z型排列，以及编码过程：

Tile切割规则、Tile Level编码规则

Tile id 32位各个Level的存储规则

一个Tile ID由32位bit组成。

根据上图的对照变，可以看出 红色的位存储Tile等级，蓝色位不存储，灰色位存储 NDS 横、纵坐标的莫顿码。

WGS84 坐标转化为Tile的过程

我们演示将NDS坐标准换成莫顿码的过程。

比如我们有一个WGS84坐标 （x: 121.00902, y: 30.88306）
根据公式讲WGS84 坐标转化为NDS坐标

nds = wgs84 * 2³²/ 360
2³² = 1 << 32 = 4294967296

nds_x = int(121.00902 * 4294967296 / 360) 
nds_y = int( 30.88306 * 4294967296 / 360)

nds_x = 1443693842 = 01010110000011010000010100010010 (二进制)
nds_y = 368449257  = 00010101111101100001011011101001

此时NDS坐标为32位坐标，
在转化为Tile ID时需对NDS坐标做位舍弃。

比如我们需要 6级Tile ID:

前 16 – level 位 存储tile等级
中间 15 – level 位存储 0 ，占位
后 2 * level + 1 位 DNS坐标的Morton Code

根据 6级 Tile ID编码规则 ： 10位（tile等级 ） + 9位（中间占位）+13位（nds坐标的Morton Code）

实际上 在13位Morton位中，X坐标占用7位，Y坐标占用6位。

那么nds位为什么不是14位，x，y坐标各占用七位？实现偶数位？
解释一下为啥 Tile ID中的 Morton Code不是 x y 不是对称存储： 因为在0级tile编码中，只有x轴分割出两段0,1，y 都是0，编码ID为后 00,01，舍去y占位，一级Tile就只有0,1两个ID。而后续无论1级还是13级Tile ，都是在该基础上切割，所以y都舍弃最前编码的0。

以下是对NDS坐标Morton编码过程：
使用nds坐标从头到尾编码顺序；

nds_x 选取7位
nds_x = 1443693842 = 0101011 0000011010000010100010010
nds_y 选取7位,舍弃头部一位(因为它永远是0)
nds_y = 368449257 = 0 001010 1111101100001011011101001

6级Tile编码：10位（tile等级 ） + 9位（中间占位）+13位（nds坐标的Morton Code）

坐标信息

NDS坐标与WGS84坐标转换

C++ 实现

int tile_level;
double wgs84_x,wgs84_y 
auto nds_x = std::floor(wgs84_x* 4294967296 / 360) >> (31 - tile_level);
auto nds_y = std::floor(wgs84_y* 4294967296 / 360) >> (31 - tile_level);

Python 实现

def wgs84_to_nds(wgs84_x, wgs84_y, tile_level):
    """
    WGS84 坐标转NDS坐标
    :param wgs84_x:
    :param wgs84_y:
    :param tile_level:
    :return: 
    """
    nds_x = int(wgs84_x * 4294967296 / 360) >> (31 - tile_level)
    nds_y = int(wgs84_y * 4294967296 / 360) >> (31 - tile_level)
    return nds_x, nds_y

Tile ID转NDS计算

C++ 实现

/**
 * 解析 TileId 为NDS坐标系的横轴坐标
 * @param tile_id 输入tileid
 * @param x 横坐标
 * @param y 纵坐标
 * @date  2024/10/17 18:58
 */
void ParseTileId(const int32_t tile_id, int32_t& x, int32_t& y)
{
    x = 0;
    y = 0;
    int32_t morton_code = tile_id;
    auto current_level = ParseTileLevel(tile_id);
    for (int32_t i = 0; i < current_level; i++)
    {
        auto move_bit = (i * 2);
        x |= (morton_code & 1 << move_bit) >> i;
        y |= (morton_code & 1 << (move_bit + 1)) >> (i + 1);
    }
}

Python 实现

def parse_tile_id_2_nds(tile_id):
    """
    解析TileID到NDS坐标
    :param tile_id
    :return: nds_x, nds_y
    """
    nds_x = 0
    nds_y = 0
    morton_code = tile_id
    current_level = parse_tile_level(tile_id)
    for i in range(current_level):
        move_bit = (i * 2)
        nds_x |= (morton_code & 1 << move_bit) >> i
        nds_y |= (morton_code & 1 << (move_bit + 1)) >> (i + 1)
    return nds_x, nds_y

NDS、WGS84坐标转Tile ID计算

C++ 实现

int32_t EncodeTileId(const double x, const double y, const int32_t level, const bool is_wgs84)
{
    auto nds_x = 0;
    auto nds_y = 0;

    if (is_wgs84)
    {
        auto _nds_x = std::floor(x * 4294967296 / 360);
        auto _nds_y = std::floor(y * 4294967296 / 360);
        nds_x = static_cast<int32_t>(_nds_x);
        nds_y = static_cast<int32_t>(_nds_y);
        nds_x = nds_x >> (31 - level);
        nds_y = nds_y >> (31 - level);
    }
    else
    {
        nds_x = static_cast<int32_t>(x);
        nds_y = static_cast<int32_t>(y);
    }

    int32_t tile_id = (1 << (level + 16));
    for (int32_t i = 0; i < level; i++)
    {
        tile_id |= (((nds_x & (1 << i)) >> i) << (2 * i));
        tile_id |= (((nds_y & (1 << i)) >> i) << (2 * i + 1));
    }
    return tile_id;
}

Python 实现

def encode_tile_id(x, y, level=13, is_wgs84=False):
    """
    NDS、WGS84坐标转Tile ID计算
    :param x: 坐标x
    :param y: 坐标y
    :param level: tile 等级
    :param is_wgs84: True wgs84坐标
    :return: Tile ID
    """
    if is_wgs84:
        nds_x = int(x * 4294967296 / 360) >> (31 - level)
        nds_y = int(y * 4294967296 / 360) >> (31 - level)
    else:
        nds_x = int(x)
        nds_y = int(y)
    tile_id = (1 << (level + 16))
    # Mordon code calculate 
    for i in range(level):
        tile_id |= (((nds_x & (1 << i)) >> i) << (2 * i))
        tile_id |= (((nds_y & (1 << i)) >> i) << (2 * i + 1))
    return tile_id

Tile ID 等级解析

C++ 实现

/**
 * 解析 TileId 等级 （0-15级）
 * @param tile_id
 * @return level
 */
int32_t ParseTileLevel(const int32_t tile_id)
{
    int32_t level = 15;
    while (true)
    {
        auto move_bit = level + 16;
        if (((tile_id & (1 << move_bit)) >> move_bit) == 1)
        {
            return level;
        }
        level--;
    }
}

Python 实现

def parse_tile_level(tile_id):
    """
    解析Tile level
    :param tile_id:
    :return:
    """
    level = 15
    while True:
        move_bit = level + 16
        if ((tile_id & (1 << move_bit)) >> move_bit) == 1:
            return level
        level -= 1

根据Tile ID获取周边相邻的8个Tile ID

C++ 实现

/**
 * 根据Tile Id获取相邻 Tile集合
 * @details 返回该tile的上下左右四个方向的tile 以及四个方向夹角的所有tile
 * @param tile_id
 * @param tiles
 */
void GetAdjacentTiles(const int32_t tile_id, std::vector<int32_t>& tiles)
{
    // 计算 tile 等级
    auto level = ParseTileLevel(tile_id);
    // 定义 nds 横纵坐标
    int32_t x, y;
    ParseTileId(tile_id, x, y);
    // 顺序不可变
    tiles.push_back(EncodeTileId(x - 1, y - 1, level)); // 左上
    tiles.push_back(EncodeTileId(x - 1, y, level)); // 上
    tiles.push_back(EncodeTileId(x - 1, y + 1, level)); // 右上
    tiles.push_back(EncodeTileId(x, y + 1, level)); // 右
    tiles.push_back(EncodeTileId(x + 1, y + 1, level)); // 右下
    tiles.push_back(EncodeTileId(x + 1, y, level)); //下
    tiles.push_back(EncodeTileId(x + 1, y - 1, level)); // 左下
    tiles.push_back(EncodeTileId(x, y - 1, level)); // 左
}

Python 实现

def get_adjacent_tiles(tile_id):
    """
    获取周边8和tile id集合
    :param tile_id:
    :return:
    """
    _tiles = []
    level = parse_tile_level(tile_id)
    nds_x, ndx_y = parse_tile_id_2_nds(tile_id)
    # 顺序不可变
    _tiles.append(encode_tile_id(nds_x - 1, ndx_y - 1, level))  # 左上
    _tiles.append(encode_tile_id(nds_x - 1, ndx_y, level))  # 上
    _tiles.append(encode_tile_id(nds_x - 1, ndx_y + 1, level))  # 右上
    _tiles.append(encode_tile_id(nds_x, ndx_y + 1, level))  # 右
    _tiles.append(encode_tile_id(nds_x + 1, ndx_y + 1, level))  # 右下
    _tiles.append(encode_tile_id(nds_x + 1, ndx_y, level))  # 下
    _tiles.append(encode_tile_id(nds_x + 1, ndx_y - 1, level))  # 左下
    _tiles.append(encode_tile_id(nds_x, ndx_y - 1, level))  # 左
    return _tiles

作者：DAKANGGO