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

这个是之前写的代码了，正好今天有空所以就和大家分享一下。我们在处理项目时，有时候需要高清底图作为辅助数据源去对比数据，所以可能会需要卫星数据。所以今天就和大家分享一下如何使用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) + '下载成功')
        # 创建一个请求对象，使用指定的URL
        requests = 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 Image


def 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:
            continue
        else:
            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 总结

""" 
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"
# 高德底图，偏移（火星坐标系）
"""