在线地图获取城市路网数据

在线地图获取城市路网数据

近期科研项目中，需要获取城市路网数据，于是详细阅读各大在线地图api接口，总结出来这么一条可行的思路：

1. 首先获取城市轮廓
2. 根据城市轮廓把城市分割成若干个小块
3. 在每个小块中根据在线地图的POI检索接口，检索小块中的道路，获取道路清单
4. 根据道路清单去获取道路经纬度数据

用到以下技术栈：

• Python程序开发
• 百度、高德在线地图接口
• 百度提供的JavaScript API GL

一、技术问题

针对这个问题，截止到当前日期（2023-09-25），有下面这几个问题

1. 获取城市轮廓需要用到地图的行政区划边界查询接口，目前百度在线地图不提供，高德在线地图提供接口，接口地址：

   URL	https://restapi.amap.com/v3/config/district?parameters
   请求方式	GET

2. 百度在线地图的多边形区域检索为高级权限，需要提交工单申请，但十有八九是不提供免费服务，高德在线地图的搜索POI接口提供多边形搜索接口服务

3. 百度在线地图的POI编码中不包含道路，高德有一项为道路

综上，高德完胜。但无奈项目中用的是百度地图，所以无形中多了一步坐标转换的过程，而且在线地图对个人开发者的免费额度越来越少，只能将就先用了

二、技术实现

1、城市行政区划边界查询

这一步其实我个人是用百度离线地图来实现的，但既然高德提供了免费的接口，我也试了一下

import requests

url = 'https://restapi.amap.com/v3/config/district'
key = '自己去申请ak'

params = {
    'key': key,
    'keywords': '无锡',
    'output': 'JSON',
    'extensions': 'all'
}
response = requests.get(url=url, params=params)
if response:
    res_json = response.json()
    boundary = res_json['districts'][0]['polyline']
    boundary_list = boundary.split(';')
    boundaries = []
    for bound in boundary_list:
        item = bound.split(',')
        boundaries.append([float(item[0]), float(item[1])])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

Python+json的处理，注意一下代码中输出数据的格式

得到的boundaries如下：

2、根据行政区划边界，切分网格

思路是根据行政区划边界的最大值和最小值，设定一个步长，根据步长来划分网格，下面是Python实现代码

def divide_region_by_step(coordinates, step):
    # 获取最小和最大经纬度
    min_lon = min(coord[0] for coord in coordinates)
    max_lon = max(coord[0] for coord in coordinates)
    min_lat = min(coord[1] for coord in coordinates)
    max_lat = max(coord[1] for coord in coordinates)

    # 计算经纬度方向上的步长
    # lon_step = step / 111.0  # 经度每度大概111公里
    # lat_step = step / 111.0  # 纬度每度大概111公里
    lon_step = step
    lat_step = step

    # 划分小区域
    regions = []
    current_lat = min_lat
    while current_lat < max_lat:
        current_lon = min_lon
        while current_lon < max_lon:
            region = {
                'min_lon': current_lon,
                'max_lon': current_lon + lon_step,
                'min_lat': current_lat,
                'max_lat': current_lat + lat_step
            }
            regions.append(region)
            current_lon += lon_step
        current_lat += lat_step

    return regions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

这样切分不太精细，可以看一下切分的效果：

黑色的是无锡市的行政区划图，蓝色的是我画的网格，可以看到，左上角和右下角有大部分湖州的地方和苏州的地方也画了网格，我觉得边界以外的应该去掉，应该有解决办法，不过不想弄了

3、获取网格内的道路清单

我使用的是高德的搜索POI中的多边形搜索，请求地址为：

接口地址中的polygon其实就是2中矩形的对角的两个点

需要的参数是接口文档中给的一个查询POI类型，我高德提供的POI分类编码表下下来，发现代码190301表示道路名，所以就直接用这个参数了

来看代码：

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# author:HP
# datetime:2023/9/22 14:18

# encoding:utf-8
import pandas as pd
import requests

# 读取行政区划边界文件
content = pd.read_excel('wuxi_positions_amap.xlsx', engine='openpyxl')
positions = content.values.tolist()


def divide_region_by_step(coordinates, step):
    # 获取最小和最大经纬度
    min_lon = min(coord[0] for coord in coordinates)
    max_lon = max(coord[0] for coord in coordinates)
    min_lat = min(coord[1] for coord in coordinates)
    max_lat = max(coord[1] for coord in coordinates)

    # 计算经纬度方向上的步长
    # lon_step = step / 111.0  # 经度每度大概111公里
    # lat_step = step / 111.0  # 纬度每度大概111公里
    lon_step = step
    lat_step = step

    # 划分小区域
    regions = []
    current_lat = min_lat
    while current_lat < max_lat:
        current_lon = min_lon
        while current_lon < max_lon:
            region = {
                'min_lon': current_lon,
                'max_lon': current_lon + lon_step,
                'min_lat': current_lat,
                'max_lat': current_lat + lat_step
            }
            regions.append(region)
            current_lon += lon_step
        current_lat += lat_step

    return regions


step = 0.1
results = divide_region_by_step(positions, step)

url = 'https://restapi.amap.com/v3/place/polygon'
key = '自己去申请高德的key'
roads = []
for i in range(len(results)):
    print(str(i) + '-------------')
    polygon = str(results[i]['min_lon']) + ',' + str(results[i]['min_lat']) + '|' + str(
        results[i]['max_lon']) + ',' + str(results[i]['max_lat'])
    params = {
        "polygon": polygon,
        "types": '190301',
        "key": key,
    }
    response = requests.get(url=url, params=params)
    if response:
        response_json = response.json()
        for j in response_json['pois']:
            print(j['name'] + '-------------')
            district = j['adname']
            road = j['name']
            location = j['location']
            roads.append([district, road, location])
df_roads = pd.DataFrame(roads, columns=['行政区划', '路名', '经纬度'])
df_roads.to_excel('wuxi_roads.xlsx', index=False)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

注意我读的文件就是第一步中获取的城市行政区划边界，看一下导出的数据：

其实效果不太好，有个坑需要注意下，接口中的参数：

params = {
        "polygon": polygon,
        "types": '190301',
        "key": key,
    }
1
2
3
4
5

只提供了秘钥、道路POI代码、网格坐标，其实还有一个offset我用的默认的，就是这个网格中返回的道路数据的条数，它默认是20，而且文档中标注强烈建议不超过25，若超过25可能造成访问报错，我看了一下我的运行结果，它最后生成的数据是1970条，而我的网格有99个，99*20=1980，也就是说，几乎每个网格中都返回了20条道路，我并没有一个个去考证网格中的数据是不是对的，但是就这个数据而言，也就是说，99个网格，总共只少了10条路，换句话说，几乎所有的网格中的数据应该是大于20条的，但因为限制，只返回了20条数据。所以返回的数据是不够的，因此，网格应该加密，加密的方法就是把划分网格的步长调小，但是由于高德限额的限制，我没有去尝试了。

至此已经获取到道路清单了

4、道路经纬度数据

按理说，这是最重要的，但是目前没有任何在线地图提供相关的接口

我最后使用离线地图获取的，给代码也没有意义了

三、JavaScript API GL的使用

技术实现的时候，我用了百度地图的JavaScript API，其实最后用到项目中也要使用JavaScript，不少小伙伴还只会Python，所以我提供一下我的html文件代码

DOCTYPE html>
<html>
<head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no"/>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
    <title>Baidu Map title>
    <style type="text/css">
        html {
            height: 100%
        }

        body {
            height: 100%;
            margin: 0px;
            padding: 0px
        }

        #container {
            height: 100%
        }
    style>
    <script type="text/javascript"
            src="https://api.map.baidu.com/api?v=1.0&type=webgl&ak=自己申请的百度地图AK">

    script>
head>
<body>
<div id="container">div>

<script>
    const points = [[[119.527015, 31.11083], [119.627015, 31.21083]]] // 这里用自己的网格数据
    // console.log(points.length)
    const map = new BMapGL.Map("container");           // 创建地图实例
    const point = new BMapGL.Point(120.3119, 31.4912);  // 创建点坐标
    map.centerAndZoom(point, 10);                 // 初始化地图，设置中心点坐标和地图级别

    map.enableScrollWheelZoom(true);     //开启鼠标滚轮缩放
    for (let i = 0; i < points.length; i++) {
        let pStart = new BMapGL.Point(points[i][0][0], points[i][0][1])
        let pEnd = new BMapGL.Point(points[i][1][0], points[i][1][1])
        let rectangle = new BMapGL.Polygon([
            new BMapGL.Point(pStart.lng, pStart.lat),
            new BMapGL.Point(pEnd.lng, pStart.lat),
            new BMapGL.Point(pEnd.lng, pEnd.lat),
            new BMapGL.Point(pStart.lng, pEnd.lat)
        ], {strokeColor: "blue", strokeWeight: 2, strokeOpacity: 0.5});  //创建矩形
        map.addOverlay(rectangle)
    }

    const boundary = [[120.141768, 31.301344], [120.107967, 31.281801]] // 这里用自己的行政区划边界数据
    const boundary_points = []
    for(let i = 0; i < boundary.length; i++){
        boundary_points.push(new BMapGL.Point(boundary[i][0], boundary[i][1]))
    }
    const polygon = new BMapGL.Polygon(boundary_points, {strokeColor:"black", strokeWeight:5, strokeOpacity:0.5})
    map.addOverlay(polygon)
script>
body>
html>


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

代码不解释了，主要是数据要换成自己的