Cesium - 地图下载器(python)

开发 需要卫星地图作为底图。百度、高德、谷歌、arcgis都提供在线服务，但在不能访问在线服务的地方就不适用了。

网上有不少地图下载器，BigeMap、水经注、太乐、91卫图、Google Maps Downloader等等，但是它们都不是完全免费的（免费版有水印，或限制功能），想要授权就得买或者帮着推广。

于是自己编程实现。

一、基础知识

1、地图瓦片

各服务商既然提供在线服务，就必然有取得瓦片的方式。

以谷歌(可能为了和国内服务商提供的道路图层叠加，谷歌的卫星影像在国内是偏移的)为例：
http://mt1.google.cn/vt/lyrs=s&hl=zh-CN&x=3&y=1&z=2&s=Gali

瓦片大小为：256x256。

1. lyrs 表示的是图层类型，即瓦片类型，具体含义如下：

    m：路线图
    t：地形图
    p：带标签的地形图
    s：卫星图
    y：带标签的卫星图
    h：标签层（路名、地名等）
   

2. x , y 是瓦片坐标系的坐标值，z代表缩放级别
   x 瓦片的横向索引，起始位置为最左边，数值为0，向右+1递增。
   y 瓦片的纵向索引，起始位置为最上面，数值为0，向下+1递增。
   z 地图的级别，以Google为例，最上一级为0，向下依次递增。

3. 地图切图方式
   一幅地图由4^n个256的正方形组成，n为级别。
   例如：第0级为4^0个，即世界地图由一个256图片表示。
   第1级世界地图应由4^1 = 4个256图片组成，即世界地图等分成4块256图片。
   往下每一级依此类推……

2、Web墨卡托投影

Web墨卡托投影是 互联网地图通用的地图投影方式，将椭圆形地图投影成平面的正方形。

• Bounds(地图范围)：
  [ -20037508.3427892, -20037508.3427892, 20037508.3427892, 20037508.3427892]，单位为米，20037508.3427892表示地图周长的一半，以地图中心点做为（0，0）坐标。

• Levels：地图的级别，例如：0……22。

• Resolutions：
  分辨率数组，与级别相对应，即一个级别对应一个分辨率，分辨率表示当前级别下单个像素代表的地理长度。
  Resolutions[n] = 20037508.3427892 * 2 / 256 / (2^n)

• Center：地图显示中心点。

• Level：地图显示级别。

• viewSize：地图控件窗口的大小。

根据已知地图中心点、显示级别可以将地图显示范围计算出来:

viewBounds = [Center.x - Resolutions[l]*viewSize.width/2,
              Center.y - Resolutions[l]*viewSize.height/2,
              Center.x + Resolutions[l].viewSize.width/w, 
              Center.y + Resolutions[l].viewSize.height/h]

二、Python实现

使用python爬取地图瓦片，核心逻辑主要分为以下几步：

1. 确定下载的级别，经纬度范围（可以通过百度坐标拾取获得）
   如果是全球范围，设为[ -180, 90, 180, -90 ] 即可。

2. 计算出这个范围内瓦片的起始和终止行列号

3. 根据行列号拼接出瓦片的url地址

4. 下载保存图片

1、初始化

先导入依赖包和初始化全局参数。

import math
from math import floor, pi, log, tan, atan, exp
from threading import Thread, Lock
import urllib.request as ur
import PIL.Image as pil
import io
import random
import os

MAP_URLS = {
   
   
    "arcgis": "https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
    "google": "http://mt1.google.cn/vt/lyrs={style}&hl=zh-CN&gl=CN&src=app&x={x}&y={y}&z={z}&s=Gali"}

agents = [
    'Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/60.0.3112.101 Safari/537.36',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/532.5 (KHTML, like Gecko) Chrome/4.0.249.0 Safari/532.5',
    'Mozilla/5.0 (Windows; U; Windows NT 5.2; en-US) AppleWebKit/532.9 (KHTML, like Gecko) Chrome/5.0.310.0 Safari/532.9',
    'Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US) AppleWebKit/534.7 (KHTML, like Gecko) Chrome/7.0.514.0 Safari/534.7',
    'Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/9.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/10.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.20 (KHTML, like Gecko) Chrome/11.0.672.2 Safari/534.20',
    'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/534.27 (KHTML, like Gecko) Chrome/12.0.712.0 Safari/534.27',
    'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.1 (KHTML, like Gecko) Chrome/13.0.782.24 Safari/535.1']

2、经纬度转瓦片坐标

# 根据WGS-84的经纬度获取地图中的瓦片坐标
def wgs84_to_tile(lon, lat, zoom):
    isnum = lambda x: isinstance(x, int) or isinstance