Python教程：Python+gdal制作一个简单的地图下载！

最近想要做一下古建筑物的样本标注，想下载一下rgb影像。。某新地球和xxmap都不能下载17以上（含17）的，急死俺了，只能去github上找一下参考，自己去做一个。这边来说一下自己做的一个思路。
github地址
https://github.com/sadassimov/tiledownload


Web 地图和 XYZ 切片格式
一般在线地图使用共享的方式来呈现地图。大多数地图是通过使用多个 256x256px 光栅图块进行可视化的。通过提供它们的 x、y 和 z 坐标来加载正确的图块。其中 z 对应于当前缩放级别。
例如,可以使用遵循以下模式的 URL 获取openstreetmap：https://tile.openstreetmap.org/{z}/{x}/{y}.png
x、y 和 z 参数是整数.
osm原文：
X goes from 0 (left edge is 180 °W) to 2^zoom − 1 (right edge is 180 °E) Y goes from 0 (top edge is 85.0511 °N) to 2^zoom − 1 (bottom edge is 85.0511 °S) in a Mercator projection
For the curious, the number 85.0511 is the result of arctan(sinh(π)). By using this bound, the entire map becomes a (very large) square.
将纬度和经度转换为对应的瓦片
为了从纬度和经度推导出瓦片名称，必须执行一些操作。
1.计算当前缩放级别的瓦片总数：tile_count = 2^z
2.将 lat 和 lon 重新投影到墨卡托投影 (EPSG:3857)。x = lon
y = log(tan(lat) + sec(lat))
3.变换 x 和 y 值以适应 0 - 1 的范围，并将原点设置为左上角。x = (lon + 180) / 360
y = (1 - (y / π)) / 2
将 x 和 y 乘以图块的数量。
#用于从 WGS84 坐标转换到指定缩放级别的相应图块
def sec(x):
    return (1 / cos(x))
def latlon_to_xyz(lat, lon, z):
    tile_count = pow(2, z)
    x = (lon + 180) / 360
    y = (1 - log(tan(radians(lat)) + sec(radians(lat))) / pi) / 2
    return (tile_count * x, tile_count * y)

接下来，我们将添加一个新函数，它允许我们获取指定矩形区域/边界框的图块范围.
def bbox_to_xyz(lon_min, lon_max, lat_min, lat_max, z):
    x_min, y_max = latlon_to_xyz(lat_min, lon_min, z)
    x_max, y_min = latlon_to_xyz(lat_max, lon_max, z)
    return(floor(x_min), floor(x_max),
           floor(y_min), floor(y_max))

瓦片下载
要下载地图切片，只需发送与切片服务器提供的 URL 规范匹配的 GET 请求。（请遵守规范）
def donwload(tile_source, output_dir, bounding_box, zoom):
    lon_min, lat_min, lon_max, lat_max = bounding_box

    # Script start:
    if not os.path.exists(temp_dir):
        os.makedirs(temp_dir)

    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    x_min, x_max, y_min, y_max = bbox_to_xyz(
        lon_min, lon_max, lat_min, lat_max, zoom)
    alltiles = (x_max - x_min + 1) * (y_max - y_min + 1)
    print(f"总共：{alltiles} 瓦片")
    i = 0
    for x in range(x_min, x_max + 1):
        for y in range(y_min, y_max + 1):
            try:
                i = i + 1
                png_path = fetch_tile(x, y, zoom, tile_source)
                print(f"{x},{y} 下载进度{i}/{alltiles}")
                georeference_raster_tile(x, y, zoom, png_path)
            except OSError:
                print(f"错误, failed to get {x},{y}")
                pass

    print("下载与转换瓦片成功")
    print("整合瓦片中请稍后。。。")
    merge_tiles(temp_dir + '/*.tif', output_dir + '/out.tif')
    print(f"整合完成！输出至{output_dir}/out.tif'")
    shutil.rmtree(temp_dir)

地理配准
为了对下载的瓦片进行地理配准，我们首先需要将我们从纬度和经度执行的转换反转为瓦片名称格式。我们需要知道瓷砖四个角中每个角的坐标。然后我们将使用gdal.Translate将图像保存为带有嵌入地理位置数据的 TIFF 文件。
对于纬度转换，我们将使用相同的方法，但添加了一个额外的函数来从墨卡托投影转换回来。

def x_to_lat_edges(x, z):
    tile_count = pow(2, z)
    unit = 360 / tile_count
    lon1 = -180 + x * unit
    lon2 = lon1 + unit
    return(lon1, lon2)
    
 def tile_edges(x, y, z):
    lat1, lat2 = y_to_lat_edges(y, z)
    lon1, lon2 = x_to_lon_edges(x, z)
    return[lon1, lat1, lon2, lat2]
 
def georeference_raster_tile(x, y, z, path):
    bounds = tile_edges(x, y, z)
    filename, extension = os.path.splitext(path)
    gdal.Translate(filename + '.tif',
                   path,
                   outputSRS='EPSG:4326',
                   outputBounds=bounds)
  
      for x in range(x_min, x_max + 1):
        for y in range(y_min, y_max + 1):
            try:
                i = i + 1
                png_path = fetch_tile(x, y, zoom, tile_source)
                print(f"{x},{y} 下载进度{i}/{alltiles}")
                georeference_raster_tile(x, y, zoom, png_path)
            except OSError:
                print(f"错误, failed to get {x},{y}")
                pass

将图块合并为tif
由于我们已经包含了对 GDAL 的依赖项，因此我选择也使用 GDAL python 脚本来合并平铺图像。我们将调用gdal_merge.py。

def merge_tiles(input_pattern, output_path):
    vrt_path = temp_dir + "/tiles.vrt"
    gdal.BuildVRT(vrt_path, glob.glob(input_pattern))
    gdal.Translate(output_path, vrt_path)
通过PysimpleGUI封装界面

设置了输出目录，地图级别，下载范围，下载影像，PysimpleGUI的控件通过  key='folder' / values['folder']传输值。

def gui(): 
    layout = [
        [sg.FolderBrowse('选择输出文件夹', key='folder', target='file'), sg.Button('开始'), sg.Button('关闭')],
        [sg.Text('输出文件夹为:', font=("宋体", 10)), sg.Text('', key='file', size=(50, 1), font=("宋体", 10))],
        [sg.Combo(['高德地图', '谷歌地图', 'arcgis地图'], key='tile_source', default_value='高德地图', size=(21, 1)),
         sg.Text('下载地图级别'),
         sg.InputText(size=(20, 1), key='zoom')],
        [sg.Text('lng_min'), sg.InputText(size=(19, 1), key='lng_min'), sg.Text('lat_min'),
         sg.InputText(size=(19, 1), key='lat_min')],
        [sg.Text('lng_max'), sg.InputText(size=(19, 1), key='lng_max'), sg.Text('lat_max'),
         sg.InputText(size=(19, 1), key='lat_max')],
        [sg.Output(size=(70, 5), font=("宋体", 10))]
    ]

    window = sg.Window('在线地图下载器 -_-', layout, font=("宋体", 10), default_element_size=(50, 1), icon='./proj/earth.ico')

    while True:
        event, values = window.read()
        if event in (None, '关闭'):  # 如果用户关闭窗口或点击`关闭`
            break
        if event == '开始':
            if values['tile_source'] == '高德地图':
                tile_source = 'https://webst01.is.autonavi.com/appmaptile?style=6&x={x}&y={y}&z={z}'
            elif values['tile_source'] == '谷歌地图':
                tile_source = 'http://www.google.cn/maps/vt?lyrs=s@189&gl=cn&x={x}&y={y}&z={z}'
            elif values['tile_source'] == 'arcgis地图':
                tile_source = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}'
            output = values['folder']
            lng_min = int(values['lng_min'])
            lat_min = int(values['lat_min'])
            lng_max = int(values['lng_max'])
            lat_max = int(values['lat_max'])
            zoom = int(values['zoom'])
            donwload(tile_source, output, [lng_min, lat_min, lng_max, lat_max], zoom)

以上就是“Python教程：Python+gdal制作一个简单的地图下载！”的详细内容，想要了解更多Python教程欢迎持续关注编程学习网。