城市行政区边界

一、工作背景

1、数据分析从一个城市视角拓展到国内外主要城市的数据比较

2、互联网城市数据的对比分析

3、 百度慧眼，提供了迁徙数据，具体看下长三角的联系

二、准备工作

1、行政区划信息

中国民政部官网提供了每年更新的行政区划，2020年12年有完整的行政区划代码表，省、市、区县、街道、镇。数据下载地址：行政区划代码_中华人民共和国民政部

2、行政区边界信息

在网上搜了一圈后，采用阿里云datav数据可视化平台提供的边界信息，该边界信息来源于高德开放平台。网站地址：datav.aliyun.com/portal/school/atlas/area_selector?spm=a2crr.23498931.0.0.4ad815ddMDdXGW）

三、边界处理

1、数据下载

阿里云datav数据可视化平台提供json数据下载，通过python遍历国内333个地级市下载，python下载代码如下：

urlstr = 'https://geo.datav.aliyun.com/areas_v3/bound/'+ str(id2) +'_full.json'

req = urllib.request.Request(urlstr)

f = urllib.request.urlopen(req)

result = f.read().decode('utf-8')

res_data = json.loads(result)

下载得到境内31个省市、333个地级市的边界json，如下：

2、市、区边界拼合

上述下载的是每个省、每个城市的json，需要拼合为一个全国的省、市、区地理文件。

拼合方法是，逐个解析json文件进行拼合，共输出为三个文件

省级边界、市级边界、区县边界

重点说一下镂空面数据的处理，通过QGIS加载json查看镂空面的wkt，发现镂空面的multipolygon表示如下

（（（），（），（））），第一个黑括号为完整边界的坐标串，第二个、第三个镂空的面坐标串；

（（（）），（（）），（（））），这个常规的面表达。

3、处理程序

try:

req = urllib.request.Request(urlstr)

f = urllib.request.urlopen(req)

result = f.read().decode('utf-8')

res_data = json.loads(result)

rdatas = res_data["features"]

for rdata in rdatas:

polygons = rdata["geometry"]["coordinates"]

ptype = rdata["geometry"]["type"]

if ptype == 'Polygon' :

#print(rdata["properties"]["adcode"])

c1 = 1

wkt2 = '"MultiPolygon((('

for polygon in polygons:

if c1 == 1:

c1 = 2

c2 = 1

for wkt1 in polygon:

if c2==1:

c2 = 2

wkt2 += str(wkt1[0]) +' '+ str(wkt1[1])

else:

wkt2 += ','+ str(wkt1[0])+' '+str(wkt1[1])

wkt2 += ')'

else:

c2 = 1

for wkt1 in polygon:

if c2==1:

c2 = 2

wkt2 += ',('+str(wkt1[0]) +' '+ str(wkt1[1])

else:

wkt2 += ','+ str(wkt1[0])+' '+str(wkt1[1])

wkt2 += ')'

wkt2 += '))"'

elif ptype == 'MultiPolygon' : #or rdata["properties"]["adcode"] in(120113,340800,341700)

c1 = 1

wkt2 = '"MultiPolygon('

for polygon in polygons:

if c1 == 1:

c1 = 2

c2 = 1

for wkt10 in polygon:

for wkt1 in wkt10:

if c2==1:

c2 = 2

wkt2 += '(('+str(wkt1[0]) +' '+ str(wkt1[1])

else:

wkt2 += ','+ str(wkt1[0])+' '+str(wkt1[1])

wkt2 += '))'

else:

c2 = 1

for wkt10 in polygon:

for wkt1 in wkt10:

if c2==1:

c2 = 2

wkt2 += ',(('+str(wkt1[0])+' '+str(wkt1[1])

else:

wkt2 += ','+ str(wkt1[0])+' '+str(wkt1[1])

wkt2 += '))'

wkt2 += ')"'