参考文章：城市公交可达圈绘制方法（一） - 知乎 (zhihu.com)
本篇文章我们聚焦于通过公共交通出行方式（包括公交、地铁、公交+地铁的组合）来获取一定时间内可以到达的范围。为了实现这一目标，我们将使用高德地图API中的公交到达圈功能，对城市某一点的公交可达圈进行详细分析。通过这一分析，我们可以更好地理解城市的公共交通网络覆盖范围，为城市规划、交通优化和个人出行提供有力的数据支持。我们将利用高德地图API提供的强大功能，结合具体的起始点和时间参数，获取并可视化从该点出发，在指定时间内可以通过公交和地铁到达的所有区域。
先讲一下方法思路，一共三个步骤；
方法思路
- 确认URL位置和数据结构，并根据配置参数构建请求URL
- 坐标转换——高德坐标系（GCJ-02） to WGS84
- 点集转成csv
老规矩，通过快捷键 Ctrl + Shift + i，找到数据源的位置，然后打开这个url的链接，可以看到数据源是一个个点坐标，我们需要把点坐标提取出来，并把点的高德坐标系（GCJ-02）转为WGS84坐标；
公交到达圈-公交信息查询-示例中心-JS API 1.4 示例 | 高德地图API (amap.com)
参数设置：location（坐标）、key（高德key）、time（出行时间）、strategy（出行方式）；
config = {'key': '你的key', # 高德的API的key'location': '118.797539,31.969705', # 坐标可以自定义'time': '30', # 时间选择范围有(1-45)都可以输入,通常我们使用(15、30、45)'strategy': '1' # 如果strategy==2,strategy_name='公交+地铁';如果strategy==1,strategy_name='地铁';如果strategy==0,strategy_name='公交';}
这些参数要提前自定义好，这里以上海虹桥火车站为例121.320011,31.19403（GCJ-02），time参数：45（分钟），strategy：2（地铁+公交，'0': '公交','1': '地铁', '2': '公交+地铁'），key的话在高德地图开放平台上注册一个账号并创建一个应用，获取API密钥即可；
完整代码#运行环境Python 3.11
import requests
import json
import logging
import math
import csv# 配置日志
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')#坐标转化函数
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626 # π
a = 6378245.0 # 长半轴
ee = 0.00669342162296594323 # 偏心率平方def gcj02_to_wgs84(lng, lat):"""将GCJ02坐标系（火星坐标系）转换为WGS84坐标系:param lng: 火星坐标系的经度:param lat: 火星坐标系的纬度:return: WGS84坐标系的经纬度列表"""if out_of_china(lng, lat):return [lng, lat]dlat = _transformlat(lng - 105.0, lat - 35.0)dlng = _transformlng(lng - 105.0, lat - 35.0)radlat = lat / 180.0 * pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)mglat = lat + dlatmglng = lng + dlngreturn [lng * 2 - mglng, lat * 2 - mglat]def _transformlat(lng, lat):"""辅助函数，用于计算纬度偏移:param lng: 经度:param lat: 纬度:return: 纬度偏移值"""ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *math.sin(2.0 * lng * pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lat * pi) + 40.0 *math.sin(lat / 3.0 * pi)) * 2.0 / 3.0ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *math.sin(lat * pi / 30.0)) * 2.0 / 3.0return retdef _transformlng(lng, lat):"""辅助函数，用于计算经度偏移:param lng: 经度:param lat: 纬度:return: 经度偏移值"""ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *math.sin(2.0 * lng * pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lng * pi) + 40.0 *math.sin(lng / 3.0 * pi)) * 2.0 / 3.0ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *math.sin(lng / 30.0 * pi)) * 2.0 / 3.0return retdef out_of_china(lng, lat):"""判断给定的经纬度是否在中国范围内:param lng: 经度:param lat: 纬度:return: 如果在中国范围内返回False，否则返回True"""return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)def lines_wgs84(coords):"""将一系列GCJ02坐标转换为WGS84坐标:param coords: GCJ02坐标列表，每个坐标为 "经度,纬度" 格式的字符串:return: WGS84坐标列表"""wgs84_coords = []for coord in coords:lng, lat = map(float, coord.split(','))wgs84_coords.append(gcj02_to_wgs84(lng, lat))return wgs84_coordsdef build_url(config):"""构建请求URL:param config: 配置字典:return: 完整的请求URL"""base_url = 'https://lbs.amap.com/AMapService/v3/direction/reachcircle'params = {'key': config['key'],'location': config['location'],'time': config['time'],'s': 'rsv3','extensions': 'all','output': 'json','strategy': config['strategy'],'callback': 'jsonp_89634_','platform': 'JS','logversion': '2.0','appname': 'https%3A%2F%2Flbs.amap.com%2Fdemo%2Fjavascript-api%2Fexample%2Fbus-info%2Farrival-range','csid': 'B6FE1E8C-C71D-4EDB-9261-2CBCD8A0C433','sdkversion': '1.4.27'}query_string = '&'.join([f'{k}={v}' for k, v in params.items()])return f"{base_url}?{query_string}"def main(config):# 构建URLurl = build_url(config)# 设置请求头headers = {'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/129.0.0.0 Safari/537.36 Edg/129.0.0.0','Referer': 'https://lbs.amap.com/demo/javascript-api/example/bus-info/arrival-range','Cookie': 'AMAPID=d365685e5855236e1224b874d3d5528f; passport_login=MTcyMzA0NzUzLGFtYXBBVGo4T0h1R1IsZW9odmR2bTJhY21iY2FjaXpjNTZkZmlqcGNobXR0Y3ksMTcyNjcxNzM2MyxNbUZtTVdVek5qWmtNRFV6T0dVNFpEUmxZVFZtTWpFNE1EYzFZalUxTXpVPQ%3D%3D; gray_auth=2; isg=BH9_D-UqJLpNOyOjsUs-Ceo0DlMJZNMGzI2-3hFNHS51IJuiGTUrVPIyZnDeY6t-; tfstk=gnDqqA64FKp4UqvQVueaY0RZSQeYd-YBSAa_jcmgcr4m5CUiaVgyDPiMkVkar2cbG-wD_RoE7tMb1jigjzga5mwDDG2ZPc2sG-AOIP4Su5G_h-2ascuV1nHZXOrijVKY5ndSDmeTIeTI7pixDQ6wUVM4muVo2kzgjpNsmxCgIeTBFaITzkyilRzFQ_ouylrcSt0gZ_48ft2go5blZk4uSR0goaquYk7GIS40E3rTrP2gobrmmcXzASxCb8QkRQd8iym0zOXEvoPD4m11BOC3DSaHpz6gVYr4gymmHhnLcuDSUWN5fUwovbg04JJVJoloY-qZBHfa75kYUzueKOE-iqku_Yt9K0w4bW2gaGXi48z8axwHKar-ZmGiHqSNIoHjd5zLahXT18DQtX0VXHnunl0T9ATdHrmivvhQLK7L0bmUUgScWuVmRAhVS1V02uzB43S6dcPx6AvtT1CTZJEzRnZf61F02uzB435O67qL4ytbc'}# 发送HTTP请求获取数据try:response = requests.get(url, headers=headers)response.raise_for_status()logging.info("请求成功")# 解析JSONP响应data = response.text.strip('jsonp_89634_(').strip(')')data = json.loads(data)# 打印数据logging.info(json.dumps(data, indent=4, ensure_ascii=False))# 处理多边形数据coordinates = []for i, polygon in enumerate(data['polylines']):outer_coords = polygon['outer'].split(';')outer_points = lines_wgs84(outer_coords)for point in outer_points:coordinates.append([f"多边形 {i + 1} 外轮廓", point[0], point[1]])logging.info(f"多边形 {i + 1} 的外轮廓点 (WGS84): {outer_points}")if polygon['inners']:for j, inner_coords in enumerate(polygon['inners']):inner_points = lines_wgs84(inner_coords.split(';'))for point in inner_points:coordinates.append([f"多边形 {i + 1} 内轮廓 {j + 1}", point[0], point[1]])logging.info(f"多边形 {i + 1} 的内轮廓点 {j + 1} (WGS84): {inner_points}")# 导出CSV文件with open('reach_circle.csv', mode='w', newline='', encoding='utf-8-sig') as file:writer = csv.writer(file)writer.writerow(['多边形名称', 'X坐标', 'Y坐标'])writer.writerows(coordinates)logging.info("CSV文件已保存为 reach_circle.csv")except requests.RequestException as e:logging.error(f"请求失败: {e}")except json.JSONDecodeError as e:logging.error(f"JSON解析失败: {e}")if __name__ == "__main__":# 用户自定义参数config = {'key': '你的key', # 高德的API的key'location': '121.320011,31.19403', # 坐标可以自定义'time': '45', # 时间选择范围有(1-45)都可以输入,通常我们使用(15、30、45)'strategy': '2' # 如果strategy==2,strategy_name='公交+地铁';如果strategy==1,strategy_name='地铁';如果strategy==0,strategy_name='公交';}main(config)
脚本运行结束，就会生成一个名为 reach_circle.csv 的CSV文件，其中包含多边形的名称、X坐标和Y坐标，并且文件编码为UTF-8；
我们把结果添加到GIS里做可视化，先添加xy坐标，即【显示xy数据】，然后在工具里检索【点集转线】，把点集转成线层；
最后在ArcToolbox，【数据管理工具】→【要素】→【要素转面】，注意：线要素中的线必须是闭合的；
这样我们就得到了从虹桥火车站出发，45分钟内地铁+公交的可达范围。由图可知，我们可以清晰地看到虹桥火车站周边的公共交通网络布局以及各个方向上的可达性情况。
- 向西最远可以到达徐泾新区，这一区域包括了多个住宅区和商业综合体。
- 向东最远可以到达世纪大道，世纪大道是浦东新区的核心地带，汇集了众多高端商务楼宇和购物中心。
- 向北可以到达七宝老街，七宝老街坐拥悠久的历史文化和独特的江南水乡风貌，可以满足旅客45分钟直达。
- 向南可以到达丰庄，丰庄可是上海市普陀区的一个重要节点社区，拥有丰富的居民生活设施和便利的交通条件。
文章仅用于分享个人学习成果与个人存档之用，分享知识，如有侵权，请联系作者进行删除。所有信息均基于作者的个人理解和经验，不代表任何官方立场或权威解读。