gdal2tiles工具

定义解释

TMS(瓦片地图服务):tiled map service

离线安装

1.安装包准备

1.下载Python Zip免安装版本(笔记下载3.8)

2.下载gdal3.4.3 whl版本（笔者前期使用GDAL3.3.2）下载链接列表http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

2.安装命令

安装whl命令（pip install *.whl安装离线库，这种安装一般是二进制包）

pip install *.whl

卸载安装包命令

pip uninstall 包名。这里包名就是gdal

pip uninstall gdal

因为是没有设置python环境变量的。在使用命令是需要指定前面python.exe路径。可能命令安装命令如下：

python.exe完整路径 -m pip install *.whl文件完整路径

卸载命令安装命令如下：

python.exe完整路径  pip uninstall 包名称（这里是gdal）

安装参考开源影像tif切图工具gdal2tiles部署以及切图

问题一：

使用gdal2tiles.py时,报错

Please convert this file to 8-bit and run gdal2tiles on the result

分析主要原因是因为该tif的深度信息不是8位。需要将tif的深度信息转换为8位。

方法1.gdal_translate.exe工具

使用

gdal_translate.exe     gdal_tanslate -of VRT -ot Byte -scale 输入的tif路径 temp.vrt

或默认是输出tif，可以使用如下语句

 gdal_tanslate -ot Byte -scale 输入的tif路径 d:\out\1.tif

再使用gdal2tiles.py 调用执行

方法2.ArcMap进行转换

主要使用ArcMap 栅格|复制栅格工具。参考arcgis如何将16bit栅格数据转换为8bit栅格数据

问题二：

如果过程中报错ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db 解决办法设置环境变量 PROJ_LIB="proj.db所在目录"

解决方案一：配置Proje_LIB环境变量

解决方案二：修改代码

Scripts/gdal2tile.py文件代码中添加如下代码（更灵活，推荐此方法）

import os
BASE_DIR = os.path.dirname(os.path.dirname(__file__))
path = BASE_DIR + '\Lib\site-packages\osgeo\data\proj'
os.environ['PROJ_LIB'] = path

问题三：

源于同事问CGS2000坐标系的tif经过gdal2tiles可以在Cesium三维球中正常加载（他说Cesium三维球使用的4326坐标），问gdal2tiles是否已经做了坐标转换?顿时一年茫然。因为gdal2tiles.py代码没看过。只是将gdal python环境搭建起来，可以运行，没有关注里面的源代码。

1.猜测

gdal2tiles.py [-p profile] [-r resampling] [-s srs] [-z zoom]
[-e] [-a nodata] [-v] [-q] [-h] [-k] [-n] [-u url]
[-w webviewer] [-t title] [-c copyright]
[--processes=NB_PROCESSES] [--xyz]
--tilesize=PIXELS
[-g googlekey] [-b bingkey] input_file [output_dir] [COMMON_OPTIONS]

gdal2tiles.py -p profile Tile切割轮廓（墨卡托，大地测量，栅格）-默认“墨卡托”（谷歌地图兼容）.默认为墨卡托（墨卡托投影为EPSG:3857）.而Cesium三维球使用EPSG:4326坐标系。（1.可能Cesium前端代码中使用投影坐标系转4326地理坐标系代码。2.可能Cesium加载墨卡托投影为EPSG:3857直接动态转换到EPSG:4326坐标系)。关于源码的理解可以参考How it works(14) GDAL2Tiles源码阅读

2.代码阅读后

1.gdal2tiles中确实有坐标转换的逻辑，且目标坐标系为EPSG:3857(墨卡托投影)和EPSG:4326（地理投影）

2.切片可以使用--processses参数传入进程多进程的切割。

if nb_processes == 1:
        single_threaded_tiling(input_file, output_folder, options)
    else:
        multi_threaded_tiling(input_file, output_folder, options)

3.create_base_tile切割顶级瓦片，create_overview_tiles切割下层瓦片

多线程方式切图报错的问题

指定 --processes参数时报错。大致是python版本与gdal版本不匹配导致的。参考的文章指出python3.7没有问题，我用的python3.8就出现如下问题。具体怎样解决需要实际的验证和尝试。

 is not going to be frozen to produce an executable.    exec(code, run_globals)

  File "C:\Users\Administrator\Desktop\03buid\Python38\Scripts\gdal2tiles.py", line 11, in <module>
    sys.exit(main(sys.argv))
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\site-packages\osgeo_utils\gdal2tiles.py", line 3260, in main
    multi_threaded_tiling(input_file, output_folder, options)
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\site-packages\osgeo_utils\gdal2tiles.py", line 3209, in multi_threaded_tiling
    pool = Pool(processes=nb_processes)
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\context.py", line 119, in Pool
    return Pool(processes, initializer, initargs, maxtasksperchild,
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\pool.py", line 212, in __init__
    self._repopulate_pool()
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\pool.py", line 303, in _repopulate_pool
    return self._repopulate_pool_static(self._ctx, self.Process,
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\pool.py", line 326, in _repopulate_pool_static
    w.start()
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\process.py", line 121, in start
    self._popen = self._Popen(self)
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\context.py", line 327, in _Popen
    return Popen(process_obj)
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\popen_spawn_win32.py", line 45, in __init__
    prep_data = spawn.get_preparation_data(process_obj._name)
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\spawn.py", line 154, in get_preparation_data
    _check_not_importing_main()
  File "C:\Users\Administrator\Desktop\03buid\Python38\lib\multiprocessing\spawn.py", line 134, in _check_not_importing_main
    raise RuntimeError('''
RuntimeError:
        An attempt has been made to start a new process before the
        current process has finished its bootstrapping phase.

        This probably means that you are not using fork to start your
        child processes and you have forgotten to use the proper idiom
        in the main module:

            if __name__ == '__main__':
                freeze_support()
                ...

        The "freeze_support()" line can be omitted if the program

Error in gdal2tiles with --processes

Gdal2tiles suffers from Python multiprocessing bug #42949 #4951

从以上链接信息好像gdal3.4后端版本没有此问题（自己当前用Python为3.8而GDAL3.3.2版本），需要亲手验证.

经过验证GDAL版本换成GDAL-3.4.3版本后，此问题解决

需要补充理解的知识

对于 Web Map 开发人员来说，最熟悉的应该是EPSG:4326 (WGS84) and EPSG:3857(Pseudo-Mercator)，

EPSG:4326 (WGS84)

前面说了 WGS84 是目前最流行的地理坐标系统。在国际上，每个坐标系统都会被分配一个 EPSG 代码，EPSG:4326 就是 WGS84 的代码。GPS是基于WGS84的，所以通常我们得到的坐标数据都是WGS84的。一般我们在存储数据时，仍然按WGS84存储。

EPSG:3857 (Pseudo-Mercator)

伪墨卡托投影，也被称为球体墨卡托，Web Mercator。它是基于墨卡托投影的，把 WGS84坐标系投影到正方形。我们前面已经知道 WGS84 是基于椭球体的，但是伪墨卡托投影把坐标投影到球体上，这导致两极的失真变大，但是却更容易计算。这也许是为什么被称为”伪“墨卡托吧。另外，伪墨卡托投影还切掉了南北85.051129°纬度以上的地区，以保证整个投影是正方形的。因为墨卡托投影等正形性的特点，在不同层级的图层上物体的形状保持不变，一个正方形可以不断被划分为更多更小的正方形以显示更清晰的细节。很明显，伪墨卡托坐标系是非常显示数据，但是不适合存储数据的，通常我们使用WGS84 存储数据，使用伪墨卡托显示数据。

Web Mercator 最早是由 Google 提出的，当前已经成为 Web Map 的事实标准。但是也许是由于上面”伪“的原因，最初 Web Mercator 被拒绝分配EPSG 代码。于是大家普遍使用 EPSG:900913（Google的数字变形）的非官方代码来代表它。直到2008年，才被分配了EPSG:3785的代码，但在同一年没多久，又被弃用，重新分配了 EPSG:3857 的正式代码，使用至今。