测试使用 Python GDAL 下载 Mapbox 瓦片地图及拼接

本教程将展示如何以编程方式从网络地图（通常称为瓦片地图）瓦片服务器下载地图图像，对其进行地理参考（设置坐标系）并将其保存为GeoTIFF。

Code

import lib

#!/usr/bin/env python
# coding: utf-8
import urllib.request
import os
from math import log, tan, radians, cos, pi,floor,degrees,atan, sinh
from osgeo import gdal
import glob

设置环境变量

temp_dir = os.path.join(os.getcwd(), 'temp')
output_dir = os.path.join(os.getcwd(), 'output')
os.environ['GDAL_DATA'] = '/home/wblong/Anaconda3/share/gdal'
os.environ['PROJ_LIB'] = '/home/wblong/Anaconda3/share/proj'

瓦片XYZ与经纬度换算

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))

def mercatorToLat(mercatorY):
    return(degrees(atan(sinh(mercatorY))))

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 x_to_lon_edges(x, z):
    tile_count = pow(2, z)
    unit = 360 / tile_count
    lon1 = -180 + x * unit
    lon2 = lon1 + unit
    return (lon1, lon2)

def y_to_lat_edges(y, z):
    tile_count = pow(2, z)
    unit = 1 / tile_count
    relative_y1 = y * unit
    relative_y2 = relative_y1 + unit
    lat1 = mercatorToLat(pi * (1 - 2 * relative_y1))
    lat2 = mercatorToLat(pi * (1 - 2 * relative_y2))
    return(lat1, lat2)

# 根据瓦片{X,Y,Z}计算经纬度范围

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]

瓦片下载、设置空间参考及合并

# 对瓦片png添加空间参考EPSG：4326
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)
#地图瓦片合并拼接
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)
# 瓦片下载
def download_tile(x, y, z, tile_server):
    url = tile_server.replace(
        "{x}", str(x)).replace(
        "{y}", str(y)).replace(
        "{z}", str(z))
    path = f'{temp_dir}/{x}_{y}_{z}.png'
    urllib.request.urlretrieve(url, path)
    return(path)
    

测试下载Mapbox瓦片

tile_server = "https://api.mapbox.com/v4/mapbox.satellite/{z}/{x}/{y}.png?access_token=xxxxxxxxxxxxxxxxxxxxxxxxxx"

zoom = 16

lon_min = 21.49147
lon_max = 21.5
lat_min = 65.31016
lat_max = 65.31688

x_min, x_max, y_min, y_max = bbox_to_xyz(
    lon_min, lon_max, lat_min, lat_max, zoom)

for x in range(x_min, x_max + 1):
    for y in range(y_min, y_max + 1):
        print(f"{x},{y}")
        png_path = download_tile(x, y, zoom, tile_server)
        georeference_raster_tile(x, y, zoom, png_path)

print("Download complete")
print("Merging tiles")
merge_tiles(temp_dir + '/*.tif', output_dir + '/merged.tif')
print("Merge complete")

合并后

参考

1. https://dev.to/jimutt/generate-merged-geotiff-imagery-from-web-maps-xyz-tile-servers-with-python-4d13
2. https://blog.csdn.net/mrbaolong/article/details/137610742?spm=1001.2014.3001.5502
3. https://blog.csdn.net/mrbaolong/article/details/137479780?spm=1001.2014.3001.5502