矢量切片下载并导出SHP文件思路及代码

背景：

最近手里有矢量切片数据（z/x/y.pbf），但是找不到原始shp数据文件了，想着如何让矢量切片合成原始的shp数据，今天把步骤和相关代码整理出来。

部分互联网还有反爬虫机制，有的可以通过模拟浏览器请求获的。

如果对你有帮助，麻烦点赞+收藏，欢迎关注。

步骤：

1. 矢量切片数据准备

无论是自己有的矢量切片还是从其他渠道获得矢量切片，存储路径保持为z/x/y.pbf。有了矢量切片可以直接进行第2步。

下载互联网矢量需要主要反爬虫设置来限制获取数据。

文件存储样例

Python代码：

直接修改main()的base_url地址就可以，代码通过请求头来了模拟浏览器可以来代码下载。（设置的是全国范围extend，所以会出现瓦片不存在，下载失败的正常，是由通过QGIS产看是否正确就可。）

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
智能GIS瓦片下载器 - 处理瓦片缺失和区域覆盖问题
"""

import requests
import os
import math
import time
from urllib.parse import urlparse
import sys
import json

class SmartGISTileDownloader:
    def __init__(self, base_url, output_dir="tiles", max_retries=3, chunk_size=8192):
        """
        初始化智能GIS瓦片下载器
        
        Args:
            base_url (str): 瓦片服务的基础URL
            output_dir (str): 输出目录
            max_retries (int): 最大重试次数
            chunk_size (int): 分块大小
        """
        self.base_url = base_url
        self.output_dir = output_dir
        self.max_retries = max_retries
        self.chunk_size = chunk_size
        
        # 设置请求头，模拟真实浏览器
        self.headers = {
            'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/120.0.0.0 Safari/537.36',
            'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,image/apng,*/*;q=0.8',
            'Accept-Language': 'zh-CN,zh;q=0.9,en;q=0.8',
            'Accept-Encoding': 'gzip, deflate, br',
            'Connection': 'keep-alive',
            'Upgrade-Insecure-Requests': '1',
        }
        
        # 创建会话
        self.session = requests.Session()
        self.session.headers.update(self.headers)
        
        # 设置适配器
        adapter = requests.adapters.HTTPAdapter(
            pool_connections=10,
            pool_maxsize=10,
            max_retries=3,
            pool_block=False
        )
        self.session.mount('http://', adapter)
        self.session.mount('https://', adapter)
        
        # 统计信息
        self.stats = {
            'total_tiles': 0,
            'successful_downloads': 0,
            'failed_downloads': 0,
            'missing_tiles': 0,
            'skipped_tiles': 0
        }
        
        # 瓦片缺失记录
        self.missing_tiles_log = []
        
        # 区域覆盖分析
        self.coverage_analysis = {}
    
    def analyze_coverage(self, min_lat, max_lat, min_lon, max_lon, z_range=(0, 15)):
        """
        分析瓦片服务的覆盖情况
        """
        print("=== 瓦片服务覆盖分析 ===")
        print(f"地理范围: 纬度({min_lat}, {max_lat}), 经度({min_lon}, {max_lon})")
        print(f"缩放级别范围: {z_range[0]} - {z_range[1]}")
        print()
        
        # 测试不同缩放级别的覆盖情况
        test_levels = [0, 3, 6, 9, 12, 15] if z_range[1] >= 15 else list(range(z_range[0], z_range[1] + 1, 3))
        
        for z in test_levels:
            if z > z_range[1]:
                break
                
            print(f"分析 z={z} 级别覆盖情况...")
            
            # 计算边界瓦片
            min_x, max_x, min_y, max_y = self.calculate_tile_bounds(min_lat, max_lat, min_lon, max_lon, z)
            
            # 采样测试（避免测试所有瓦片）
            sample_size = min(20, (max_x - min_x + 1) * (max_y - min_y + 1))
            sample_coords = self.generate_sample_coordinates(min_x, max_x, min_y, max_y, sample_size)
            
            available_tiles = 0
            missing_tiles = 0
            
            for x, y in sample_coords:
                if s