给定经纬度计算距离_栅格瓦片坐标计算

GIS日刊开篇想和大家聊聊谷歌的栅格瓦片坐标计算；瓦片地图诞生以前，电子地图通常是根据需要的地图范围和屏幕大小计算后得到的一张大图，这种模式有两个缺点：1 地图范围较大，增加服务端渲染压力，2 渲染后地图占用空间大不利于于网络传输；这样会直接导致客户端请求延迟，浏览器会长时间白屏，用户体验效果很差。瓦片地图采用分而治之的思路，每次处理一小部分，最终按照一定规则组织拼接成一张图，这样能有效提升用户体验。

Web墨卡托投影

Web墨卡托投影最早是谷歌提出的，基于墨卡托投影，将地球投影到正方形平面，投影后地图最大范围-20037508.342789244, 20037508.342789244, 20037508.342789244, -20037508.342789244,切图原点是坐标为-20037508.342789244, 20037508.34278924。目前大部分互联网地图都采用的是这个投影。下边是js实现的wgs84和web墨卡托投影转换。

//wgs84->epsg:3857const R = 6378137;const MAX_LATITUDE= 85.0511287798;class Project{ project (lonlat) { let d = Math.PI / 180, max = MAX_LATITUDE, lat = Math.max(Math.min(max, lonlat.lat), -max), sin = Math.sin(lat * d); return new Point( R * lonlat.lon * d, R * Math.log((1 + sin) / (1 - sin)) / 2); } unproject (point) { var d = 180 / Math.PI; return new LonLat(         point.x * d / R,         (2 * Math.atan(Math.exp(point.y / R)) - (Math.PI / 2)) * d); }}

                  瓦片组织规则        Web墨卡托投影的瓦片以金字塔形式分布，顶部为地球全貌(两极部分除外)投影到256*256大小的0级瓦片，第1级投影到4张256*256大小的瓦片，依次类推，金字塔层级和瓦片总数呈4^n分布，行列数呈2^n分布，如图所示，图中方格中心标注即为瓦片坐标。                  瓦片坐标解析         想要准确展现一幅完整的地图，需要将地图范围内瓦片按照上述的规则计算相应坐标并拼接成图。下边以leaflet.js引擎加载openstreetmap tms瓦片服务为例，剖析地图瓦片坐标计算。做过gis开发的同学一定熟悉这个方法setCenter或者setView，就是设置地图层级和中心点的方法。但是这个方法做了什么事情呢？首先是计算中心瓦片坐标，其次计算视图范围内所有瓦片坐标，然后将瓦片坐标换算到屏幕坐标。        解析瓦片坐标前先介绍两个概念，resolution每像素的地图单位(本文中可以理解为每像素代表多少米)，origin切图原点从投影后地图什么坐标开始切图(-20037508.342789244,20037508.34278924)。        map.setView([39.85, 116.3], 10);  将地图层级设置为10级，中心点定位到39.85(纬度),116.3(经度)，接下来开始计算地图中心瓦片坐标。         1 计算resolution        开发过程中resolution是不需要计算的，一般情况下，瓦片服务会给定每一个对应的resolution值；为了大家方便理解这个概念，在此计算一下web墨卡托投影瓦片第10级的resolution。首先计算10级有多少行或多少列瓦片，由于投影后地图为正方形，全幅地图行数和列数是相等的，Web墨卡托投影计算简单也体现在这；其次计算全幅多少像素；然后计算全幅多少米，最后得到  resolution = 全幅米 / 全幅像素。       列数cols = Math.pow(2,10) = 1024;       全幅像素pixels = cols* 256 = 262144;       全幅米 meters = 20037508.342789244+ 20037508.342789244 = 40075016.68557849;       resolution  = meters / pixels = 152.8740565703525;         2 计算经纬度转换到web墨卡托投影系的坐标xy = new Project().project({
       lon: 116.3,       lat: 39.85}) ;结果为：xy = { x: 12946456.779257717, y: 4844168.570470275}         3 计算中心瓦片横坐标       当前位置x方向距离原点多少米mxs = 12946456.779257717– (-20037508.342789244) = 32983965.122046962;       当前位置x方向距离原点多少像素pxs = mxs/ resolution = 215759.07555555555;       瓦片横坐标X = Math.floor(pxs/ 256) = 842;         4 计算中心瓦片纵坐标       当前位置y方向距离原点多少米mys = 20037508.342789244- 4844168.570470275 = 15193339.772318969;       当前位置y方向距离原点多少像素pys = mys / resolution = 99384.68379248513;       瓦片纵坐标Y= Math.floor(pys/ 256) = 388;    计算后中心瓦片坐标为：    {        x:  842,        y:  3 88    }