「GIS算法」Bing Maps中的Quadkeys算法解析

小编最近在研究 ArcGIS Pro 加载 Bing Maps 的时候，对 Bing Maps 进行了简单的抓包，结果发现一个有点反常的现象，Bing Maps 在请求瓦片数据的时候其请求格式与常见的有明显不同。常见的比如高德地图这种会传一个xyz参数，这个比较好理解就是对应xy坐标和绽放层级z。而Bing Maps则传了一个整数参数，并没有出来xyz三个数，这是怎么一回事呢？下图是Bing Maps请求遥感影像的时候的请求抓包。

下图是高德地图请求的请求抓包。

很显然，Bing Maps的请求格式与高德地图的请求格式有区别，Bing Maps的请求格式是整数参数，而高德地图的请求格式是xyz参数。那么Bing Maps的请求格式是怎么生成的呢？

Quadkeys算法原理

说到这不得不说Bing Maps的 Quadkeys算法了，根据官方文档(https://learn.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system?redirectedfrom=MSDN)的介绍，其核心将之前的xyz三维参数转换成了一维参数，这一点和GoeHash的有点像。那么是怎么转换的呢？

我们先以官方的这张图为例：

假设我们要表示小编出的这个点位，使用xyz的方式的话应该是（3，5，3），即x=3，y=5，z=3。使用 Quadkeys算法的话，应该如何表示呢？

第一步：将x y转换成二进制数，先将x转换成二进制数：011，为什么是011不是11？这里要保持其位数与z相等，不够使用0补全，再将y表示成二进制数：101。

第二步：交叉编码，先y后x，得到编码 10011

第三步：将编码 10011 转换成4进制数，得到编码 213。

Quadkeys算法特点

Quadkeys算法有三个重要特征：

1. Quadkey的长度等于瓦片对应的zoom等级，比如这个213的长度等于3，即对应zoom等级为3。
2. Quadkey的字符串的开头是父瓦片对应编码，还是看官方那个图，瓦片 2 是瓦片 20 到 23 的父瓦片，瓦片 13 是瓦片 130 到 133 的父瓦片。
3. 相邻 XY 坐标的两个瓦片通常具有相对靠近的编码，这对于优化数据库性能很重要，因为通常批量形式请求相邻的瓦片，将这些瓦片存储在相同的磁盘块上，可以以最大限度地减少磁盘读取次数。

Quadkeys算法实现

理解了算法本身就可以了，目前官方给了C#实现，小编使用AI实现了一个Python版本，代码如下：

import math

class TileSystem:
    EARTH_RADIUS = 6378137
    MIN_LATITUDE = -85.05112878
    MAX_LATITUDE = 85.05112878
    MIN_LONGITUDE = -180
    MAX_LONGITUDE = 180

    @staticmethod
    def clip(n: float, min_value: float, max_value: float) -> float:
        """Clips a number to the specified minimum and maximum values."""
        return min(max(n, min_value), max_value)

    @staticmethod
    def map_size(level_of_detail: int) -> int:
        """Determines the map width and height (in pixels) at a specified level of detail."""
        return 256 << level_of_detail

    @staticmethod
    def ground_resolution(latitude: float, level_of_detail: int) -> float:
        """Determines the ground resolution (in meters per pixel) at a specified latitude and level of detail."""
        latitude = TileSystem.clip(latitude, TileSystem.MIN_LATITUDE, TileSystem.MAX_LATITUDE)
        return (math.cos(latitude * math.pi / 180) *
                2 * math.pi * TileSystem.EARTH_RADIUS / TileSystem.map_size(level_of_detail))

    @staticmethod
    def map_scale(latitude: float, level_of_detail: int, screen_dpi: int) -> float:
        """Determines the map scale at a specified latitude, level of detail, and screen resolution."""
        return TileSystem.ground_resolution(latitude, level_of_detail) * screen_dpi / 0.0254

    @staticmethod
    def latlong_to_pixel_xy(latitude: float, longitude: float, level_of_detail: int):
        """Converts latitude/longitude WGS-84 coordinates into pixel XY coordinates."""
        latitude = TileSystem.clip(latitude, TileSystem.MIN_LATITUDE, TileSystem.MAX_LATITUDE)
        longitude = TileSystem.clip(longitude, TileSystem.MIN_LONGITUDE, TileSystem.MAX_LONGITUDE)

        x = (longitude + 180) / 360
        sin_latitude = math.sin(latitude * math.pi / 180)
        y = 0.5 - math.log((1 + sin_latitude) / (1 - sin_latitude)) / (4 * math.pi)

        map_size = TileSystem.map_size(level_of_detail)
        pixel_x = int(TileSystem.clip(x * map_size + 0.5, 0, map_size - 1))
        pixel_y = int(TileSystem.clip(y * map_size + 0.5, 0, map_size - 1))

        return pixel_x, pixel_y

    @staticmethod
    def pixel_xy_to_latlong(pixel_x: int, pixel_y: int, level_of_detail: int):
        """Converts pixel XY coordinates into latitude/longitude WGS-84 coordinates."""
        map_size = float(TileSystem.map_size(level_of_detail))
        x = (TileSystem.clip(pixel_x, 0, map_size - 1) / map_size) - 0.5
        y = 0.5 - (TileSystem.clip(pixel_y, 0, map_size - 1) / map_size)

        latitude = 90 - 360 * math.atan(math.exp(-y * 2 * math.pi)) / math.pi
        longitude = 360 * x

        return latitude, longitude

    @staticmethod
    def pixel_xy_to_tile_xy(pixel_x: int, pixel_y: int):
        """Converts pixel XY coordinates into tile XY coordinates."""
        tile_x = pixel_x // 256
        tile_y = pixel_y // 256
        return tile_x, tile_y

    @staticmethod
    def tile_xy_to_pixel_xy(tile_x: int, tile_y: int):
        """Converts tile XY coordinates into pixel XY coordinates."""
        pixel_x = tile_x * 256
        pixel_y = tile_y * 256
        return pixel_x, pixel_y

    @staticmethod
    def tile_xy_to_quadkey(tile_x: int, tile_y: int, level_of_detail: int) -> str:
        """Converts tile XY coordinates into a QuadKey."""
        quadkey = []
        for i in range(level_of_detail, 0, -1):
            digit = 0
            mask = 1 << (i - 1)
            if (tile_x & mask) != 0:
                digit += 1
            if (tile_y & mask) != 0:
                digit += 2
            quadkey.append(str(digit))
        return ''.join(quadkey)

    @staticmethod
    def quadkey_to_tile_xy(quadkey: str):
        """Converts a QuadKey into tile XY coordinates and level of detail."""
        tile_x = tile_y = 0
        level_of_detail = len(quadkey)
        for i in range(level_of_detail, 0, -1):
            mask = 1 << (i - 1)
            digit = quadkey[level_of_detail - i]
            if digit == '1':
                tile_x |= mask
            elif digit == '2':
                tile_y |= mask
            elif digit == '3':
                tile_x |= mask
                tile_y |= mask
            elif digit != '0':
                raise ValueError("Invalid QuadKey digit sequence.")
        return tile_x, tile_y, level_of_detail

使用方法：

if __name__ == "__main__":
    lat, lon = 47.60357, -122.32945  # 西雅图
    level = 15

    px, py = TileSystem.latlong_to_pixel_xy(lat, lon, level)
    tx, ty = TileSystem.pixel_xy_to_tile_xy(px, py)
    quadkey = TileSystem.tile_xy_to_quadkey(tx, ty, level)

    print(f"LatLon=({lat}, {lon})")
    print(f"Pixel=({px}, {py})  Tile=({tx}, {ty})")
    print(f"QuadKey={quadkey}")

    # 反向转换验证
    txx, tyy, lvl = TileSystem.quadkey_to_tile_xy(quadkey)
    print(f"QuadKey -> Tile=({txx}, {tyy}), Level={lvl}")

PS：上述代码由AI生成，请自行测试，需要其他版本的可以自行使用AI来实现。

总结

其实通过上述的总结可以看出这个 QuadKey 算法，其实和xyz这种算法是可以相互转换，这也是小编之前的文章《ArcGIS Pro 添加 Bing Maps 底图 》中能够实现功能核心因素。

各位在项目中有用到 QuadKey 的吗？欢迎分享更多的使用经验。