MapBox 加载本地离线地形

一、参考文章：

（1）https://beginor.github.io/2021/01/11/publish-dem-data-with-mapbox-terrain-rgb.html

（2）https://blog.csdn.net/semian7633/article/details/120027672

二、相关概念：

        Mapbox Terrain-RGB 包含以 PNG 栅格切片编码的全球数字高程数据， 这些颜色值可以解码为以米为单位的原始高度。 您可以将 Terrain-RGB 数据用于各种视觉和分析应用程序， 从设计地形坡度和山体阴影样式到生成用于视频游戏的 3D 地形网格。

① 有关Terrain-RGB的一些细节：

水平精度为 5 米 对于 256x256 分辨率的切片， 缩放级别到 15 级， 512x512 分辨率的切片， 缩放级别到 14 级；垂直精度为 0.1 米 数据以 0.1 米的高度增量进行映射；坐标系为 WGS84 Web 墨卡托 wkid 为 EPSG:3389 的 Web 墨卡托坐标系， 是 WebGIS 的事实标准；

② Terrain-RGB 数据解码：

Terrain-RGB 使用每个颜色通道以 256 进制数的来表示高度，从而允许有 16,777,216 个唯一值。 收到图块后， 将需要获取各个像素的红色（R），绿色（G）和蓝色（B）值。 您可以使用浏览器中的画布层（示例）或使用诸如 get-pixel 之类的工具来执行此操作。

使用下面的公式可以将像素值解码得到高度值， 高度值以米为单位：

height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

三、操作步骤：

1、首先需要安装三个Python程序用于对tif文件进行格式转换（tif文件默认是单通道的灰度图，mapbox的Terrain-RGB仅支持rgb三通道的数据，是以每个颜色通道的 256 进制数来表示高度，因此需要进行通道转换操作）：

（1）gdal 开源的空间数据处理程序：GDAL — GDAL documentation

（2）rasterio： MapBox 在 gdal 基础上开发的栅格工具：

rasterioRasterio reads and writes geospatial raster datasets - GitHub - rasterio/rasterio: Rasterio reads and writes geospatial raster datasetshttps://github.com/mapbox/rasterio/（3）rio-rgbify： MapBox 发布的将 dem 栅格编码为 rgb 栅格的 rasterio 插件：  

rio-rgbifyEncoded arbitrary bit depth rasters in psuedo base-256 - GitHub - mapbox/rio-rgbify: Encoded arbitrary bit depth rasters in psuedo base-256https://github.com/mapbox/rio-rgbify/

注： 一定要按顺序安装，并且前两个python程序一定要安装指定版本，否则安装第三个程序时会报错。

2、GeoTiff 的坐标系必须是 WGS84 Web 墨卡托 (EPSG:3857)，可以用 rasterio 提供的命令行工具来获取tif文件信息：

命令：rio info --indent 2 3857_gd_dem.tif

3、如果不是3857坐标系，则需要进行坐标转换；此外还需要对”nodate”行的负值进行清除，因为Terrain-RGB 无法表示负值，使用 gdal 提供的 gdalwarp 命令行工具进行坐标转换的同时清除负值：

命令：gdalwarp -t_srs EPSG:3857 -dstnodata None -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_NEEDED 3857_gd_dem.tif 3857_gd_dem_n.tif

4、此时会生成一个新的tif文件，再次使用rasterio 提供的命令行工具来查看新的tif文件信息，确保上一步执行的命令生效：

5、接下来就是使用 rgbify工具 将灰度数据转换成 rgb 数据，计算高度的公式是：

height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

因此设置 ribify 的参数 base value 的参数为 -10000 ， interval 为 0.1 ， 继续输入下面的命令：

命令：rio rgbify -b -10000 -i 0.1 3857_gd_dem_n.tif 3857_gd_dem_n_rgb.tif

如果提示rgbify工具未安装的话，需要先执行" pip install rio-rgbify "命令进行安装。

完成之后， 得到的 3857_gd_dem_n_rgb.tif 看起来是这样子的：（Terrain-RGB 用 3 个 byte 通过 rgb 三通道来表示高程， 因此会比原来的灰度 tiff 要小很多）

6、把新生成的rgb三通道的tif文件进行切片：（这里牵扯到一个切片算法的概念，大多数切片工具包括geoserver自带的切片工具默认都是谷歌的xyz切片算法，但是该算法生成的切片mapbox是不识别的，mapbox支持的是osm tiles的切片算法），目前是发现使用Global Mapper工具可以使用该算法进行切片。

 （1）把新生成的tif文件用Global Mapper工具打开，选择“文件”>“输出”>“输出web格式”,在这里选择OSM Tiles:

（2）选择层级，层级越高，生成的地形精细程度越高，当然文件也就越大；取消勾选底部的“高级：总是创建8位调色板的PNG文件”（因为我们的地形默认是16或32位的，强行转换成8位的话会丢失掉很多细节），单击“确定”选择输出位置，等待切片完成：

（3）生成好的切片文件格式时这样的，第一层代表Z，Z表示缩放层级，Z=zoom；第二层代表X，第三层代表Y，XY的原点在左上角，X从左向右，对应坐标系的经度，Y从上向下，对应坐标系的纬度：

7、把切片文件目录放到本地或任意的 http 服务器就可以测试调用了（必须是mapbox-gl-js 和css 2.0以上版本才能支持加载地形服务）。

8、成功调用~

9、现存问题：

（1）精度较低的地形处理完成后可能会出现切片连接处有毛刺的现象。最好是精度10m或更高的地形数据，目前测试30m和90m精度的地形都有可能会出现上述问题。