深圳详规“一张图”数据获取【代码更新】

之前写了一篇深圳详规“一张图”数据获取的文章，后台不断有规划行业的小伙伴询问我代码运行不了的问题，或者索求数据。趁着这次得闲，更新一下代码。

距上次获取数据已经过了一年多了，网站的数据入口还是没有变化。但是在首次打开网站详规“一张图” (sz.gov.cn)的查询须知小窗口中有注明“每月10日更新一版详细规划‘一张图’公众版数据”，所以可以把代码保存，需要更新的时候再跑一次。

具体过程参考原来的文章，新的代码如下：

import requests
import pandas as pd
import geopandas as gpd
from shapely import Polygon
import time


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

def getData(rect):
    try:
        url = 'http://pnr.sz.gov.cn/d-suplicmap/dynamap_1/rest/services/NET_XGYZT/MapServer/0/query'
        header = {
            # 'Cookie': 'AD_insert_cookie=72221312',
            'Host': '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/110.0.0.0 Safari/537.36 Edg/110.0.1587.63',
            'X-Requested-With': 'XMLHttpRequest'
            }
        params = {
            'geometry': rect,
            'geometryType': 'esriGeometryEnvelope',
            'spatialRel': 'esriSpatialRelIntersects',
            'returnGeometry': 'true',
            'outFields': '*',
            'f': 'pjson',
            'where': '1=1'
            }
        r = requests.post(url, headers=header, data=params, timeout=(3,20))
        jdata = r.json()
        f = jdata['features']
        if len(f) > 0:
            df = None
            for a in f:
                dt = pd.DataFrame([a['attributes']])
                print(dt)
                rings = a['geometry']['rings'] # polygon存在无孔和有孔两种情况，其数据结构不同，需分开处理
                if len(rings) == 1: # polygon无孔的情形
                    l = []
                    for m in rings[0]:
                        l.append((m[0],m[1]))
                    dt['geometry'] = Polygon(l)
                else: # polygon有孔的情形
                    a = []
                    for m in rings[0]:
                        a.append((m[0],m[1]))
                    b = []
                    for n in rings[1:]:
                        c = []
                        for k in n:
                            c.append((k[0],k[1]))
                        b.append(tuple(c))
                    dt['geometry'] = Polygon(a,b)
                if df is None:
                    df = dt
                else:
                    df = pd.concat([df,dt],ignore_index=True)
            return df
        else:
            print('该图幅无数据')
    except:
        print('error...try again')
        time.sleep(3)
        getData(rect)

if __name__ == '__main__':
    print('开始爬取深圳详规一张图数据')
    print('==========================================================')
    loc = LocaDiv('473934.2251,2478026.4717,564199.7618,2529560.8922', 1500)
    locs_to_use = loc.ls_row()
    n = 1
    df = None
    for loc_to_use in locs_to_use:
        print('获取第{}图幅的数据'.format(n))
        print(loc_to_use)
        dt = getData(loc_to_use)
        if df is None:
            df = dt
        else:
            df = pd.concat([df,dt],ignore_index=True)
        n += 1
        time.sleep(1)
    dd = df.drop_duplicates() # 去重
    output = gpd.GeoDataFrame(dd, crs='EPSG:4547') # 返回数据中有指明坐标系为4547
    output.set_geometry('geometry')
    output.to_file('./shp/深圳详规一张图.shp', driver='ESRI Shapefile', encoding='utf-8')
    print('数据爬取完毕')

代码运行结束后，即可得到完整的深圳市详细规划一张图公众版数据。在ArcGIS中打开，发现无法显示属性表。没关系，将数据再导出一次即可。

 共有45470个面数据。