2023-02-19--------瓦片投影与行列号转换一、原理墨卡托适用于二维，从0级开始，为正方形，列号从左往右，

一、原理

1.1 经纬度投影（4326）

影像目前只有天地图提供4326瓦片

特点：从第1级开始，左右各一张,为长方形，列号从左往右，行号从上到下，每级瓦片:

WGS84经纬度投影瓦片的第0级无法构成一张可以将一个256x256图片填充满的合法瓦片，因此第0级为无效瓦片

行数=Math.pow(2,zoom-1) 列数=Math.pow(2.zoom)

第1级

第2级

1.2 Web墨卡托投影（3857）

特点：从0级开始，为正方形，每级瓦片:

行数=列数=Math.pow(2,zoom)

1.2.1 TMS瓦片（从下往上）

TMS(Tile Map Service) 是 OSGeo (开源地理基金会) 提出的一种地图瓦片服务。额外补充一句，WMTS、WMS、WFS这些是 OGC（开放地理空间信息联盟）提出的。

TMS 地图瓦片有如下特点：

瓦片编号从左下角开始，x轴为 -85.0511° 纬线，y轴为 180° 经线，第一个瓦片编号为 (0, 0)；
x 轴编号（行号）从左到右依次递增，y 轴编号（列号）从下到上依次递增；
单个地图瓦片为 256x256 大小的正方形图片。

此外，地图缩放等级 z 和每列（或每行）瓦片的数量(记为 n )有如下关系：

TMS 的优点是地图瓦片可存放在本地，类似于静态文件，可使用 Nginx 等 Web 服务器直接代理，然后通过一定规则进行访问，便于地图瓦片服务离线化部署。

TMS 通常采用类似于如下 url 进行访问：

http://xxx/xxx/{z}/{x}/{y}.png

1.2.2 Google/Bing/OSM/ESRI 地图瓦片（从上往下）

这种地图瓦片的组织方式为：

原点在左上角
x轴在 85.0511° 纬线，y轴为 180° 经线；
y轴编号从上到下递增，其他特点和 TMS 无异。

简单说就是 y 轴方向和 TMS 相反

加载这类瓦片需要预先将所有瓦片的 y 轴编号转置一下，然后再加载，转置公式如下：

y_2=2^z-y_1-1

1.3 cesium中的地形划分

目前标准的Cesium地形特点：

从0(不同于影像从1级开始)级开始，划分为左右两张；
列号从左往右，行号从下到上

请求url一般为以下形式： http://xxx/xxx/{zoom-1}/{x}/{y}.terrain

二、js代码

2.1 经纬度转墨卡托

function lonlatToMercator(lonlat){
    var mercator = {};
    
    let x = lonLat.lon * 20037508.34 / 180;
    let y = Math.log(Math.tan((90 + lonLat.lat) * Math.PI / 360)) / (Math.PI / 180);
    
    y = y * 20037508.34 / 180;
    
    mercator.x = x;
    mercator.y = y;
    
    return mercator;
}

2.2 墨卡托转经纬度

function mercatorTolonlat(mercator){
    let lonlat={lon:0,lat:0};
    
    let x = mercator.x/20037508.34*180;
    let y = mercator.y/20037508.34*180;
    
    y= 180/Math.PI*(2*Math.atan(Math.exp(y*Math.PI/180))-Math.PI/2);
    
    lonlat.lon = x;
    lonlat.lat = y;

    return lonlat;
};

2.3 经纬度投影下，经纬度和瓦片编号互转

// 经度转瓦片列号（适用于4326和3857，因为经度是均匀的，因此行号一致）
function lon2tile(lon, zoom) {
    return (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)));
}

// 瓦片列号转经度
function tile2long(x, z) {
    return (x / Math.pow(2, z) * 360 - 180);
}



//纬度转瓦片行号（4326）
function lat2tile4326(lat, zoom) {
    return (Math.floor((90 - lat) / 180 * Math.pow(2, zoom)));
}
//瓦片行号转纬度
function tileGeo2Lat(row, zoom) {
    return 90-180*row/Math.pow(2, zoom);
}

2.3 web墨卡托行列号与经纬度互转

// 经度转瓦片列号（适用于4326和3857，因为经度是均匀的，因此行号一致）
function lon2tile(lon, zoom) {
    return (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)));
}

// 瓦片列号转经度
function tile2long(x, z) {
    return (x / Math.pow(2, z) * 360 - 180);
}



//纬度转瓦片行号（3857）,纬度不均匀，因此需要先转为web墨卡托再计算
function lat2tile3857(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)));
}

//3857下的行号转纬度
function tileWeb2lat(y, z) { 
    var n = Math.PI - 2 * Math.PI * y / Math.pow(2, z); 
    return (180 / Math.PI * Math.atan(0.5 * (Math.exp(n) - Math.exp(-n)))); 
}



/**
 * web墨卡托坐标转瓦片行列号
 * @param {*} zoom 
 * @param {*} x 
 * @param {*} y 
 * @returns 
 */
function web2XY(zoom, x, y) {
    //计算总共多少列瓦片
    const cols = Math.pow(2, zoom + 1)
    //计算总像素值
    const pixels = cols * 256;
    //计算总米数
    const meters = 20037508.342789244 + 20037508.342789244
    //计算分辨率
    const resolution = meters / pixels;

    //当前位置x方向距离原点多少米
    const rangeX = x - (-20037508.342789244);
    //当前位置x方向距离原点多少像素
    const pxs = rangeX / resolution;

    //
    //瓦片横坐标
    const i = Math.floor(pxs / 256);


    //当前位置y方向距离原点多少米
    const rangeY = 20037508.342789244 - y;

    //当前位置y方向距离原点多少像素
    const pys = rangeY / resolution;
    //瓦片纵坐标
    const j = Math.floor(pys / 256)
    return { i, j }


}

//计算经纬度范围对应的4326瓦片编号
function computeTileArea(bounds, zoom) {

    const startCol = lon2tile(bounds.left, zoom + 1)

    const endCol = lon2tile(bounds.right, zoom + 1)
    const startRow = lat2tile4326(bounds.top, zoom)
    const endRow = lat2tile4326(bounds.bottom, zoom)
    return {
        startCol, endCol, startRow, endRow
    }
}