背景

Webgis的加载服务端发布的瓦块地图数据多种格式，其中发布的数据通常可以按照服务协议（WMTS、WFS、TMS)，按照投影坐标系有经纬度投影（椭球，即WGS84，EPSG:4326坐标系)、球面墨卡托投影（WebMercator,Web墨卡托、圆球、正球,EPSG:3857），按照数据格式(mvt、geojson、kml、png等等)。

当前现有的瓦片地图下载工具有很多，免费的有locaspace viewer ,收费的有水经注等。其中实现原理并不复杂，网上也有一些开源代码的实现,我在Github上找到过相关python实现，不过没有进行验证。

本文主要是尝试给定地理矩形范围内，通过计算包围盒对应点的瓦块坐标，循环拼取URI地址，调用libcurl相关接口，将对应的图片按照z/x/y形式保存到文件夹中，最后使用Cesium加载该数据进行查看验证。代码量不大，源码没有上传。

相关原理

谷歌卫星影像地图数据是一系列组织起来图片，这些图片上像素点的坐标球面墨卡托投影坐标系，这个坐标系描述了地球椭球经纬度坐标(WGS84)简化成圆球，投射（映射、投影）到二维平面对应的坐标。具体描述可以参考 https://epsg.io/3857 ,可以看到该平面坐标系原点对应经纬度为（0,0）,东向为+X方向，北向为+Y方向，包含的地理范围是(-180.0 -85.06),(180.0 85.06),该坐标系对应的范围是（-20026376.39 -20048966.10）、（20026376.39 20048966.10)。

网络请求获取这些图片并保存下来，通常对应有xyz描述，其中z描述的是地图的切片层级，xy是对应的瓦块坐标（切片、瓦块、瓦片的索引），谷歌的这套坐标是TMS切分算法的一种变种，即在墨卡托坐标系的左上角起算，向右（东）为x方向，向下（南）为y方向，与TMS(Tile Map Service)的y轴方向相反。

计算单个经纬度点坐标对应的瓦块坐标，主要思路是如下：

  1. 经纬度坐标转换成WebMercator坐标系
  2. 计算WebMercator坐标系在指定地图层级（分辨率）下的像素坐标
  3. 计算该像素坐标针对当前瓦块大小(通常 256 * 256)对应的瓦块坐标

这方面的更多信息可以查看 https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/

其中算法部分可以参考的代码有https://github.com/datalyze-solutions/globalmaptiles/blob/master/globalmaptiles.js

我把这部分代码简化整理成脑图文件，图片如下：

downloadGoogleImage

实现

windows环境下,使用Cmake和vcpkg配置了libcurl以后就可以了(或者其他http请求技术)。

其中谷歌地图影像地址模板为 http://mt{s}.google.cn/vt/lyrs=s&hl=zh-CN&x={x}&y={y}&z={z}&s=Gali) ，其中s 对应的域名可选0,1,2,3 ,注意国内由于政策原因，无法直接访问谷歌地图服务器，需要科学上网，其中libcurl部分代码配置参考：

	CURL* _handle = curl_easy_init();
	curl_easy_setopt(_handle, CURLOPT_PROXY, "127.0.0.1:1080");
	struct curl_slist *list;
	list = curl_slist_append(NULL, "Accept-Language: zh-CN");
	list = curl_slist_append(list, "Accept:*/*");
	list = curl_slist_append(list, "Cache-Control:no-cache");
	list = curl_slist_append(list, "Accept-Encoding:gzip,deflate");

	curl_easy_setopt(_handle, CURLOPT_HTTPHEADER, list);

在代码实现中，容易出错的部分是瓦块坐标的遍理，即在瓦块坐标(maxx,maxy),(minx,miny)范围内的循环，原理上不理解很可能出错。

windows环境下，在libcurl请求回调写入二进制文件，如果使用fopen_s函数,写入模式上应当指明二进制方式写入。

std::filesystem::path imgSavePath = zxFolder/(std::to_string(y)+".jpg");
FILE *fp;
fopen_s(&fp, imgSavePath.string().c_str(), "wb+");

由于没有多线程，下载速度并不快，过程截图如下：

image-20211106204018312

验证

Cesium引擎支持经纬度投影和 WebMercator投影坐标系。由于谷歌影像数据是Web Mercator投影，可

viewer.imageryLayers.removeAll();

// TileCoordinatesImageryProvider构造函数默认使用经纬度投影
viewer.imageryLayers.addImageryProvider(
  new Cesium.TileCoordinatesImageryProvider({
    tilingScheme: new Cesium.WebMercatorTilingScheme(),
  })
);

viewer.imageryLayers.addImageryProvider(
  new Cesium.UrlTemplateImageryProvider({
    url: "http://127.0.0.1:8080/googleEarthImages/{z}/{x}/{y}.jpg",
    rectangle: Cesium.Rectangle.fromDegrees(
      86.011963,
      38.822591,
      116.63268,
      40.05064
    ),
    maximumLevel: 10,
    minimumLevel: 0,
    enablePickFeatures: false,
  })
);

image-20211106203406416

image-20211106203502286