最近想要做一下古建筑物的样本标注,想下载一下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)
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教程欢迎持续关注编程学习网。
扫码二维码 获取免费视频学习资料
- 本文固定链接: http://www.phpxs.com/post/11782/
- 转载请注明:转载必须在正文中标注并保留原文链接
- 扫码: 扫上方二维码获取免费视频资料
查 看2022高级编程视频教程免费获取