东莞控制性详细规划数据获取

Atom数据

已于 2023-06-01 20:16:02 修改

阅读量801

点赞数 2

分类专栏：城市规划 GIS 文章标签： arcgis python 经验分享其他

于 2023-06-01 20:08:09 首次发布

本文链接：https://blog.csdn.net/silvertao123/article/details/130994444

版权

城市规划同时被 2 个专栏收录

5 篇文章 14 订阅

订阅专栏

GIS

3 篇文章 0 订阅

订阅专栏

有网友后台私信了东莞控规成果查询系统网站http://120.86.191.153/DG.GisClient/default.aspx，希望能获得该网站的控规数据。

首先按F12，分析并找到控规数据查询接口。如下图，点击一个地块，即出现该地块的基础信息面板，同时在Network里，可以看到数据接口http://120.86.191.153:6080/arcgis/rest/services/KGCX/KGCX/MapServer/identify?......（后接各种参数）。

接口返回的数据里，有地块的基本信息attributes，也有地块的几何信息geometry，通过循环，把整个东莞区域内的地块“点击”一遍，就可以获得全部数据。

然后分析接口里各个参数，通过删减的方法，得到必须保留的参数，如下：

(1)geometry=12710257,2593523 # 请求（点击）的点坐标

(2)geometryType=esriGeometryPoint

(3)layers=visible:1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35 # 图层ID，经分析，15为行政区划，20为控规图层，这里只需要控规图层，仅保留20即可

(4)tolerance=2

(5)mapExtent=12708661,2592799,12710891,2594657

(6)imageDisplay=1115,929,96 # 不知什么含义

(7)f=json # 获得json数据

(8)token=smLb_QfiJj_0mBbVF32JNFpCvn2PlHdRTRNRkwpRV_HB26u84puTMkLZWYjXTKbk8cucgJVlATcJUX9KvqM5Ew..

(9)sr=3857

上面9个参数里，只有两个参数是变动的：geometry和mapExtent。我们直接把mapExtent设为整个东莞市的范围，则只有一个变动的参数geometry。

通过点击地图，可以大致得到东莞市的范围mapExtent:12621272,2591650,12715908,2649830。

接着就是解决geometry参数了，我们将整个东莞市mapExtent划分为大小均等的网格，以500米为例，然后获得每个网格中心的坐标，作为我们需要的geometry参数。

代码如下：

首先导入所需模块

import pandas as pd
import geopandas as gpd
import numpy as np
import requests
from shapely.geometry import Polygon
import time
import random

然后通过Numpy划分网格并获得网格中心坐标

def loc():
    x1, y1 = 12621272, 2591650
    x2, y2 = 12715908, 2649830
    delta = 500 # 500米间隔
    x = np.arange(x1, x2 + delta, delta)
    y = np.arange(y1, y2 + delta, delta)
    x_center = (x[:-1] + x[1:]) / 2
    y_center = (y[:-1] + y[1:]) / 2
    center_list = []
    for i in range(len(x_center)):
        for j in range(len(y_center)):
            center_list.append((x_center[i], y_center[j]))
    return center_list

通过网格中心坐标获取数据。可能由于属性表字数限制的原因，直接将属性保存到shapefile中一直失败，因此这里把属性表和shapefile分别保存，最后在arcgis中连接属性表即可。

def getshp(p):
    try:
        global gl_df, gl_dt
        url = 'http://120.86.191.153:6080/arcgis/rest/services/KGCX/KGCX/MapServer/identify?geometry={},{}&geometryType=esriGeometryPoint&layers=visible:20&tolerance=2&mapExtent=12621272,2591650,12715908,2649830&imageDisplay=1269,929,96&f=json&token=smLb_QfiJj_0mBbVF32JNFpCvn2PlHdRTRNRkwpRV_HB26u84puTMkLZWYjXTKbk8cucgJVlATcJUX9KvqM5Ew..&sr=3857'.format(p[0],p[1])
        headers = {
            'Cookie': 'AGS_ROLES="LEqVyEWpMB/+iFTxb6RIUQ=="; ASP.NET_SessionId=sfdxtj5lrm2z5zk4pvttpy4o',
            'Referer': 'http://120.86.191.153/',
            'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/113.0.0.0 Safari/537.36'
        }
        r = requests.get(url, headers=headers, timeout=(3,7))
        data = r.json()['results']
        if len(data) > 0:
            print('有数据')
            for i in data:
                att = i['attributes']
                coords = i['geometry']['rings']
                if len(coords) == 1:
                    geo = Polygon(coords[0])
                else:
                    geo = Polygon(coords[0], [coords[1]])
                dt = pd.DataFrame([att])
                print(dt)
                if gl_df is None:
                    gl_df = dt
                else:
                    gl_df = pd.concat([gl_df, dt], ignore_index=True)
                dm = dt[['OBJECTID']]
                dm['geometry'] = geo
                print(dm)
                if gl_dt is None:
                    gl_dt = dm
                else:
                    gl_dt = pd.concat([gl_dt, dm], ignore_index=True)
        else:
            print('无数据')
    except:
        print('error')
        time.sleep(random.randint(5,10)) # 频繁的请求会对网站造成压力，有时没有反应，暂停一段时间再继续
        getshp(p)

最后执行几个函数

if __name__ == "__main__":
    gl_df = None
    gl_dt = None
    l = loc()
    m = 1
    for p in l:
        print('{}/{}'.format(m,len(l))) # 查看进度
        getshp(p)
        m += 1
        time.sleep(0.2)
    df = gl_df.drop_duplicates()
    df.to_csv('./dgkg.csv', index=False, encoding='utf-8-sig') # 设置好保存路径
    dt = gl_dt.drop_duplicates()
    output = gpd.GeoDataFrame(dt, crs='EPSG:3857') # 返回数据中有指明坐标系为3857
    output.set_geometry('geometry')
    output.to_file('./shp/dgkg.shp', driver='ESRI Shapefile', encoding='utf-8') # 设置好保存路径

以500米间隔为例，获得的数据共36972条。