规划技术

绘制一个真实的交通等时圈

使用百度批量算路 API 和 ArcGIS Python 工具绘制等时圈

近日因为导师要求,要做「交通时距」的分析,所以需要根据「真实的交通情景」来绘制从某一点出发到各个位置的时间。换言之,需要绘制真实的交通等时圈(这是规划常见的分析图之一)。

考虑到这篇文章的读者可能是对规划有兴趣的非专业人士,这里首先介绍一下等时圈可以用来干什么,为什么要画这个分析。

等时圈(isochrone),是指以某地点作为中心(行程的起点,origin),以某种出行方式(步行,骑车,公共交通,驾车等),在一定时间内(比如,15 分钟)能够到达的范围。也可以理解为出行时间的等高线(contour)。

它是交通耗时在地理空间上的反映,可以在一定程度上反映了该地区的交通便利程度。一般来说,做等时圈是做周边公服设施可达性分析,类似社区 15 分钟生活圈,看周边有无超市,买菜,吃饭,娱乐……以及可以用于商业选址,规划路网结构,规划公共交通站点,等等。

政苑社区的等时圈,三个颜色层级分别为步行 5 分钟、10 分钟、15 分钟。
政苑社区的等时圈,三个颜色层级分别为步行 5 分钟、10 分钟、15 分钟。

大三给政苑社区做调研分析的时候,我就绘制过该社区出发步行 15 分钟的等时圈。可以看出,这个等时圈实际不是标准的「圆圈」,而是沿着道路伸出畸形触角的图形。可以说,这「畸形」就是等时圈分析的意义所在,因为它比以某个点出发画一个半径(比如 1.5 公里)的圆更加真实地反映了交通路径的实际情况。

尤其是,比如小区边上有一条大河或者一条快速路,大河或马路对面就不那么好过去了。可能直线距离离你只有 100 米,但为了过去,得花上二十多分钟。

下面,让我来介绍等时图的几种画法。

构建路网计算法:使用网络分析工具手动计算

有一种原始的美,或许很久很久以前,人们就是这么做的。

或可参考这篇文章,我不细致介绍了。

现成法:使用现成网页分析工具

21 世纪 20 年代了,我们坚信总是有好人会写好现成的工具的。确实。城规在这点上是幸运的,还有一些地图厂商,以及想赚钱的企业,做了一点工作为我们提供帮助。

首先要感谢 Mapbox 提供 Isochrone API,这是大部分目前在线版工具的基础。上面的等时图就是我使用 微思 Layers ,调用 Mapbox 的 API 做的。

可以在 Mapbox 的 Isochrone API Playground 中轻松绘制等时圈
可以在 Mapbox 的 Isochrone API Playground 中轻松绘制等时圈

类似的工具我列在此处,希望所有网站都能一直活下去:

其中,有几个工具也提供了腾讯地图 API 调用或高德地图 API 调用,后者可以做所谓的公共交通出行(公交+地铁)等时圈。(注:事实上,百度地图也提供等时圈 API,但是并不给普通开发者免费使用,需要额外联系。)

是啊,那么简单,只要点一个起点,然后选一选时间和出行方式就好了。五分钟完成出图!

一般的读者(比如城市研究爱好者,对我的杂文感兴趣的小可爱🥰)看到这里已经可以结束了。你已经知道了最重要的概念和做法,已经可以更好地使用工具认识城市了!

但是对于更严苛更严谨、思考更加深刻的人来说……这件事真的这么简单吗?代价呢?

  1. Mapbox 的路网数据主要来源是 OSM(OpenStreetMap),众所周知,它的路网数据在国内往往并不准确,会有缺失。
  2. Mapbox 根据路网的不同「等级」(来自 OSM 的道路标注)来生成的所谓的道路通行速度,然而这个速度并不公开,而且也并非实际道路限速。(注:Mapbox 最新的 API 提供了「driving-traffic」的方法,然而它用的交通数据对于国内……好像也不准确。)
  3. 如果,你做的研究,还想考虑交通拥堵,某个特定时刻的交通状况……那它也无法满足。

总结起来,就是这个交通时间,还是不够「真实」。

可是……怎么得到最真实准确的交通时间呢?难道真的要开车/走路,走一遍吗?

算路法:使用导航提供商的 API 批量算路

你想到,我们在用百度(高德)地图导航的时候,选一个目的地,总会得到一个「预计时间」,诶,这个时间,其实就是相对准确的真实交通时间。尤其是,它真的可以考虑到交通状况、限速、红绿灯……

那么,我就可以在研究区域内选取一堆终点,记录它们的交通时间,最后把 x 分钟的所有点连起来,就能画出等时圈了!

是的,百度地图和高德地图也分别提供了路径规划(用于导航方案)的 API(百度的高德的),但是他们返回的都是非常完整的导航方案,固然包含交通时间数据,但是有点大材小用。而且这个请求返回耗时久,百度还把它定为高级功能,个人用户也无法使用。

但好消息是,百度提供了这个计算时间的 API:批量算路,根据起点和终点坐标计算路线规划距离和行驶时间。嗯,这本意并不是给我们规划研究者提供的,其实是「适用于高并发场景,如网约车派单、物流配送派单场景,同时发起多个起终点之间的算路,筛选所需要的订单起终点」。当然,既然能用,我们也拿来用。这样可以自动得得到交通时间(而不用手动在导航软件中测试 qwq)

有了这个 API,虽然不如 Mapbox 方便,但也有很多手段来具体操作。

以下,我使用 Python + arcpy 库(需要 ArcGIS Pro),通过调用百度批量算路 API 的方式,来示范绘制一个真实的等时圈。

提前的 Acknowledgement:
后续大部分代码来源于该文章: 利用 ArcGIS_Python 制作考虑路况的交通等时圈。非常感谢原作者 Renhai实验室。我针对不规则研究区域和大批量终点的多次调用改进重写了代码,增加了一些针对初学者的解释。
另外,我后续了解到可以 QGIS 里或者有其他的插件直接在 GIS 分析软件内执行 http 请求的调用,可以参考文末的【参考阅读3】视频。我认为使用 QGIS 可能是一个更加简明的方式。
或者,如果读者知晓 FME(一个空间数据转化处理软件),也可以考虑用它来实现 GIS+空间数据处理。

提前的准备

  • 基本了解 Python,以及虚拟环境、conda、pip 的相关概念;
  • Windows 电脑,安装好 ArcGIS Pro;
  • 克隆 ArcGIS 的默认 python 环境 arcgispro-py3,并在你的 IDE 和终端中切换到你克隆的新环境。使用这个新环境作为内核激活 notebook (另外,我推荐使用 vscode)。

绘制研究区域

先在 ArcGIS 中绘制面要素和点要素也是可以的,但是要先创建要素啊,编辑绘制啊,比较麻烦,这里出图之前,直接全程使用 python 工具完成编辑和绘制。

首先确定出发点和框定一个研究区域(区域内的所有点作为出行终点)。这一步需要在 notebook 中完成,需要交互式的环境。(或者,先在 ArcGIS 中画好,后续再通过 python 导入。)

1
2
3
4
5
6
7
8
9
10
11
# 可能需要先在终端中执行 pip install leafmap
import leafmap

zoom = 14
preview_map = leafmap.Map(
center=[30.252714, 120.20979],
zoom=zoom,
draw_control=True,
)

preview_map

执行应当会打开一个地图 widget。选取地图左侧的标注工具,使用 Marker 工具,点出出发点;使用 Polyline 工具(矩形工具也行),绘制需要的研究区域。

执行这个 cell 会打开 leafmap 的 widget,在这里画范围比较快捷方便
执行这个 cell 会打开 leafmap 的 widget,在这里画范围比较快捷方便

绘制多边形,而不是一个简单的矩形,是因为我个人需要绘制最远60分钟左右的驾驶等时圈,而这,会产生大量的终点;但是这个百度的 API 是限额 5000 次/天的,为了节省费用,这里预先测试一下大概的等时圈边界,只把研究范围缩减到一个较精确的范围。另外,考虑到后续的研究实际需求,我也删去了钱塘江另一侧的区域。

我没有截我本来在这个 widget 中绘制时的图,不过实际在 ArcGIS 中,研究范围和起点大概是这样的
我没有截我本来在这个 widget 中绘制时的图,不过实际在 ArcGIS 中,研究范围和起点大概是这样的

输出并记录一下 points 点位的坐标和 polygon 所有顶点的坐标。

1
2
3
4
5
6
7
8
9
10
11
12
data = preview_map.draw_control.data

points = []
polygon = None

for i in data:
if i['geometry']['type'] == 'Point':
points.append(i['geometry']['coordinates'])
print("points:", points)
if i['geometry']['type'] == 'Polygon':
polygon = i['geometry']['coordinates'][0]
print("polygon:", polygon)

由于 notebook 重启后所有中间变量数据会丢失,我不想重新绘制研究区域,我会手动记下这个输出。

1
2
3
4
5
points = [[120.20979, 30.252714]]  # 根据上一个 cell 的实际输出修改
polygon = [[120.139618, 30.19796], [120.115585, 30.255509], [120.109406, 30.294054], [120.110779, 30.325471], [120.247421, 30.408413], [120.324669, 30.401603], [120.397453, 30.351843], [120.396423, 30.294646], [120.363464, 30.259067], [120.337029, 30.264701], [120.314369, 30.286642], [120.285187, 30.298796], [120.238495, 30.28012], [120.201073, 30.229111], [120.157814, 30.203894], [120.139618, 30.19796]] # 根据上一个 cell 的实际输出修改

print("Final points:", points)
print("Final polygon:", polygon)

确认一下输出是对的。以后再次打开激活 notebook 时,可以从这个 cell 开始运行。

使用 arcpy 库批量生成终点

先在当前的工作文件夹(notebook 或 python 文件所处的文件夹)下新建一个 resource 文件夹。下面的代码中,我们会在里面新建(如果已有,则是打开)一个叫 data.gdb 的 地理数据库文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 导入arcpy
import arcpy

# 设置工作空间
home_dir = os.path.join(os.getcwd(), "resource")

# 创建一个地理数据库
if not arcpy.Exists(os.path.join(home_dir, "data.gdb")):
arcpy.CreateFileGDB_management(home_dir, "data.gdb")

arcpy.env.workspace = os.path.join(home_dir, "data.gdb")
arcpy.env.overwriteOutput = True

# 设置空间参考对象
sr = arcpy.SpatialReference("WGS 1984") # 百度api可以直接使用WGS84坐标
# 设置输出坐标系
arcpy.env.outputCoordinateSystem = sr

根据 polygon 的顶点坐标,在 data.gdb 里生成叫做 study_area 的要素。

1
2
3
4
5
6
7
8
9
10
# 把 polygon 转化为 arcgis 的几何对象
coords = polygon if polygon[0] == polygon[-1] else polygon + [polygon[0]]
ring = arcpy.Array([arcpy.Point(x, y) for x, y in coords])
study_area = arcpy.Polygon(ring, sr)

# 将其转换为要素类并复制到数据库
arcpy.management.CopyFeatures(study_area, "study_area")

# 检查一下 study_area 形状对不对
study_area
notebook 中大概是这样一个输出,小小的也很可爱
notebook 中大概是这样一个输出,小小的也很可爱

下一步生成渔网。为了让终点覆盖到最靠近研究范围边界的地方,生成渔网时,范围要比研究范围稍微大些。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 计算多边形的边界框
# 提取所有经度和纬度
lons = [point[0] for point in polygon]
lats = [point[1] for point in polygon]

# 计算边界框
x_min = min(lons)
x_max = max(lons)
y_min = min(lats)
y_max = max(lats)

# 偏移量
offset = 0.001

# extent用来指定渔网的范围
extent = " ".join([str(coord) for coord in [x_min - offset, y_min - offset, x_max + offset, y_max + offset]])

out_fcs = "study_area_fishnet"
# 避免输出文件已存在
if arcpy.Exists(out_fcs):
# 生成唯一的输出文件名
unique_name = arcpy.CreateUniqueName(out_fcs)
out_fcs = os.path.basename(unique_name)
print(f"输出要素类已存在,已经重命名为: {out_fcs}")

下面的 cell_width 会控制渔网的大小,也就是改变生成终点的数量多少。可以多次修改这个值,将终点的数量调整到合适,满足自己 API 用量的需求。

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
# 修改间距大小
cell_width = cell_height = 0.0028 # 0.001度约为100米

# 创建渔网
out_fishnet = arcpy.management.CreateFishnet(out_feature_class = out_fcs, # 包含由矩形像元组成的渔网的输出要素类(的名称)。
origin_coord = str(x_min-offset) + " " + str(y_min-offset), # 矩形框的左下端点为原点
y_axis_coord = str(x_min-offset) + " " + str(y_max-offset), # 此点与原点的连线用于判断旋转的角度 我们不用旋转所以定义为原点正上方的点
cell_width = cell_width,
cell_height = cell_height,
number_rows = "", # 留空,由cell_width和cell_height决定
number_columns = "", # 留空,由cell_width和cell_height决定
corner_coord = None, # 对角坐标不填写
labels = "LABELS", # 输出标签 即输出渔网中心点的要素类
template = extent, # 以空格分隔的坐标字符串 - 将使用指定渔网的范围。坐标以 x-min,y-min,x-max,y-max 的顺序表示。
geometry_type = "POLYGON" # 生成面
)

# 定义新产生的点要素的名称
out_label = out_fcs + "_label"

# 用研究区面裁剪 label 点,只保留研究区内部的格心。
# 如果是矩形的研究区域,这一步其实可以省略。
# 此时应当赋值 out_label_clip = out_label
out_label_clip = out_label + "_clip"

arcpy.analysis.Clip(in_features = out_label,
clip_features = study_area,
out_feature_class = out_label_clip)

检查一下要素的信息(主要是看数量,如果过多或过少则返回来修改上一个 cell 中的cell_width):

1
2
3
4
5
6
desc = arcpy.Describe(out_label_clip)
print(f"要素类名称: {desc.name}")
print(f"要素数量: {arcpy.management.GetCount(out_label_clip)[0]}")
print("\n字段列表:")
for field in desc.fields:
print(f" - {field.name} ({field.type})")

输出应该类似:

1
2
3
4
5
6
要素类名称: study_area_fishnet_label_clip
要素数量: 4705

字段列表:
- OBJECTID (OID)
- Shape (Geometry)

一直到这一步,如果读者对 ArcGIS Pro 的操作更加熟悉,其实也是可以直接在软件中完成的。可以后续再在 python 中读取 out_label_clip,进行下一步 API 的调用。

构造请求 url

百度批量算路 API 存在一些限制,但对于规划的城市研究情景来说,不是很重要。具体可以参考它的开发文档

在我们的代码中,构造的请求 url 形式如:

1
2
url = f'https://api.map.baidu.com/routematrix/v2/{trans_type}?
output=json&origins={origins}&destinations={des}&tactics={tactics}&coord_type=wgs84&ak={ak}'

其中,trans_type 可选 drivingwalkingriding(注:riding 可以在 url 请求参数中进一步区分 riding_type 是自行车还是电动车,但摩托车需要向百度申请权限)tactics 对于驾车情景,是策略选择,具体的说明为:

1
2
3
4
10:不走高速;  
11:常规路线,即多数用户常走的一条经验路线,满足大多数场景需求,是较推荐的一个策略;
12:距离较短(考虑路况):即距离相对较短的一条路线,但并不一定是一条优质路线。计算耗时时,考虑路况对耗时的影响;
13:距离较短(不考虑路况):路线同以上,但计算耗时时,不考虑路况对耗时的影响,可理解为在路况完全通畅时预计耗时。

百度注明,除 13 外,其他策略的耗时计算都考虑实时路况。也就是说,在一天的早上、中午、晚上不同时间调用 API,由于每次都是按请求时的路况来计算,做出来的等时圈结果也是不一样的。

我们这次试试推荐策略 11,我是在下午 3 点左右运行的代码,此时应当并不会因为早晚高峰改变结果的普适性。

ak 则是你的 access key,可以前往百度地图开放平台的应用管理页面,创建一个应用,得到 ak。代码中不建议明文存储这类访问密钥,我们在工作文件夹下新建一个 .env 文件,写入 baidu_ak = "xxxxxx……(你的ak)",保存。

让 python 读入这个 ak。顺便,定一下我们的交通工具和策略吧。

1
2
3
4
5
6
7
8
from dotenv import load_dotenv
import os

load_dotenv(".env") # 读取环境变量文件
ak = os.getenv("baidu_ak") # 读取百度api

trans_type = 'driving' # 驾车出行
tactics = '11' # 采用常规路线的策略

然后,我们根据 out_label_clip 中每个点的坐标,批量构造我们待会儿要请求的 url,把它们保存到字典里:

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
import urllib.parse
import arcpy

data_dict = {}
des = f"{points[0][1]},{points[0][0]}"
# out_label_clip = "study_area_fishnet_label_clip"

# 读取点坐标
with arcpy.da.SearchCursor(out_label_clip, ["OBJECTID", "SHAPE@XY"]) as cursor:
for oid, xy in cursor:
lat, lon = xy[1], xy[0]
coord_str = f"{lat},{lon}"

data_dict[oid] = {"coord": coord_str}

# 为每个点构造 URL
base_url = f"https://api.map.baidu.com/routematrix/v2/{trans_type}"

for oid, info in data_dict.items():
params = {
"output": "json",
"origins": info["coord"],
"destinations": des,
"tactics": tactics,
"coord_type": "wgs84",
"ak": ak
}
query_string = urllib.parse.urlencode(params)
info["url"] = f"{base_url}?{query_string}"

print(f"构造了 URL 数量: {len(data_dict)}")

批量调用 API

我们可以把调用得到的通行时间存进一张表格,这样方便在出错的时候也保留一部分数据,或者,分成好几天来请求百度给我数据。

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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
import asyncio
import aiohttp
import json
import csv
import os

# ---------------------------
# 工具函数
# ---------------------------

async def perform_request(url):
"""发送请求并返回文本内容;失败则返回 response.text 供错误输出"""
async with aiohttp.ClientSession() as session:
try:
async with session.get(url) as response:
text = await response.text()
if response.status == 200:
return text
else:
return {"error": True, "msg": text}
except Exception as e:
return {"error": True, "msg": str(e)}


async def limited_fetch(url, semaphore):
"""限制并发"""
async with semaphore:
result = await perform_request(url)
await asyncio.sleep(1)
return result


def get_time(content):
"""从 JSON 内容解析 duration(分钟)"""
try:
item = json.loads(content)
if item.get("status") == 0:
duration = item['result'][0]['duration']['value'] / 60
return round(duration, 1)
else:
# 不打印重复消息,只返回 service 的 message
return None
except:
return None


# ---------------------------
# 并发请求 -> 写入CSV(断点续跑 + 覆写更新)
# ---------------------------

async def main_fetch_to_csv(csv_path, data_dict):
"""并发请求,写CSV,支持断点续跑"""

# ---- 读取已有 CSV ----
csv_map = {}
order_oids = []

if os.path.exists(csv_path):
with open(csv_path, "r", encoding="utf-8") as f:
reader = csv.DictReader(f)
for row in reader:
oid_str = (row.get("oid") or "").strip()
dur_str = (row.get("duration_min") or "").strip()
if not oid_str:
continue
oid = int(oid_str)
if oid not in csv_map:
order_oids.append(oid)
csv_map[oid] = dur_str

def has_value(v):
return v not in ("", None, "None")

done_oids = {oid for oid, v in csv_map.items() if has_value(v)}
print(f"已完成:{len(done_oids)} / {len(data_dict)}")

# ---- 确保 CSV 存在表头 ----
if not os.path.exists(csv_path):
with open(csv_path, "w", newline='', encoding="utf-8") as f:
writer = csv.writer(f)
writer.writerow(["oid", "duration_min"])

semaphore = asyncio.Semaphore(10)
csv_lock = asyncio.Lock()

# ---- CSV 覆写函数 ----
def flush_csv():
with open(csv_path, "w", newline='', encoding="utf-8") as f:
writer = csv.writer(f)
writer.writerow(["oid", "duration_min"])
for oid in order_oids:
writer.writerow([oid, csv_map.get(oid, "")])
# 避免遗漏新增 oid
for oid in csv_map:
if oid not in order_oids:
writer.writerow([oid, csv_map[oid]])

# ---- 单个请求任务 ----
async def request_task(oid, url):
if oid in done_oids:
return

result = await limited_fetch(url, semaphore)

if isinstance(result, dict) and result.get("error"):
print(f"[请求失败] oid={oid}{result['msg']}")
duration_min = None
else:
duration_min = get_time(result)

async with csv_lock:
if oid not in csv_map:
order_oids.append(oid)
csv_map[oid] = "" if duration_min is None else str(duration_min)
flush_csv()

print(f"[写入] oid={oid}, duration_min={duration_min}")

# ---- 创建任务 ----
tasks = []
for oid, info in data_dict.items():
if oid not in done_oids:
tasks.append(asyncio.create_task(request_task(oid, info["url"])))

if not tasks:
print("没有需要处理的 oid。")
return

await asyncio.gather(*tasks)
print("CSV 更新完成")


# ---------------------------
# 调用
# ---------------------------

# 在当前文件夹下创建 / 使用 durations.csv
csv_path = "durations.csv"
await main_fetch_to_csv(csv_path, data_dict)

输出可能类似

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
已完成:0 / 4705
[写入] oid=6, duration_min=29.7
[写入] oid=5, duration_min=20.8
[写入] oid=8, duration_min=30.0
[写入] oid=10, duration_min=21.6
[写入] oid=2, duration_min=28.2
[写入] oid=1, duration_min=28.2
[写入] oid=7, duration_min=30.3
[写入] oid=9, duration_min=27.4
[写入] oid=3, duration_min=22.4
[写入] oid=4, duration_min=22.2

......(省略)

CSV 更新完成

你会注意到,oid 的更新并不是按序的。这是因为请求的返回时间并不一致,这并不影响结果。

这段代码是支持当再次调用时,继续从未得到结果的点位开始继续请求的。也就是说,如果你的研究区域比较大,突破了百度一天 5000 次请求的限制,可以分成好几天来运行……(我建议每天在同一个时刻运行)

另外百度可能会发邮件或短信提示你「使用的批量算路服务并发量已接近约定上限」,同时,有些 url 请求结果因此被限制而报错。前者不用理会,后者的话,这段代码结束运行之后,再多运行几次,会对没有得到结果的点位继续请求,最后总会得到完整的、包含所有终点的时间的结果的表格的。

使用请求结果更新 ArcGIS

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
def update_row(out_feature_class, oid, duration_min):
"""更新游标,将duration_min填入time字段"""
field_names = ["OBJECTID","SHAPE@XY", "time"] # 定义字段名称
where_clause= f"OBJECTID = {oid}"

with arcpy.da.UpdateCursor(out_feature_class, field_names, where_clause) as cursor:
for row in cursor:
row[2] = duration_min # 将时间填入time字段
cursor.updateRow(row)
print(f"成功更新{oid}{duration_min}分钟")

def main_update_from_csv(csv_path):
fcs = out_label_clip
ws = arcpy.env.workspace

if not arcpy.ListFields(fcs, "time"):
arcpy.management.AddField(fcs, "time", "DOUBLE")

# 开启编辑并逐行更新
with arcpy.da.Editor(ws) as edit: # 对工作空间开启编辑,防止锁占用
with open(csv_path, "r", encoding="utf-8") as f:
reader = csv.DictReader(f)
for row in reader:
oid = int(row["oid"]) if row["oid"] else None
duration_min = float(row["duration_min"]) if row["duration_min"] not in (None, "", "None") else None
if oid is not None and duration_min is not None:
update_row(fcs, oid, duration_min)

main_update_from_csv(csv_path)

绘制等时圈

我们得到的是许许多多个点的时间,然后,我们要进行插值,得到面状的时间分布(好抽象,其实意思是,我们要为没有测量值的点来「预测」值,这个预测是使用周边测量点的已有值来预测的)。

插值的算法还蛮多的,我们就用反距离权重法吧。

1
2
3
4
5
6
7
8
Idw_raster = arcpy.sa.Idw(
in_point_features=out_label_clip,
z_field="time",
cell_size=7.99999999999272E-05,
power=2,
search_radius="VARIABLE 12",
in_barrier_polyline_features=None
)

我们会得到一张漂亮的栅格。我们把它裁剪到研究区域。

1
2
3
4
5
6
7
8
arcpy.management.Clip(
in_raster=Idw_raster,
out_raster="Idw_raster_clip",
in_template_dataset="study_area",
nodata_value="3.4e+38",
clipping_geometry="NONE",
maintain_clipping_extent="NO_MAINTAIN_EXTENT"
)

大功告成!现在,我们可以打开 ArcGIS Pro。找到地理数据库中剪裁后的栅格,把它拖到地图上,然后在符号系统中调整一下不同交通时间段的颜色,就可以看到漂亮的等时图了~

后续,我们就可以在等时圈的基础上做进一步的研究。

结语

大部分规划同学的代码基础真的太弱了……写这篇文章是因为,老师同时让我们好几个师门同学来画真实的等时圈,结果都不知道怎么画,跑来问我。我做了一个分步阐述做法和成果的 PPT,并把我的 notebook 给了他们,然而……最后他们依然画不出来。

当然一部分原因是 ArcGIS 的 arcpy 库太恶心了,并不公开提供(商业软件可以理解),然后 ArcGIS 自己又集成了一套 conda 工具,包管理与调用异常麻烦,和主流 python 世界完全脱轨。

但是最重要的……我想大家找不到手把手的用于规划的分析教程才是主要原因。或者,能找到的都还不够手把手(?)

希望我的文章能教会读者叭。不管是不是专业人士,都能看懂就最好了。如果看完还不会,可以留言或者邮箱联系我😭


参考阅读与参考资料:

  1. 位置分析那些事|如何从0到1进行等时圈分析?
  2. FME 与等时圈
  3. QGIS系列视频(八):使用百度API,利用QIGS计算15分钟步行等时圈,并进行可视化和提取等时线
  4. 利用 ArcGIS_Python 制作考虑路况的交通等时圈
  • 觉得本文不错?
    本文采用 CC BY-NC-SA 4.0 协议共享,欢迎分享与转载。