PYTHON 实现爬取天地图瓦片生成GeoTIFF

一、引言

在地理信息系统（GIS）相关的开发和应用中，天地图的瓦片数据是非常有用的资源。通过 Python，我们可以方便地编写程序来爬取天地图的瓦片数据，并将其保存和合并。本文将详细介绍如何使用 Python 实现这一功能。

二、实现思路

整个程序的实现主要分为以下几个步骤：

1. 经纬度与瓦片编号的转换：将 WGS84 经纬度转换为瓦片编号，以及反向转换。
2. 瓦片下载：根据瓦片编号，构造请求 URL，下载瓦片数据。
3. 多线程下载：使用多线程技术，提高下载效率。
4. 瓦片合并：将下载的瓦片合并为 GeoTIFF 文件。

三、代码实现

import math
import os
import requests
import random
import time
import threading
from queue import Queue
import rasterio
from rasterio.merge import merge
from rasterio.crs import CRS
from rasterio.transform import Affine
import numpy as np

# 天地图API密钥
TIANDITU_KEY = "<填写密钥，从天地图官网注册后可获取自己的访问密钥>"


# 定义WGS84经纬度转瓦片编号的函数
def lonlat_to_tile(lon, lat, zoom):
    n = 2.0**zoom
    xtile = int((lon + 180.0) / 360.0 * n)
    lat_rad = math.radians(lat)
    ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
    return xtile, ytile


# 定义下载瓦片的函数
def download_tile(zoom, x, y, save_path, tk=TIANDITU_KEY):
    t = random.randint(0, 7)
    url = f"https://t{t}.tianditu.gov.cn/img_w/wmts?tk={tk}&Service=WMTS&Request=GetTile&Version=1.0.0&Style=Default&Format=tiles&serviceMode=KVP&layer=img&TileMatrixSet=w&TileMatrix={zoom}&TileRow={y}&TileCol={x}"
    headers = {
        "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/133.0.0.0 Safari/537.36",
        "accept": "text/html,application/xhtml+xml,application/xml;q=0.9,image/avif,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3;q=0.7",
        "accept-encoding": "gzip, deflate, br, zstd",
        "accept-language": "zh-CN,zh;q=0.9",
    }
    max_retries = 5
    for attempt in range(max_retries):
        try:
            response = requests.get(url, headers=headers)
            if response.status_code == 200:
                with open(save_path, "wb") as f:
                    f.write(response.content)
                print(f"Downloaded tile {zoom}-{x}-{y}")
                return
            else:
                print(
                    f"Attempt {attempt + 1} failed to download tile {zoom}-{x}-{y}, status code: {response.status_code}"
                )
        except Exception as e:
            print(f"Attempt {attempt + 1} error downloading tile {zoom}-{x}-{y}: {e}")
        if attempt < max_retries - 1:
            time.sleep(5)
    print(f"Failed to download tile {zoom}-{x}-{y} after {max_retries} attempts.")


# 线程工作函数
def worker(queue):
    while True:
        item = queue.get()
        if item is None:
            break
        zoom, x, y, save_path = item
        download_tile(zoom, x, y, save_path)
        # 随机等待一段时间，避免被识别为爬虫
        time.sleep(random.uniform(0.5, 2))
        queue.task_done()


# 定义下载指定区域瓦片的函数
def download_tiles(
    min_lon,
    min_lat,
    max_lon,
    max_lat,
    zoom_levels,
    output_dir,
    num_threads=4,
    tk=TIANDITU_KEY,
):
    tile_queue = Queue()
    threads = []

    # 创建线程
    for _ in range(num_threads):
        t = threading.Thread(target=worker, args=(tile_queue,))
        t.start()
        threads.append(t)

    for zoom in zoom_levels:
        min_x, min_y = lonlat_to_tile(min_lon, max_lat, zoom)
        max_x, max_y = lonlat_to_tile(max_lon, min_lat, zoom)

        zoom_folder = os.path.join(output_dir, f"zoom_{zoom}")
        if not os.path.exists(zoom_folder):
            os.makedirs(zoom_folder)

        for x in range(min_x, max_x + 1):
            for y in range(min_y, max_y + 1):
                save_path = os.path.join(zoom_folder, f"{x}_{y}.png")
                tile_queue.put((zoom, x, y, save_path))

    # 等待所有任务完成
    tile_queue.join()

    # 停止线程
    for _ in range(num_threads):
        tile_queue.put(None)
    for t in threads:
        t.join()


# 新增函数：瓦片坐标转经纬度（与lonlat_to_tile对应）
def tile_to_lonlat(xtile, ytile, zoom):
    n = 2.0 ** zoom
    lon_deg = (xtile / n) * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return lon_deg, lat_deg


# 修改合并函数
def merge_tiles_to_tiff(zoom, output_dir, min_lon, min_lat, max_lon, max_lat):
    zoom_folder = os.path.join(output_dir, f"zoom_{zoom}")
    tile_files = []
    for root, dirs, files in os.walk(zoom_folder):
        for file in files:
            if file.endswith(".png"):
                tile_files.append(os.path.join(root, file))

    src_files_to_mosaic = []
    for f in tile_files:
        # 从文件名提取x和y
        filename = os.path.basename(f)
        x_str, y_str = os.path.splitext(filename)[0].split('_')
        x = int(x_str)
        y = int(y_str)
        
        # 计算瓦片左上角和右下角经纬度
        ul_lon, ul_lat = tile_to_lonlat(x, y, zoom)
        lr_lon, lr_lat = tile_to_lonlat(x + 1, y + 1, zoom)
        
        # 计算分辨率（确保像素高度为负）
        res_lon = (lr_lon - ul_lon) / 256  # 假设瓦片为256x256像素
        res_lat = (lr_lat - ul_lat) / 256  # res_lat自动为负
        
        # 创建Affine变换矩阵
        transform = Affine(res_lon, 0, ul_lon,
                          0, res_lat, ul_lat)
        
        # 打开文件并手动设置变换矩阵和CRS
        with rasterio.open(f) as src:
            array = src.read()
            profile = src.profile.copy()
        
        profile.update({
            'transform': transform,
            'crs': CRS.from_epsg(4326),
            'dtype': array.dtype
        })
        
        # 使用内存文件加载数据
        memfile = rasterio.MemoryFile()
        # with rasterio.MemoryFile() as :
        with memfile.open(**profile) as dataset:
            dataset.write(array)
        src = memfile.open()
        src_files_to_mosaic.append(src)

    # 合并瓦片
    mosaic, out_trans = merge(src_files_to_mosaic)
    
    # 输出元数据设置
    out_meta = src_files_to_mosaic[0].meta.copy()
    out_meta.update({
        "driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": out_trans,
        "crs": CRS.from_epsg(4326)
    })
    
    output_file = os.path.join(output_dir, f"merged_zoom_{zoom}.tif")
    with rasterio.open(output_file, "w", **out_meta) as dest:
        dest.write(mosaic)
    
    for src in src_files_to_mosaic:
        src.close()
    
    print(f"Merged tiles for zoom level {zoom} into {output_file}")


if __name__ == "__main__":
    # 指定区域范围（WGS84经纬度）
    min_lon = 120.295
    max_lon = 120.341
    min_lat = 32.745
    max_lat = 32.7783

    # 指定层级
    zoom_levels = [18] #支持1-18级，级号越大图片越清晰，但是生成的文件也会越大

    # 指定线程数
    num_threads = 6

    # 指定输出目录
    output_dir = "e:\\outputGeoTiff"
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 下载瓦片
    download_tiles(min_lon, min_lat, max_lon, max_lat, zoom_levels, output_dir, num_threads)

    # 合并每个层级的瓦片为GeoTIFF文件
    for zoom in zoom_levels:
        merge_tiles_to_tiff(zoom, output_dir, min_lon, min_lat, max_lon, max_lat)

本代码旨在为技术交流与学习提供示例，严禁将其用于任何违法违规活动。使用者需自行承担因不当使用引发的一切法律责任。