写点什么

通过 Python 来获取北京市乡镇、街道行政区划数据

用户头像
Puran
关注
发布于: 2020 年 06 月 28 日
通过Python来获取北京市乡镇、街道行政区划数据

在查找北京市地图时,发现精确到乡镇、街道的数据基本很难找到。偶然间看到有人分享关于天地图的一些信息,里面包含有精确到乡镇、街道的行政区划数据。于是便尝试去查找、处理下,并学习一下Python。

1. 北京市行政区划的几个有意思的点

首先,把我在北京市乡镇、街道行政区划数据处理中发现的几个有意思的地方,在这里分享一下。这部分时参考了「北京·天地图」和「国家卫健委疫情风险等级查询系统」。

  1. 有两个「八里庄街道」

1. 海淀一个,北京市海淀区北洼路64号

2. 朝阳一个,北京市朝阳区柴家湾12号

  1. 好几个乡镇、街道没有在「天地图」数据中标出

1. 「长辛店镇」 (只有一个「长辛店街道」)

2. 连接在一起的「田村路街道」、「永定路街道」和「万寿路街道」

  1. 北京存在好几个 Multiple Polygon

1. 「首都机场」是属于朝阳区的

2. 「燕园街道」等是multiple polygon

  1. 「南苑乡」在地理位置上是包含「南苑街道」、「大红门街道」的,但是从行政区划上,这两者又都是并列的关系。这个有点搞不清楚是怎么划分和管理的。



下面来详细记录下数据获取和处理的过程。

2. 数据源

「北京·天地图」有一个入口,提供北京市行政区划的查询,里面有精确到乡镇、街道的数据。通过浏览器的开发者模式来查看网络请求,能发现在请求的数据中,有关于各街道的详细的 json 信息,当然还包括各区域的边界线的点坐标。



通过简单的分析能发现完整的 json 信息请求的 URL 为http://beijing.tianditu.gov.cn:8080/dfc/services/sgssfs/2220?request=getfeature。从 json 数据中,我们可以了解到:

  • 各区域的数据都记录在 rows

  • 各区域的行政区划的边界点位信息都记录在 points



根据数据特点,接下来要做的是:

  1. 请求数据源并得到完整的 json 文件

  2. 从 json 文件中得到各区域的数据,即提取出 rows 字段

  3. 在从各区域的数据中提取出其边界点位信息,即 points 字段

  4. 最后将提取出的信息输出为 shapefile



3. 数据处理

由于我们的数据是 json 格式的数据,于是要用到 Python 的 json 模块。最终要得到的数据格式是 shapefile 等地理数据,这里就需要借助 Geopandas。



3.1. 获取 json 文件并将其转换为 DataFrame

  • 通过 requests 模块来请求数据,并将其序列化为 json:

beijing_street_info = requests.get(online_service).json()

  • 从 json 中提取出 rows 字段

detail_info = beijing_street_info["rows"]

  • 将 json 序列化成 Pandas 的 DataFrame,这里需要from pandas.io.json import json_normalize

street_df = json_normalize(detail_info)



3.2. 提取 points 中的坐标

这部分是数据处理的重点。由于 points 是字符串,以逗号,分隔单个点位的经纬度,以分号;分隔不同的点。

另外,由于北京市存在好几个行政区都包含多个不相邻的面,即存在Multiple Polygon的情况,在points中是以竖线 | 来分隔的。这部分一开始没能想到,是在后续处理数据时才发现的。

这是一个例子(该数据并不是真实数据):

"116.4535,39.9614;116.4529,39.9609;116.4519,39.9600;116.4512,39.9595;116.4535,39.9614|116.4508,39.9590;116.4502,39.9586;116.4498,39.9581;116.4480,39.9566;116.4508,39.9590"



这部分的数据提取最主要是需要将作为一长串的 points(string)提取出经纬度坐标对(tuple)出来。主要有这么几步:

  • 根据数据特点,对包含多面的数据先以竖线 | 来分割points字符串,得到一个包含两个字符串的列表

>>> points_polygon = points.split("|")
>>> points_polygon
['116.4535,39.9614;116.4529,39.9609;116.4519,39.9600;116.4512,39.9595;116.4535,39.9614', '116.4508,39.9590;116.4502,39.9586;116.4498,39.9581;116.4480,39.9566;116.4508,39.9590']



  • 再对列表中的字符串以分号 ; 分割,得到一个包含两个列表的列表,这两个列表分别包含了两个多边形内的点。这里需要注意的是,列表是没有 split() 函数的,因此这里需要用一个列表推导式(List Comprehensions)来对列表中的字符串进行分割。

>>> points_list = list(points_region.split(";") for points_region in points_polygon)
>>> points_list
[['116.4535,39.9614', '116.4529,39.9609', '116.4519,39.9600', '116.4512,39.9595', '116.4535,39.9614'], ['116.4508,39.9590', '116.4502,39.9586', '116.4498,39.9581', '116.4480,39.9566', '116.4508,39.9590']]



  • 现在所有的点位都已经被分割出来了,但还是字符串的形式,接下来要做的就是将字符串的点位转为经纬度坐标对(tuple)。

3.3. 通过eval 和 map 将字符串的坐标转换为坐标对

  • eval() 函数

我们可以先看看Python内置的 eval() 函数,该函数可以把字符串当作表达式执行并返回一个结果,这一特性可以直接拿来用在我们的点位字符串转为经纬度坐标对上。

>>> eval("2 + 3")
5
>>> eval('116.4535,39.9614')
(116.4535, 39.9614)



当只有一个列表时,可以直接应用列表推导式。

>>> list1 = ['116.4535,39.9614', '116.4529,39.9609', '116.4519,39.9600', '116.4512,39.9595', '116.4535,39.9614']
>>> coor_pair = [eval(coor_string) for coor_string in list1]
>>> coor_pair
[(116.4535, 39.9614), (116.4529, 39.9609), (116.4519, 39.96), (116.4512, 39.9595), (116.4535, 39.9614)]



但是我们上面的列表内部又有列表,由于 eval() 的参数只能是string、char和code,如果直接使用列表推导式的话会出错,这里需要做一些特殊处理。引入 Python 内置的 map() 函数。



  • map() 函数

map() 会根据提供的函数对指定序列做映射,其语法为:map (function, iteration, …)

def addition(n):
return n + n
# We double all numbers using map()
numbers = (1, 2, 3, 4)
result = map(addition, numbers)
print(list(result))



因此,需要对我们上面的 points_list 列表中的列表应用 eval() 函数:

>>> coor_pair_list = [list(map(eval, polygon_list)) for polygon_list in points_list]
>>> coor_pair_list
[[(116.4535, 39.9614), (116.4529, 39.9609), (116.4519, 39.96), (116.4512, 39.9595), (116.4535, 39.9614)], [(116.4508, 39.959), (116.4502, 39.9586), (116.4498, 39.9581), (116.448, 39.9566), (116.4508, 39.959)]]



3.4. 将坐标对转换为shapely.geometry

使用shapely时需要引入from shapely.geometry import Point, LineString, Polygon, MultiPolygon,点、线、面、多面分别为:

point = Point(2.2, 4.2)
line = LineString([point1, point2, point3])
poly = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])
multi_poly = MultiPolygon([poly1, poly2])



我们之前已经把点位数据转成了 shapely.geometry 需要的格式,直接使用就好。

  • 对 polygon,直接用 Polygon(coor_pair_list)

  • 若是 multiple polygon,则需要列表推导式:MultiPolygon(Polygon(polygon_list) for polygon_list in coor_pair_list)



以上都是对一个区域进行的数据处理,对所有的区域,我们需要用到循环,并将各个区域的数据追加到总的geometry 中。

polygon_geometry = []
for i in range(total):
points_string = street_df.points[i]
if "|" in points_string:
...
polygon_geometry.append(MultiPolygon(Polygon(polygon_list) for polygon_list coor_pair_list))
else:
...
polygon_geometry.append(Polygon(coor_pair_list))



3.5. 将 DataFrame 转为 GeoDataFrame

得到了 geometry 后,就可以将 DataFrame 转换为 GeoDataFrame 了:

street_gdf = gpd.GeoDataFrame(street_df, geometry=polygon_geometry)



4. 保存

最后,根据数据的特点,按照某一字段进行分类后导出,得到完整的北京市街道图。

在分类时,使用的是 groupby() 函数,导出为 shapefile 时需要注意编码。天地图采用的坐标系统是CGCS 2000,在导出时可以直接指定坐标系统,也可以在后续的 GIS 软件中指定。

groupby_PRODATE = street_gdf.groupby("PRODATE")
for key, group in groupby_PRODATE:
...
group.to_file(output_path, encoding='utf-8')



5. 数据后处理

在 Python 中,也是可以直接通过 matplotlib 来展示 shapely.geometry 数据的,但是能支持的功能有限,在我看来还不是很习惯,便还是用传统的 GIS 工具来对导出后的数据进行一些可视化的处理。

数据导出后,通过 QGIS 来查看,发现了这么几个问题:

  1. 部分街道的数据缺失,包括:「田村路街道」、「永定路街道」和「万寿路街道」。这部分需要将街道数据与北京市区县数据进行 difference 处理,再将 difference 的结果与街道数据进行 union,最后根据大致的街道界限手动来对 difference 那部分进行分割。

  2. 部分街道的数据被占用,这部分是「长辛店镇」的地理位置是被「长辛店街道」给完全占用了。需要根据大致的街道界限手动来对「长辛店街道」进行分割。



6. 参考



发布于: 2020 年 06 月 28 日阅读数: 231
用户头像

Puran

关注

GIS从业者,正在往开发的路上小跑。 2018.03.29 加入

从业4年的GIS开发小白,work@esri。

评论

发布
暂无评论
通过Python来获取北京市乡镇、街道行政区划数据