Python+gdal制作一个简单的地图下载（支持高德、arcgis、google）

原创 mu 在小岛学gis的穆 2022-04-26 14:55

最近想要做一下古建筑物的样本标注，想下载一下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)




下载高德18级 3gb，大概总共40min。

确定

• 不看此公众号