测试使用 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")
合并后
使用python_GDAL测试地图瓦片下载及拼接测试