深圳详规“一张图”数据获取

2021年4月19日，深圳市规划和自然资源局推出了“详细规划一张图”公众版（试运行）详规“一张图”。“公众版”上线内容包括深圳市已批法定图则及局部调整、城市更新、土地整备等规划成果，具体信息包含规划地块的位置、用地面积、规划功能、容积率、配套设施等。

 通过F12开发者工具调试页面一番摸索，一下子就找到了数据的入口http://pnr.sz.gov.cn/d-suplicmap/dynamap_1/rest/services/NET_XGYZT/MapServer/0/query?callback=jQuery331046554915921307893_1633935763478

 该入口为POST请求，通过点选地图上不同图斑，可以发现每次只需要改变请求参数的geometry，该参数为请求的矩形范围，如“512622.92051969835,2494047.3300745627,512622.92052169837,2494047.3300765622”，其中“512622.92051969835,2494047.3300745627”为左下角坐标，“512622.92052169837,2494047.3300765622”为右上角坐标。

返回的数据包含了图斑的所有信息，图斑的属性信息在attributes中，如地块编号DKBH，规划名称GHMC，用地性质YDXZ，容积率RJL，用地面积YDMJ等，图斑的几何信息在geometry中，图斑的坐标系信息在spatialReference中。

 如果我们把矩形范围调大，返回的数据量也随着变多，即该POST实际上是用一个矩形去框选图斑数据，返回在该矩形范围内或相交的图斑。而当返回的图斑数量到达1000个时，即使矩形范围继续增大 ，数据量也不会增加，且返回的数据在结构上也有变化。因此需要将单次请求的数据量控制在1000以内。再通过点选网站上深圳范围的左下角和右上角，可以得到涵盖深圳的矩形范围为：473934.2251,2478026.4717,564199.7618,2529560.8922。

经过一定测试，当矩形的边长为2000时，我们得到的数据量不会超过1000条，也不会因矩形过于小而造成需要很多次请求的结果。搞清楚了数据的来龙去脉，下面开始编写代码：

（1）导入需要使用的库

import requests
import json
import shapefile
import osr
import time

（2）编写划分网格的函数

class LocaDiv(object):
    def __init__(self, loc_all, divd):
        self.loc_all = loc_all
        self.divd = divd

    def lat_all(self):
        lat_sw = float(self.loc_all.split(',')[0])
        lat_ne = float(self.loc_all.split(',')[2])
        lat_list = [str(lat_ne)]
        while lat_ne - lat_sw >= 0:
            m = lat_ne - self.divd
            lat_ne = lat_ne - self.divd
            lat_list.append('%.2f' % m)
        return sorted(lat_list)

    def lng_all(self):
        lng_sw = float(self.loc_all.split(',')[1])
        lng_ne = float(self.loc_all.split(',')[3])
        lng_list = [str(lng_ne)]
        while lng_ne - lng_sw >= 0:
            m = lng_ne - self.divd
            lng_ne = lng_ne - self.divd
            lng_list.append('%.2f' % m)
        return sorted(lng_list)

    def ls_com(self):
        l1 = self.lat_all()
        l2 = self.lng_all()
        ab_list = []
        for i in range(0, len(l1)):
            a = str(l1[i])
            for i2 in range(0, len(l2)):
                b = str(l2[i2])
                ab = a + ',' + b
                ab_list.append(ab)
        return ab_list

    def ls_row(self):
        l1 = self.lat_all()
        l2 = self.lng_all()
        ls_com_v = self.ls_com()
        ls = []
        for n in range(0, len(l1) - 1):
            for i in range(len(l2) * n, len(l2) * (n + 1) - 1):
                a = ls_com_v[i]
                b = ls_com_v[i + len(l2) + 1]
                ab = a + ',' + b
                ls.append(ab)
        return ls

（3）编写获取数据的函数，变量为获取数据的矩形区域rect和该矩形的编号num

def getData(rect,num):
    url = 'http://pnr.sz.gov.cn/d-suplicmap/dynamap_1/rest/services/NET_XGYZT/MapServer/0/query?callback=jQuery331046554915921307893_1633935763478'
    header = {
        'Cookie': 'AD_insert_cookie=72221312',
        'Host': 'pnr.sz.gov.cn',
        'Origin': 'http://pnr.sz.gov.cn',
        'Referer': 'http://pnr.sz.gov.cn/d-xgmap/',
        'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/89.0.4389.114 Safari/537.36',
        'X-Requested-With': 'XMLHttpRequest'
        }
    params = {
        'geometry': rect,
        'geometryType': 'esriGeometryEnvelope',
        'spatialRel': 'esriSpatialRelIntersects',
        'returnGeometry': 'true',
        'outFields': '*',
        'f': 'pjson',
        'where': '1=1'
        }
    try:
        r = requests.post(url, headers=header, data=params, timeout=(3,20))
        jdata = json.loads(r.text.replace('\n','').replace(' ','').split('"features":')[1].split('})')[0])
        if len(jdata)>0:
            w = shapefile.Writer('./shp/shenzhen_{}.shp'.format(num))
            w.field('序号','C')
            w.field('来源','C')
            w.field('代码','C')
            w.field('规划名称','C')
            w.field('地块编号','C')
            w.field('用地性质','C')
            w.field('容积率','C')
            w.field('用地面积','C')
            for item in jdata:
                oid = item['attributes']['OBJECTID']
                source = item['attributes']['SOURCE']
                maincode = item['attributes']['MAINCODE']
                ghmc = item['attributes']['GHMC']
                dkbh = item['attributes']['DKBH']
                ydxz = item['attributes']['YDXZ']
                rjl = item['attributes']['RJL']
                ydmj = item['attributes']['YDMJ']
                rings = item['geometry']['rings']
                w.poly(rings)
                w.record(序号=oid,
                         来源=source,
                         代码=maincode,
                         规划名称=ghmc,
                         地块编号=dkbh,
                         用地性质=ydxz,
                         容积率=rjl,
                         用地面积=ydmj)
            w.close()
            #添加投影信息
            proj = osr.SpatialReference()
            proj.ImportFromEPSG(4547)
            wkt = proj.ExportToWkt()
            projfile = './shp/shenzhen_{}.prj'.format(num)
            f = open(projfile, 'w')
            f.write(wkt)
            f.close()
        else:
            print('该图幅无数据')
    except:
        print('error')

（4）将整个深圳范围划分为2000米边长的网格，获得整个深圳的详规“一张图”数据

if __name__ == '__main__':
    print('开始爬取深圳详规一张图数据')
    print('==========================================================')
    loc = LocaDiv('473934.2251,2478026.4717,564199.7618,2529560.8922', 2000)
    locs_to_use = loc.ls_row()
    n = 1
    for loc_to_use in locs_to_use:
        print('获取第{}图幅的数据'.format(n))
        getData(loc_to_use,n)
        n += 1
        time.sleep(1)
    print('数据爬取完毕')

 

最后将得到的所有图幅数据在ArcGIS软件中进行合并，并删除重复的数据，即可得到整个深圳市的详规“一张图”数据。