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

发布于 2020-9-3 18:07
浏览
0收藏

北京的详细地理位置数据我用Python获取了-开源基础软件社区在查找北京地图时,我发现很难找到准确的城镇和街道数据。 我偶尔看到有人分享有关天地的一些信息,其中包含精确到城镇和街道的行政区划数据。 因此,我尝试查找,处理和学习Python。

 

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

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

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

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

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

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

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

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

 

北京的详细地理位置数据我用Python获取了-开源基础软件社区

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

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

 

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

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

北京的详细地理位置数据我用Python获取了-开源基础软件社区

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

请求数据源并得到完整的 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")
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 来查看,发现了这么几个问题:

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

 

标签
已于2020-9-3 18:07:03修改
收藏
回复
举报
回复
添加资源
添加资源将有机会获得更多曝光,你也可以直接关联已上传资源 去关联
    相关推荐