北京的详细地理位置数据我用Python获取了

在查找北京地图时，我发现很难找到准确的城镇和街道数据。 我偶尔看到有人分享有关天地的一些信息，其中包含精确到城镇和街道的行政区划数据。 因此，我尝试查找，处理和学习Python。

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

有两个「八里庄街道」
1. 海淀一个，北京市海淀区北洼路64号

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

好几个乡镇、街道没有在「天地图」数据中标出
1. 「长辛店镇」 （只有一个「长辛店街道」）

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

北京存在好几个 Multiple Polygon
1. 「首都机场」是属于朝阳区的

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

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

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

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

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

各区域的数据都记录在 rows 中
各区域的行政区划的边界点位信息都记录在 points 中

根据数据特点，接下来要做的是：

请求数据源并得到完整的 json 文件
从 json 文件中得到各区域的数据，即提取出 rows 字段
在从各区域的数据中提取出其边界点位信息，即 points 字段
最后将提取出的信息输出为 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")
>>> 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 来查看，发现了这么几个问题：

部分街道的数据缺失，包括：「田村路街道」、「永定路街道」和「万寿路街道」。这部分需要将街道数据与北京市区县数据进行 difference 处理，再将 difference 的结果与街道数据进行 union，最后根据大致的街道界限手动来对 difference 那部分进行分割。
部分街道的数据被占用，这部分是「长辛店镇」的地理位置是被「长辛店街道」给完全占用了。需要根据大致的街道界限手动来对「长辛店街道」进行分割。