ArcGIS REST API 逆向实战：Python 脚本批量获取东莞控规 36972 条数据

东莞作为粤港澳大湾区的重要节点城市，其城市规划数据对研究者、开发商和政府机构都具有重要参考价值。控制性详细规划（简称"控规"）作为城市规划体系中的核心内容，包含了地块用途、容积率、建筑高度等关键指标。本文将详细介绍如何通过逆向工程获取东莞控规数据，并构建完整的Python自动化采集方案。

1. 技术实现原理与接口分析

1.1 ArcGIS REST API 工作机制

ArcGIS REST API 是Esri公司提供的一套基于HTTP协议的GIS服务接口标准，采用RESTful架构设计。其核心特点包括：

• 资源导向 ：每个GIS服务（如地图、要素、几何等）都有唯一的URL标识
• 无状态通信 ：每个请求包含完整参数，服务器不保存客户端状态
• 多种返回格式 ：支持JSON、HTML、图片等数据格式

在东莞控规查询系统中，关键接口为 identify 操作，其功能是返回地图上指定位置处要素的属性信息。该接口的典型请求URL结构如下：

http://[server]:6080/arcgis/rest/services/[serviceName]/[mapName]/MapServer/identify?[parameters]

1.2 接口参数深度解析

通过分析原始请求，我们提取出9个核心参数，其中两个为动态变化：

参数名	示例值	说明	是否必需
geometry	12710257,2593523	查询点坐标（Web墨卡托投影）	是
geometryType	esriGeometryPoint	几何类型（固定为点）	是
layers	visible:20	查询图层ID（控规为20）	是
tolerance	2	点击容差（像素）	是
mapExtent	12621272,2591650,12715908,2649830	地图范围（xmin,ymin,xmax,ymax）	是
imageDisplay	1115,929,96	图像显示参数（宽,高,DPI）	可选
f	json	返回格式（固定为JSON）	是
token	[动态令牌]	访问令牌	是
sr	3857	空间参考（Web墨卡托）	是

关键发现 ：实际测试表明， imageDisplay 参数对数据获取影响不大，可以简化请求URL。而 token 参数虽然看似复杂，但在当前系统中保持固定不变。

2. 空间网格划分策略

2.1 东莞全域范围确定

通过多次点击测试，我们确定了东莞市的Web墨卡托坐标范围：

x_min, y_min = 12621272, 2591650  # 西南角坐标
x_max, y_max = 12715908, 2649830  # 东北角坐标

该范围覆盖东莞全部行政区划，总面积约2460平方公里。在Web墨卡托投影（EPSG:3857）下，1单位≈1米，因此实际地理跨度约为94.6km×58.2km。

2.2 网格划分算法实现

为确保数据采集无遗漏，我们采用等间距网格法覆盖整个区域。核心算法如下：

import numpy as np

def generate_grid_centers(x_min, y_min, x_max, y_max, step=500):
    """
    生成等间距网格中心点坐标
    :param step: 网格间距（单位：米）
    :return: 中心点坐标列表[(x1,y1), (x2,y2), ...]
    """
    x_coords = np.arange(x_min, x_max + step, step)
    y_coords = np.arange(y_min, y_max + step, step)
    
    # 计算每个网格的中心点
    x_centers = (x_coords[:-1] + x_coords[1:]) / 2
    y_centers = (y_coords[:-1] + y_coords[1:]) / 2
    
    # 生成所有中心点组合
    centers = [(x, y) for x in x_centers for y in y_centers]
    return centers

当设置500米间隔时，将生成约36972个查询点（实际数量会根据精确范围略有变化）。网格密度选择需要考虑：

• 密度过大 ：增加请求次数，可能触发反爬机制
• 密度过小 ：导致数据遗漏，尤其对小地块影响明显

经过实测，在东莞地区：

• 500米间隔：约36,000次请求，可能遗漏部分小地块
• 300米间隔：约100,000次请求，覆盖率>99%
• 200米间隔：约225,000次请求，覆盖率接近100%

3. 完整Python实现方案

3.1 核心代码结构

我们构建了一个模块化的数据采集方案，主要包含以下组件：

dg_kgc/
├── config.py       # 配置参数
├── spatial.py      # 空间计算工具
├── api_client.py   # API请求处理
├── data_processor.py # 数据处理
└── main.py         # 主程序

3.1.1 API请求封装

# api_client.py
import requests
import time
import random

class KGCApiClient:
    def __init__(self, base_url, token, headers=None):
        self.base_url = base_url
        self.token = token
        self.headers = headers or {
            'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36',
            'Referer': 'http://120.86.191.153/'
        }
        
    def identify(self, point, layers='visible:20', tolerance=2, retry=3):
        """
        执行identify请求
        :param point: (x,y)坐标元组
        :param layers: 查询图层
        :param tolerance: 容差像素
        :param retry: 重试次数
        :return: JSON响应数据或None
        """
        params = {
            'geometry': f"{point[0]},{point[1]}",
            'geometryType': 'esriGeometryPoint',
            'layers': layers,
            'tolerance': tolerance,
            'mapExtent': '12621272,2591650,12715908,2649830',
            'f': 'json',
            'token': self.token,
            'sr': '3857'
        }
        
        for attempt in range(retry):
            try:
                resp = requests.get(
                    self.base_url,
                    params=params,
                    headers=self.headers,
                    timeout=(3, 10)
                )
                resp.raise_for_status()
                return resp.json()
            except Exception as e:
                if attempt == retry - 1:
                    print(f"请求失败: {e}")
                    return None
                sleep_time = random.randint(5, 15)
                time.sleep(sleep_time)

3.1.2 数据处理模块

# data_processor.py
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon

class DataProcessor:
    def __init__(self):
        self.attributes_df = None
        self.geometries_gdf = None
        
    def process_response(self, response_data):
        """处理API返回的JSON数据"""
        if not response_data or 'results' not in response_data:
            return False
            
        results = response_data['results']
        if not results:
            return False
            
        # 处理属性数据
        att_list = [r['attributes'] for r in results]
        temp_att = pd.DataFrame(att_list)
        
        # 处理几何数据
        geo_list = []
        for r in results:
            rings = r['geometry']['rings']
            if len(rings) == 1:
                geo = Polygon(rings[0])
            else:
                geo = Polygon(rings[0], rings[1:])
            geo_list.append({'OBJECTID': r['attributes']['OBJECTID'], 'geometry': geo})
            
        temp_geo = gpd.GeoDataFrame(geo_list, crs='EPSG:3857')
        
        # 合并到主数据集
        self.attributes_df = temp_att if self.attributes_df is None else pd.concat([self.attributes_df, temp_att])
        self.geometries_gdf = temp_geo if self.geometries_gdf is None else pd.concat([self.geometries_gdf, temp_geo])
        return True

3.2 主程序流程

# main.py
from config import API_URL, TOKEN
from spatial import generate_grid_centers
from api_client import KGCApiClient
from data_processor import DataProcessor
import time

def main():
    # 初始化组件
    client = KGCApiClient(API_URL, TOKEN)
    processor = DataProcessor()
    
    # 生成查询点网格
    centers = generate_grid_centers(12621272, 2591650, 12715908, 2649830, 500)
    total = len(centers)
    
    # 遍历所有网格点
    for idx, center in enumerate(centers, 1):
        print(f"进度: {idx}/{total} ({idx/total:.1%})", end='\r')
        
        # 发送请求并处理响应
        data = client.identify(center)
        if data:
            processor.process_response(data)
            
        # 礼貌性延迟
        time.sleep(0.2 + random.random())
    
    # 保存结果
    if processor.attributes_df is not None:
        processor.attributes_df.drop_duplicates().to_csv('dg_kg_attributes.csv', index=False, encoding='utf-8-sig')
    if processor.geometries_gdf is not None:
        processor.geometries_gdf.drop_duplicates().to_file('dg_kg_geometries.shp', encoding='utf-8')

if __name__ == '__main__':
    main()

4. 工程化优化与注意事项

4.1 性能优化策略

针对大规模数据采集，我们实施了多项优化：

1. 请求批处理 ：将多个点坐标合并为单个请求（需服务端支持）
2. 失败重试机制 ：指数退避算法处理临时故障
3. 本地缓存 ：避免重复获取已成功数据
4. 并行处理 ：使用 concurrent.futures 实现可控并发

# 优化后的请求处理示例
from concurrent.futures import ThreadPoolExecutor, as_completed

def batch_identify(points, client, max_workers=4):
    """批量查询多个点"""
    with ThreadPoolExecutor(max_workers=max_workers) as executor:
        futures = {executor.submit(client.identify, p): p for p in points}
        for future in as_completed(futures):
            point = futures[future]
            try:
                data = future.result()
                yield point, data
            except Exception as e:
                print(f"查询失败 {point}: {e}")

4.2 数据质量控制

获取的原始数据需要经过严格验证：

1. 几何完整性检查 ：

   • 验证多边形是否闭合
   • 检查自相交几何体
   • 修复无效坐标

2. 属性一致性检查 ：

   • 关键字段非空验证
   • 值域范围检查
   • 逻辑关系验证

3. 空间拓扑检查 ：

   • 地块之间不应有重叠
   • 不应存在空白区域
   • 边界对齐检查

4.3 反爬虫策略应对

为避免触发服务端防护机制，建议：

• 设置合理请求间隔（0.2-1秒）
• 随机化User-Agent
• 使用代理IP池（如需高频访问）
• 遵守robots.txt协议
• 限制总请求量，避免对服务器造成过大压力

5. 数据应用与可视化

5.1 属性关联与数据整合

由于Shapefile的属性表限制，我们采用了分离存储策略：

1. 几何数据保存为Shapefile
2. 属性数据保存为CSV
3. 通过QGIS或ArcGIS进行关联

在QGIS中的关联步骤：

1. 加载Shapefile和CSV文件
2. 右键图层 → 属性 → 连接
3. 设置连接字段（如OBJECTID）
4. 验证连接结果

5.2 空间分析与可视化

获取的数据可用于多种分析场景：

• 用地结构分析 ：统计各类用地比例
• 开发强度分析 ：容积率空间分布
• 公共服务设施覆盖分析 ：学校、医院等服务半径
• 三维城市建模 ：结合建筑高度数据

使用Python进行快速可视化的示例：

import matplotlib.pyplot as plt

def plot_landuse(gdf):
    """绘制用地类型分布图"""
    fig, ax = plt.subplots(figsize=(12, 10))
    
    # 按用地类型分类绘制
    for category, color in [
        ('R', 'yellow'),    # 居住
        ('C', 'red'),       # 商业
        ('M', 'brown'),     # 工业
        ('G', 'green'),     # 绿地
        ('S', 'gray')       # 道路
    ]:
        subset = gdf[gdf['LANDUSE'] == category]
        subset.plot(ax=ax, color=color, label=category)
    
    ax.legend()
    ax.set_title('东莞市控制性详细规划用地分布')
    plt.tight_layout()
    plt.show()

5.3 数据导出与共享

为方便不同平台使用，可将数据导出为多种格式：

• GeoJSON ：Web地图应用
• KML ：Google Earth展示
• PostGIS ：数据库存储
• FlatGeobuf ：高效传输格式

使用GDAL进行格式转换的示例命令：

# 转换为GeoJSON
ogr2ogr -f GeoJSON dg_kg.geojson dg_kg_geometries.shp

# 转换为KML
ogr2ogr -f KML dg_kg.kml dg_kg_geometries.shp