【PythonRS】基于矢量范围批量下载遥感瓦片高清数据（天地图、高德、谷歌等）

        这个是之前写的代码了，正好今天有空所以就和大家分享一下。我们在处理项目时，有时候需要高清底图作为辅助数据源去对比数据，所以可能会需要卫星数据。所以今天就和大家分享一下如何使用Python基于矢量范围批量下载高清遥感瓦片数据。

1 读取矢量边界

        这里我们使用osgeo中的osr、ogr库去读取矢量的地理范围。之前也分享过，感兴趣的可以去Python&GIS专栏里面看一看。注意我这里只算了四至范围而不是整个矢量的边界范围，因为边界范围计算行列数不太好计算。正因如此我们下载的瓦片数据也不具备坐标系，所以也裁剪不了。

def Open_Vector(path_shp):""":param path_shp: 输入84坐标矢量:return: 返回四至范围"""ds = ogr.Open(path_shp, True)# True表示以读写方式打开layer = ds.GetLayer(0)# 获取图层feature = layer.GetFeature(0)geom = feature.GetGeometryRef()# 获取该要素的地理空间范围left, right, down, up = geom.GetEnvelope()# 获取图层的地理范围return left, right, down, up

2 通过经纬度计算航带数

        这里没什么好说的，就是基础的公式，直接计算行列数即可。这个函数在整个函数作为辅助函数，主程序会自己调用它。

def calculation_tile(lat, lon, zoom):""":param lat: 84坐标纬度:param lon: 84坐标经度:param zoom: 缩放级别:return: 瓦片的行列号"""# 将经纬度从WGS84坐标系转换为GCJ02坐标系# lon_deg,lat_deg = WGS84_To_GCJ02(lon_deg,lat_deg)# 根据缩放级别计算格网数量n = 2.0 ** zoom# 将纬度从度转换为弧度lat_radio = math.radians(lat)# 计算瓦片中的x坐标tile_x = int((lon + 180.0) / 360.0 * n)# 计算瓦片中的y坐标tile_y = int((1.0 - math.log(math.tan(lat_radio) + (1 / math.cos(lat_radio))) / math.pi) / 2.0 * n)# 返回计算得到的瓦片坐标（行和列）return tile_x, tile_y 

3 获取瓦片下载链接

        这里使用了基础的反爬虫方法，随机调用请求头，也可以自己添加或减少请求头。即使如此有时还会爬取失败，感兴趣的可以自己改进一下反爬方法，如果不会就多运行几次。

def Get_image(url, x, y):agents = ['Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.1 (KHTML, like Gecko) Chrome/13.0.782.24 Safari/535.1','Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/534.27 (KHTML, like Gecko) Chrome/12.0.712.0 Safari/534.27','Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/60.0.3112.101 Safari''/537.36','Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/532.5 (KHTML, like Gecko) Chrome/4.0.249.0 Safari''/532.5','Mozilla/5.0 (Windows; U; Windows NT 5.2; en-US) AppleWebKit/532.9 (KHTML, like Gecko) Chrome/5.0.310.0 Safari''/532.9','Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US) AppleWebKit/534.7 (KHTML, like Gecko) Chrome/7.0.514.0 Safari''/534.7','Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/9.0.601.0 ''Safari/534.14','Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/10.0.601.0 ''Safari/534.14','Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.20 (KHTML, like Gecko) Chrome/11.0.672.2 ''Safari/534.20',]try:# 打印下载成功的消息，显示瓦片的位置和下载状态print("瓦片" + str(x) + '_' + str(y) + '下载成功')# 创建一个请求对象，使用指定的URLrequests = urllib.request.Request(url)# 为请求添加一个随机的User-Agent头，以模拟不同的浏览器或客户端requests.add_header('User-Agent', random.choice(agents))  # 换用随机请求头# 使用指定的请求打开URL，并设置超时时间为60秒image = urllib.request.urlopen(requests, timeout=60)# 读取返回的图像数据image_io = image.read()# 使用BytesIO将图像数据转换为可处理的字节流对象image_bytes = io.BytesIO(image_io)# 使用PIL库打开图像image = Image.open(image_bytes)# 将图像从RGB格式转换为BGR格式（OpenCV需要的格式）image = cv2.cvtColor(np.asarray(image), cv2.COLOR_RGB2BGR)except EOFError:# 如果发生EOFError（可能由于网络问题、超时等），打印下载失败的消息并尝试重试print("瓦片" + str(x) + '_' + str(y) + '下载失败,正在重试......')Get_image(url, x, y)  # 递归调用Get_image函数进行重试# 返回处理后的图像数据return image

4 主程序

        这里就不过多解释了，我的代码注释非常完善，如果有什么不懂的，直接留言即可。整体就是筛选出最大的缩放级别，再调用行列数函数计算出瓦片的编号，再调运瓦片的下载链接，最后拼接起来即可。

# -*- coding: utf-8 -*-
"""
@Time ： 2023/4/9 14:37
@Auth ： RS迷途小书童
@File ：Vector Download Remote Sensing Tile Data.py
@IDE ：PyCharm
@Purpose：根据矢量范围下载三方地图瓦片
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import io
import cv2
import math
import random
import numpy as np
from osgeo import ogr
import urllib.request
from PIL import Imagedef Write_image(lat1, lon1, lat2, lon2, out_path):""":param lat1: 左上角纬度:param lon1: 左上角经度:param lat2: 右下角纬度:param lon2: 右下角经度:return: 返回瓦片影像"""zooms = list()# 创建一个空列表zooms，用于存储所有的缩放级别for i in range(1, 19):# 循环缩放级别col = calculation_tile(lat1, lon1, i)# 将经纬度转换为对应的地图瓦片编号，结果存储在col中row = calculation_tile(lat2, lon2, i)if col[0] - row[0] == 0 or col[1] - row[1] == 0:continueelse:zooms.append(i)# 如果差值不为0，将当前的缩放级别i添加到zooms列表中zoom = zooms[-1]# 获取zooms列表中的最后一个元素，即最大的缩放级别，并存储在zoom变量中left_up = calculation_tile(lat1, lon1, zoom)# 使用最大的缩放级别和第一个经纬度范围，调用函数获取左上角的地图瓦片编号，存储在left_up中right_down = calculation_tile(lat2, lon2, zoom)# 使用最大的缩放级别和第二个经纬度范围，调用函数获取右下角的地图瓦片编号，存储在right_down中images_columns = list()# 创建一个空列表images_columns，用于存储所有的地图瓦片图像列print("当前瓦片行数：", right_down[0]-left_up[0])print("当前瓦片列数：", right_down[1] - left_up[1])print("--------------------------------------数据获取--------------------------------------")for x in range(left_up[0], right_down[0]):# 循环行images_rows = list()# 创建一个空列表images_rows，用于存储所有的地图瓦片图像行for y in range(left_up[1], right_down[1]):# 循环列tile_url = 'http://t4.tianditu.com/DataServer?T=img_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom) + \'&tk=45c78b2bc2ecfa2b35a3e4e454ada5ce'image = Get_image(tile_url, x, y)cv2.imwrite(out_path + "/%s.jpg" % (str(x)+"_"+str(y)), image)images_rows.append(image)# 将获取到的瓦片图像添加到images_rows列表中，用于后续的图像合成img_column_new = np.vstack(images_rows)# 使用NumPy的v stack函数，将images_rows列表中的所有图像竖直堆叠起来，形成一个新的图像列images_columns.append(img_column_new)# 将这个新的图像列添加到images_columns列表中，用于后续的图像合成print("正在拼接瓦片数据......")result = np.hstack(images_columns)# 使用NumPy的h stack函数，将images_columns列表中的所有图像水平堆叠起来，形成一个最终的大图像print("正在保存瓦片数据......")cv2.imwrite(out_path + "/result.jpg", result)return result

5 总结

下列有一些三方底图链接，也可以查看【Python&GIS】第三方地图服务汇总。

""" 
tile_url = 'http://www.google.cn/maps/vt/pb=!1m4!1m3!1i'+str(zoom)+'!2i'+str(x)+'!3i'+str(y)+'!2m3!1e0!2sm!3i345013117!3m8!2szh-CN!3scn!5e1105!12m4!1e68!2m2!1sset!2sRoadmap!4e0'
# Google地图瓦片
tile_url = 'http://mt3.google.cn/vt/lyrs=s@110&hl=zh-CN&gl=cn&src=app&x='+str(x)+'&y='+str(y)+'&z='+str(zoom)+'&s=G'
# Google影像瓦片
tile_url = 'http://t4.tianditu.com/DataServer?T=img_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom)+'&tk=45c78b2bc2ecfa2b35a3e4e454ada5ce'
# 天地图卫星数据，vec_w电子地图（2000坐标系）
"http://wprd01.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scl=1&style=6&x=" + str(x) + "&y=" + str(y) + "&z=" + str(zoom) + "&ltype=3"
# 高德底图，偏移（火星坐标系）
"""

        这里输入的矢量需要是WGS84坐标系的经纬度，不能是投影坐标系哦。此外如果使用高德、百度等底图可能会有一定的偏移，因为我国需要加密成火星坐标系，但是还是可以用的，略微有偏移而已作为对比图够用了。天地图就无所谓，它的坐标是准的。