Python批量实现MapBox等时圈

一、说明

最近做项目用到了等时圈，出了很多张图，准备记录一下。

做等时圈有很多种办法，一是用高德地图出行时间的API接口，原理大概是通过产生等距离栅格渔网点，计算目标点到周边渔网点的距离或时间，最后利用ArcGIS中的连接，将获得的数据和之前产生的渔网点关联起来，再以距离或时间为区分字段调整显示效果；或者将渔网点进行插值，将会产生更平滑的效果，但最终生成的等时圈都是栅格图，一是效果不好，二是计算量特别大，根据自己的精度要求和需要计算的点的数量，其需要调用API的次数是很多的，特别是高德现在更改了接口调用的限额，所以导致做起来比较麻烦。

还有一种方法是直接用ArcGIS内置的等时圈计算，需要用到路网数据，而且需要对路网数据进行处理，拓扑检查、修复等，还是比较麻烦，但是优点就是路网可以自己进行修补调整，路网精度越高，最终结果越准确，特别是针对一些小范围地块，高德、百度、OpenStreetMap等地图商没有地块内部道路数据，在这种情况下，得出的结果会更准确。

前面两种方式都有自己的优劣势，利用MapBox等时圈API接口来直接计算虽然或多或少有一定缺点，但毕竟简单省事！（MapBox也很大方，提供接口的额度都比较大，不像是某德 ）

二、大概出图逻辑

由于没有提供POI点，所有数据都是以图片形式给我的，所以需要先在ArcGIS中进行地理配准（好在给我的图片带有地图，有路网信息，可以通过对照路网来进行地理配准），配准结束后需要手动矢量化信息点（好累），矢量化后利用ArcGIS中计算几何工具计算点的经纬度（注意MapBox需要WGS84坐标系），最后利用表到Excel工具将生成的点的经纬度输出到Excel表中（记得勾选....作为列标题，忘记叫什么了，反正是第一个选项，后续有空再补图），数据处理内容大致结束。

转到MapBox中，使用代码获取等时圈shp数据，大概就完成啦！

三、代码

借鉴了Python|批量获取Mapbox等时圈shapefile格式数据 - Weltㅤ - 博客园，在此基础上进行了一定程度的修改。

import time
import requests
import pandas as pd
import shapefile
from osgeo import osr
import openpyxl
import os

# contours_minutes
# 文件读取路径
path = r'C:\xxxx\小型体育场_TableToExcel.xlsx'


def getdata(file_name, OBJECTID, lon, lat):
    print('获取{}的等时圈'.format(OBJECTID))
    # 注意根据需要修改url中运动时间、计算方式（距离、时间）、出行形式
    url = 'https://api.mapbox.com/isochrone/v1/mapbox/walking/{},{}?contours_minutes=5&polygons=true&denoise=1&access_token=你的token'.format(lon, lat)
    try:
        r = requests.get(url, timeout=(3, 7)).json()
        #  修改文件输出路径
        w = shapefile.Writer('.\{}_{}.shp'.format(file_name, OBJECTID), encoding='gbk')
        w.field('name', 'C')
        w.field('transit', 'C')
        w.field('contour', 'N')
        w.field('metric', 'C')
        w.field('color', 'C')
        w.field('opacity', 'N', decimal=2)

for i in r['features']:
            coords = i['geometry']['coordinates']
            w.poly(coords)
            w.record(
                name=OBJECTID,
                contour=i['properties']['contour'],
                metric=i['properties']['metric'],
                color=i['properties']['color'],
                opacity=i['properties']['opacity']
            )
        w.close()
        # 添加投影信息
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(4326)
        wkt = proj.ExportToWkt()
        projfile = '.\{}_{}.prj'.format(file_name,OBJECTID)
        f0 = open(projfile, 'w')
        f0.write(wkt)
        f0.close()
        print('获取成功！')
    except Exception as e:
        print(e)

if __name__ == '__main__':
    print('开始爬取数据')
    file_name = os.path.basename(path)[:-17]
    wb = openpyxl.load_workbook(path)
    # 获取当前所有的sheet
    sheets = wb.worksheets
    # 读取第一个sheet表格
    sheet1 = sheets[0]
    # 获取行数
    max_row_num = sheet1.max_row
    max_col_num = sheet1.max_column
    for i in range(2, max_row_num + 1):
        value_list = []
        for j in range(1, max_col_num + 1):
            cell = sheet1.cell(i, j).value
            value_list.append(cell)
        getdata(file_name, value_list[0], value_list[1], value_list[2])
        time.sleep(0.5)
    print('All Done!!!')