深入浅出学习切片LOD——ArcGIS server模拟缓存切片（影像快显）

一、第一次实践

原理

免切片实现影像服务的模拟切片，主要原理是接收前端传过来的xyz(行列层级)以及切片方案，计算出该请求的切片的四至经纬度信息，通过mapserver的exportImage接口，传入每个模拟切片的四至经纬度信息得到图片返回。好似前端请求了http://xxx/mapserver/z/y/x,后台返回给前端一张图片（切片），让客户端无感知的就像是请求了缓存切片。

而这问题的关键在于获取目标图片的四至经纬度信息。怎么计算？

首先，咱们需要知道细节层次模型LOD。LOD技术即Levels of Detail的简称，意为多细节层次（百度解释的描述会越描越黑）。实际上就是世界地图展示的时候想看某个细节，那么就需要缩放（因为屏幕大小不会变），但缩放不能直接把整个图缩放，数据细节太大，放不下。那只需要把感兴趣的细节部分数据返回给客户拼接到一起，查看不同级别就返回不同级别的图片。

这里模拟切片，第0级是一张图片，经度从-180到180，纬度从-90到90。流行方式是以左上角（-180,90）为坐标原点（origin），计算每个片的最小经纬度，最大经纬度（xmin，ymin，xmax，ymax）。

如计算第1级的四张瓦片中的右上角，直观答案是经度0到180，纬度0到90。这个数据从公式中如何计算得来呢？如上图中x，y的计算方式：最小经度xmin（0）= 原点的经度（-180）+偏移量（偏移的瓦片数量（1）*第一级瓦片宽度（以角度为单位，也就是180）），这样就可以算出每个瓦片的最小最大经纬度。

上图中用的resolution*256是什么意思，这是为了理解常说的切片方案中的参数，其中256是一张瓦片的宽度（256px像素），经过反算得出来的参数结果与天地图常用的切片方案一致。0级的resolution是0.703125，每个层级的resolution是上一层级的½。

lodinfos: [
          { 'level': 0, 'resolution': 0.703125, 'scale': 295497598.570834 },
          { 'level': 1, 'resolution': 0.3515625, 'scale': 147748799.285417 },
          { 'level': 2, 'resolution': 0.17578125, 'scale': 73874399.6427087 },
          { 'level': 3, 'resolution': 0.087890625, 'scale': 36937199.8213544 },
          { 'level': 4, 'resolution': 0.0439453125, 'scale': 18468599.9106772 },
          { 'level': 5, 'resolution': 0.02197265625, 'scale': 9234299.95533859 },
          { 'level': 6, 'resolution': 0.010986328125, 'scale': 4617149.97766929 },...
        ]

获得了这些参数信息，再输入/z/y/x，我们就获得了层级/行/列，通过公式便可计算出图片四至。这里有第一个坑，每个正方形瓦片的经度和纬度的边长是不一样的，宽（经度）是长（纬度）的两倍，可以结合第0级正方形理解，经度总共是360度，纬度总共是180度。

var lt_x = -180 + col * lodinfos[level].resolution * tilecols;
var lt_y = 90 - row * lodinfos[level].resolution * tilerows / 2;
var rb_x = -180  + (col + 1) * lodinfos[level].resolution * tilecols;
var rb_y = 90 - (row + 1) * lodinfos[level].resolution * tilerows / 2;

貌似大功告成，写了个本地服务测试两个切片结果，拉出来检验一下：用ArcGIS的切片服务看看

https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/2/1/2

两个图片的最小经度和最大经度值是一样的，但纬度差异很大。一开始我以为是公式错误，或者切片方案的问题，实际上这里是投影的问题。从一开始我的计算方案就是以地理坐标的经纬度计算的，瓦片切分下来就是正轴墨卡托投影，而ArcGIS server使用的是横轴墨卡托投影（关于投影可以看看这个一文了解地图投影 - 知乎），所以两者并不对等，当细节级别更高时两者的结果就相差十万八千里了。

二、第二次实践

思路

既然知道了原因，那么我们可以直接计算投影坐标下的四至，而不是经纬度四至。需要先获取0级的四至，即全球基准四至，然后每一级都可以通过上一级二分来计算。这里的基准四至可以从切片方案中获取也可以从ArcGIS服务中获取

      let deltX = (xmax -xmin)/Math.pow(2,level); //level级别的瓦片宽度 0级一张图，1级两张图
      let deltY = (ymax -ymin)/Math.pow(2,level); //level级别的瓦片高度 0级一张图，1级两张图
      let lt_x = xmin + col * deltX ;
      let lt_y = ymax - row * deltY ;
      let rb_x = xmin + (col + 1) * deltX ;
      let rb_y = ymax - (row + 1) * deltY ;

再拉出来看看结果

怎么回事？？？为什么还是不一样？？？

原来是取错数值了，我取原图的四至的时候使用了initial extent计算的，这个范围是超出了图本身的范围full extent，所以在计算的时候有偏差（initial不是正方形也有影响）。

更新使用fullextent参数试试12层级的图片，大体上接近了，还有少许经度差异，可以再深入研究研究。

三、最终实践

正解

第二次实践是自作聪明了。因为切片方案定义的就是以原点为起点，加偏移量去计算值。使用二分法必然不符合他原始的切片（尽管上面看起来只有微小差距，到了13级也是天差万别...）规则，纯属搞笑**行为。老老实实使用原点+切片方案计算偏移：

      const resolution = 156543.03392800014; //0级resolution（一张图）
      const originX = -2.0037508342787E7;
      const originY = 2.0037508342787E7;
      const delta = resolution*256/Math.pow(2,level);
      let lt_x = originX + col * delta ;
      let lt_y = originY - row * delta ;
      let rb_x = originX + (col + 1) * delta ;
      let rb_y = originY - (row + 1) * delta ;

检验：