• GIS中XYZ瓦片的加载流程解析与实现


    1. 什么是XYZ瓦片

    XYZ瓦片是一种在线地图数据格式,常见的地图底图如Google、OpenStreetMap 等互联网的瓦片地图服务,都是XYZ瓦片,严格来说是ZXY规范的地图瓦片

    ZXY规范的地图瓦片规则如下:将地图全幅显示时的图片从左上角开始,往下和往右进行切割,切割的大小默认为 256*256 像素,左上角的格网行号为 0,列号为 0,往下和往右依次递增,如下图所示:

    image

    从整体来说,XYZ瓦片数据结构是一种影像金字塔,如下图所示:

    image

    对于用户端的软件来说,所谓浏览XYZ格式的地图,就是根据当前的缩放等级和屏幕显示的地理范围,去服务端加载对应的XYZ瓦片(通常是PNG图片)

    2. XYZ瓦片与经纬度的计算以及原理

    首先给出经纬度与XYZ行列号之间的计算公式:

    image

    现在解释一下原理

    下面是一张OpenStreetMap在zoom等级为2时的瓦片示意图

    image

    z 是当前的瓦片等级,就是缩放等级,由上面的图可以看出:z 等级时,共有2z个瓦片,x范围为0-2z1,y范围也是0-2z1

    首先 x 的计算很简单:

    • 目的:将经度从-180度到180度,映射到0到2z之间的整数列号上
    • 过程:先将经度加180度,使其从0到360度,然后除以360(归一化)再乘以2z得到行号,最后向下取整数部分,得到最终的行号

    y 的计算就复杂多了:

    • 目的:将纬度从-90度到90度,映射到0到2z之间的整数行号上

    • 存在的问题:纬度分布不均匀,XYZ瓦片试图将地图展开为一个正方形(参考上图,本质上就是Web墨卡托投影),然而纬度是中间(赤道)长两极短,如果只是像 x 一样简单的映射,会导致两极的紧凑,赤道附近稀疏

    • 解决方案:将纬度通过一种映射,使其能均匀一点,然后就采用了下面的函数

      y=(1ln(tan(x)+1/cos(x))/π)2

      这个函数图像如下图所示:

    image

    • 过程:在采取上面的这个纬度的映射函数以后(归一化),再乘以2z,最后向下取整数部分,得到最终的列号

    3. 在浏览器端实现XYZ瓦片的加载示例

    3.1 计算公式实现

    根据上面的公式,很容易就把根据经纬度算行列号的函数写出来

    function lon2tile(lon, zoom) {
    return (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)))
    }
    function lat2tile(lat, zoom) {
    return (Math.floor((1 - Math.log(Math.tan(lat * Math.PI / 180) + 1 / Math.cos(lat * Math.PI / 180)) / Math.PI) / 2 * Math.pow(2, zoom)))
    }

    事实上,这个网站已经给出了这个公式的各种编程语言的实现:Slippy map tilenames - OpenStreetMap Wiki

    3.2 核心代码

    根据经纬度计算XYZ瓦片的URL,并加载到浏览器上,核心代码如下

    function lon2tile(lon, zoom) {
    return (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)));
    }
    function lat2tile(lat, zoom) {
    return (Math.floor((1 - Math.log(Math.tan(lat * Math.PI / 180) + 1 / Math.cos(lat * Math.PI / 180)) / Math.PI) / 2 * Math.pow(2, zoom)));
    }
    const loadMapByBounds = (minLon, minLat, maxLon, maxLat, zoom) => {
    const minTileX = lon2tile(minLon, zoom);
    const minTileY = lat2tile(maxLat, zoom); // Y轴是反的,自上而下
    const maxTileX = lon2tile(maxLon, zoom);
    const maxTileY = lat2tile(minLat, zoom);
    for (let x = minTileX; x <= maxTileX; x++) {
    for (let y = minTileY; y <= maxTileY; y++) {
    loadTile(x, y, zoom); // 加载瓦片
    }
    }
    }

    3.3 完整实现

    为了简单,这里使用img标签来加载瓦片图,并根据瓦片编号排列,设置对应的偏移值

    为了能拖动以浏览全图实现简单的交互,这里还设置了根据鼠标按压后拖动的偏移值来添加对应的偏移值

    实现效果如下:

    image

    完整代码如下:

    html>
    <html lang="en">
    <head>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <title>Documenttitle>
    <style>
    img {
    width: 256px;
    height: 256px;
    }
    html,
    body {
    margin: 0;
    padding: 0;
    overflow: hidden;
    height: 100%;
    width: 100%;
    }
    #map {
    position: absolute;
    height: 100%;
    width: 100%;
    overflow: hidden;
    border: 1px solid #000;
    }
    style>
    head>
    <body>
    <div id="map">div>
    <script>
    function lon2tile(lon, zoom) {
    return (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)));
    }
    function lat2tile(lat, zoom) {
    return (Math.floor((1 - Math.log(Math.tan(lat * Math.PI / 180) + 1 / Math.cos(lat * Math.PI / 180)) / Math.PI) / 2 * Math.pow(2, zoom)));
    }
    // 根据鼠标滚轮缩放地图
    let zoom = 2;
    document.body.onwheel = (e) => {
    if (e.deltaY > 0) {
    zoom--;
    } else {
    zoom++;
    }
    if (zoom < 0) {
    zoom = 0;
    return;
    }
    document.querySelector('#map').innerHTML = "";
    // EPSG:3857(Web墨卡托投影) 对应的 WGS84范围:-180.0 ,-85.06,180.0, 85.06,不在这个经纬度范围内,地图会显示异常(没有这个瓦片)
    const x1 = lon2tile(-179, zoom);
    const y2 = lat2tile(-80, zoom);
    const x2 = lon2tile(179, zoom);
    const y1 = lat2tile(80, zoom);
    const centerX = (x1 + x2) / 2;
    const centerY = (y1 + y2) / 2;
    for (let y = y1; y <= y2; y++) {
    for (let x = x1; x <= x2; x++) {
    const img = document.createElement("img");
    img.src = `https://a.tile.openstreetmap.org/${zoom}/${x}/${y}.png`;
    img.alt = `${zoom}-${x}-${y}`;
    img.style.position = "absolute";
    img.draggable = false;
    // img.style.left = `${(x - x1) * 256}px`;
    // img.style.top = `${(y - y1) * 256}px`;
    img.style.left = `${(x - centerX) * 256 + 256}px`;
    img.style.top = `${(y - centerY) * 256 + 256}px`;
    document.querySelector('#map').appendChild(img);
    }
    }
    }
    const event = new Event("wheel")
    document.body.dispatchEvent(event);
    document.body.onmousedown = (e) => {
    document.body.style.cursor = "grabbing";
    document.querySelector('#map').onmousemove = (e) => {
    // 移动地图
    const x = e.movementX;
    const y = e.movementY;
    const map = document.querySelector('#map');
    map.childNodes.forEach((img) => {
    img.style.left = `${parseInt(img.style.left) + x}px`;
    img.style.top = `${parseInt(img.style.top) + y}px`;
    });
    }
    }
    document.body.onmouseup = (e) => {
    document.body.style.cursor = "default";
    document.querySelector('#map').onmousemove = null;
    }
    script>
    body>
    html>

    4. 参考资料

    [1] Slippy map tilenames - OpenStreetMap Wiki

    [2] ZXY标准瓦片 (supermap.com.cn)

    [3] 瓦片底图:在线地图的下载和使用 | Mars3D开发教程

  • 相关阅读:
    SAP MM学习笔记38 - 入库/请求自动决济(ERS - Evaluated Receipt Settlement)
    优秀智慧园区案例 - 三亚市崖州湾科技城智慧园区,先进智慧园区建设方案经验
    CloudCompare 二次开发(16)——计算点云最值
    《动手学深度学习 Pytorch版》 4.4 模型选择、欠拟合和过拟合
    C++入门新手村-文件读写知识讲解
    【Java程序员面试专栏 算法思维】五 高频面试算法题:贪心算法
    Oracle数据泵导入和导出命令
    22-08-18 西安 尚医通(07)支付流程、rabbitmq的使用
    递归求解斐波那契数列
    在Linux系统启动java程序(jar包)
  • 原文地址:https://www.cnblogs.com/jiujiubashiyi/p/18160827