基于Python的地理空间分析（八）：Mapbox与Python

发表于 分类于 阅读次数： Disqus：

引言

Mapbox已经成为移动地图和数据可视化的代名词。除了已经广泛被程序开发人员和制图师所使用的底图样式工具集以外，他们还提供了Python和JavaScript编写的绘图工具。

通过将这两种语言组合在一个包中，Mapbox发布了新的"MapboxGL-Jupyter" Python模块，允许在Jupter Notebook环境中创建数据即时可视化。此外，Mapbox还提供了允许API访问帐户服务的Mapbox Python SDK。

总之，在Python环境中可以很方便地将Mapbox的工具和服务添加到地理空间应用当中。

Mapbox简介

All things cartographic

Mapbox由Eric Gunderson于2010年创立，其发展迅速，已经成为制图复兴浪潮的领导者。 他们的MapboxGL JavaScript API是一个非常有用的库，用于创建交互式Web地图和数据可视化。 他们为地理空间社区提供了多种开放式地图规范，包括矢量图块。 Mapbox专注于为地图和应用程序开发人员提供自定义底图图块，他们将自己定位为Web地图和移动应用程序的领先软件公司。

Mapbox 工具介绍

Mapbox积极参与开源运动，例如Mapbox开源的负责人Sean Gillies，就是Shapely、Fiona和Rasterio的主要开发者，Mapbox为Python的地理空间分析和可视化贡献颇多。他们的MapboxGL-Jupyter库代表了一种利用其工具套件与其他Python模块（如Pandas / GeoPandas）和多种数据类型（如GeoJSON，CSV，甚至shapefile）相结合的新方法。

除了MapboxGL-Jupyter以外，Mapbox的开源工具还包括基于Web图形库（WebGL）构建的MapboxGL JavaScript库和Mapbox Python SDK。

MapboxGL.js

MapboxGL基于Leaflet.js——一个非常知名的JS地图库构建。Leaflet于2011年发布，支持各种着名的Web地图应用程序，包括Foursquare，Craigslist和Pinterest。 Leaflet的开发人员Vladimir Agafonkin自2013年以来一直为Mapbox工作。

在Leaflet基础上，MapboxGL.js整合了WebGL库，它利用HTML5 canvas标签来支持Web图形而无需插件。MapboxGL.js支持矢量切片，以及平滑缩放和平移的3D环境。 它支持GeoJSON叠加以及标记和形状。事件（包括点击，缩放和平移）可用于触发数据处理功能，使其成为交互式Web地图应用程序的理想选择。

Mapbox Python SDK

Mapbox Python SDK用于访问大多数Mapbox服务，包括方向、地理编码、分析和数据集。 对云端服务的底层访问支持数据编辑和上传、系统管理和基于位置的查询，允许企业将其与本地GIS集成和扩展。

Mapbox Python SDK的安装非常简单，虽然MapboxGL-Jupter工具并不依赖此SDK，但是它在数据上传和查询方面非常有用：

1
pip install mapbox

Mapbox与GIS集成

得益于JavaScript库和新的MapboxGL-Jupyter Python模块，Mapbox工具比以往任何时候都更容易使用。地理空间开发人员和程序员可以将他们的工具集成到现有的GIS工作流程中，也可以创建利用Mapbox产品套件的新地图和应用程序。

Mapbox与CARTO（另一个非常流行的地理空间分析云服务）一样，允许基于帐户的云数据存储。然而，他们的重点不在于分析工具，而在于制图或者说数据可视化。对于大型和小型地图团队，使用Mapbox工具可以降低为交互式Web地图创建和支持自定义底图的成本，并且比其他地图图块选项（如Google Maps API）更加经济。

Mapbox Studio可以轻松创建具有与公司或部门品牌相匹配的制图外观的地图。可以使用现有样式构建底图，并将其覆盖在组织的图层上，也可以设计全新的底图。它甚至允许样式基于图像像素颜色分布的直方图为地理要素分配颜色。

Mapbox入门

注册Mapbox账号

首先，我们需要注册一个账号，然后进入Account Dashboard，从这里可以进入Mapbox Studio, 还可以看到访问API所需要的token（可以创建多个token，也可以刷新token）。此外从Account Dashboard上还可以看到关于不同服务调用PAPI的数据统计等。

access-tokens.png

关于token的用法，在JS代码中，API token被传递给MapboxGL对象，从而获得访问地图和工具的权限，详见 https://www.mapbox.com/help/how-access-tokens-work/ 。通过一个例子简单演示一下，保存下面的代码到html文件中，用自己的token代替“mapboxgl.accessToken =”后面的字符串，然后用浏览器打开：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
<html><head>
<script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.44.1/mapbox-gl.js'>
</script>
<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.44.1/mapbox-gl.css'
rel='stylesheet' />
</head><body>
<div id='map' style='width: 400px; height: 300px;'></div>
<script>
mapboxgl.accessToken =
'pk.eyJ1IjoibG9raXByZXNpZGVud0.8S8l9kH4Ws_ES_ZCjw2i8A';
var map = new mapboxgl.Map({
container: 'map',
style: 'mapbox://styles/mapbox/streets-v9'
});
</script></body></html>

向Mapbox账号添加数据

Mapbox支持用户使用自己的数据。不仅可以自己设置底图图层的样式，还可以将自己的数据添加到图层中，这些都可以使用Mapbox Python SDK以及上传和数据集API以编程方式进行管理。

basic_template

需要注意的是，如果要上传数据，必须创建一个secret API token，其创建方法如前所述，但是在范围（Scope）时，需要选择以下范围以允许数据集和tileset读写功能：

  • DATASETS:WRITE
  • UPLOADS:READ
  • UPLOADS:WRITE
  • TILESETS:READ
  • TILESETS:WRITE

关于如何添加命令，可以查看帮助： https://www.mapbox.com/help/how-uploads-work/

切片集（Tilesets)

切片集是许多栅格图像的集合，可以叠加到底图之上，创建滑动贴图。我们可以里用自己的数据要素，从矢量数据生成切片集，创建自定义的底图。而且，还可以使用Mapbox Python SDK中的Uploader类，将GeoJSON和shapefile文件以切片集的形式上传到个人的云账户中。

数据集（Datasets）

Mapbox中的Datasets是Geojson图层，比切片集更加便于操作。尽管我们可以从Account Dashboard上传数据，但是如果一个数据集大小超过5MB，必须使用datasets API。 具体参见： https://www.mapbox.com/api-documentation/#datasets

示例1：读写GeoJSON数据集

使用MapboxGL-Jupyter和Mapbox Python SDK可以在帐户中存储数据集并将它们连接到其他表格数据。需要注意的是，使用GeoJSON文件需要使用Secret token，但是由于secret token可以具有不同的应用范围（scope），所以要确保secret token的应用范围合理，并且针对不同的应用使用不同的secret tokens。

1
2
3
4
5
6
7
8
9
10
11
12
13
from mapbox import Datasets
import json
# be careful about the secret token
datasets = Datasets(access_token=secret_token)
create_resp = datasets.create(name="Beijing Zone",
                              description = "The area of Beijing")
listing_resp = datasets.list()

dataset_id = [ds['id'] for ds in listing_resp.json()][0]
data = json.load(open(r'data/bei_jing_geo.json', encoding='UTF-8'))
for feature in data['features']:
    resp = datasets.update_feature(dataset_id, feature['id'], feature)
    print(resp.status_code)
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200

beijing_zone

从数据集中读取数据也非常方便，可以直接使用read_dataset()方法读取JSON数据：

1
datasets.read_dataset(dataset_id).json()
{'owner': 'signup-zhang',
 'id': 'cjsorry6y18yd33mnfws1uzjg',
 'name': 'Beijing Zone',
 'description': 'The area of Beijing',
 'created': '2019-02-28T15:14:25.391Z',
 'modified': '2019-02-28T15:15:38.703Z',
 'features': 18,
 'size': 23847,
 'bounds': [115.4237, 39.4395, 117.507, 41.0607]}

还可以删除JSON中的一行数据：

1
resp = datasets.delete_feature(dataset_id, 0)

示例2：以切片集（Tilesets）形式上传数据

通过将切片集添加到自定义底图样式，可以快速加载数据图层和。接下来我们看一下，如何用Mapbox Python SDK将GeoJSON上传为切片集（更多参见https://www.mapbox.com/api-documentation/?language=Python ） ：

1
2
3
4
5
6
7
8
# Use the Uploader class
from mapbox import Uploader
import uuid
set_id = uuid.uuid4().hex
service = Uploader(access_token=secret_token)
with open(r'data/bei_jing_geo.json', 'rb') as src:
    response = service.upload(src, set_id)
print(response)

Http状态码“201 Created”，请求已经被实现，而且有一个新的资源已经依据请求的需要而建立，且其URI已经随Location头信息返回。

upload_as_tileset

Mapbox Studio

即使对于经验丰富的地图作者而言，创建自定义底图也是一个耗时的过程。 为了简化此过程，Mapbo使用Open Street Map（OSM）数据生成预构建的自定义底图，可用于商业和非商业应用程序。通过使用Mapbox Studio，还可以调整这些样式以添加更多自定义的内容。 当然，如果为了特定的需要，也可以从头开始构建底图。

如果仅仅是为了在浏览器中进行手动操作，也没有必要多费篇幅，关键还是要在Jupyter下面操作Mapbox Studio，首先需要安装MapboxGL-Jupyter包：

1
pip install mapboxgl

还需要用到geopandas、pandas等一些包。此外，如果以环境变量的方式将API的token保存起来，可以这样读取token：

1
token = os.getenv("MAPBOX_ACCESS_TOKEN")
1
2
3
4
5
6
7
8
9
10
11
12
# importing data using GeoPandas
import geopandas as gpd
import pysal.viz.mapclassify as mapclassify
import pandas as pd
import os
from mapboxgl.utils import *
from mapboxgl.viz import *
import matplotlib.pylab as plt
import matplotlib as mpl
token = public_token

mpl.rcParams['figure.figsize'] = [15,12]
1
2
import pysal
help(pysal.viz)
Help on package pysal.viz in pysal:

NAME
    pysal.viz

PACKAGE CONTENTS
    mapclassify (package)
    splot (package)

FILE
    c:\programdata\anaconda3\lib\site-packages\pysal\viz\__init__.py

将点添加到Mapbox地图中

1
2
tract_points = gpd.read_file(r'data/tract_points.geojson')
tract_points.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x2b2fe91eb00>

output_41_1

1
tract_points.head()
Total Population Male Population Female Population Percent Male Percent Female geometry
0 4893 2395 2498 0.489475 0.510525 POINT (-122.2925253413085 38.00295823758501)
1 6444 3097 3347 0.480602 0.519398 POINT (-121.7468504262252 36.95049874202599)
2 3736 1757 1979 0.470289 0.529711 POINT (-122.2594490259385 37.8916373721222)
3 4347 2301 2046 0.529331 0.470669 POINT (-122.3459716556987 37.97616227082723)
4 1952 984 968 0.504098 0.495902 POINT (-122.303938838024 37.8657801574032) 
1
2
3
4
5
6
7
8
9
10
11
12
# Generate data breaks and color stops from colorBrewer

# Create the viz from the dataframe
viz = CircleViz('http://127.0.0.1:8888/files/data/tract_points.geojson',
                access_token=public_token,
                center = (-122.5, 37.75),
                zoom = 8,
                color_default = 'red')

viz.show()
with open('circle_viz.html', 'w') as f:
    f.write(viz.create_html())

circle_viz

为了区分不同的数据，我们可以根据某种属性设置color_stops，从而获得颜色渐变的效果.
如果使用GraduatedCircleViz还可以根据某种属性设置radius_stops，从而取得点的大小渐变的效果。例如，我们根据Total Population设置半径大小，根据’Percent Male’设置点的颜色，具体代码如下：

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
radius_measure = 'Total Population'
color_measure = 'Percent Male'

radius_breaks =  [round(tract_points[radius_measure].quantile(q=x*0.1), 2) for x in range(1,10,2)]
radius_stops = create_radius_stops(radius_breaks, 2, 12)
color_breaks = [round(tract_points[color_measure].quantile(q=x*0.2), 2) for x in range(1, 5)]
color_stops = create_color_stops(color_breaks, colors='YlOrRd')

viz2 = GraduatedCircleViz(data='http://127.0.0.1:8888/files/data/tract_points.geojson',   
                          access_token=public_token,
                          center = (-122.5, 37.75),
                          zoom = 8)

viz2.color_stops = color_stops
viz2.color_property = color_measure
viz2.color_function_type='interpolate'
viz2.color_default = 'grey'
viz2.radius_stops = radius_stops
viz2.radius_property = radius_measure
viz2.radius_default=1
viz2.stroke_color='black'
viz2.stroke_black=0.5
viz2.radius_function_type = 'intepolate'
viz2.style='mapbox://styles/mapbox/dark-v9'

viz2.show()
with open('grad_circle.html', 'w') as f:
    f.write(viz2.create_html())

grad_viz

关于viz的更多用法可以参见：

https://github.com/mapbox/mapboxgl-jupyter/blob/master/docs/viz.md

关于数据处理的更多工具参见:

https://github.com/mapbox/mapboxgl-jupyter/blob/master/docs/utils.md

关于颜色渐变的内容参见：

https://github.com/mapbox/mapboxgl-jupyter/blob/master/mapboxgl/colors.py

创建分级统计图（choropleth map）

1
2
3
4
tract_poly= gpd.read_file('data/tracts_bayarea.geojson')
tract_poly['Male Population'] = tract_poly['ACS_15_5YR_S0101_with_ann_Male; Estimate; Total population']
tract_poly = tract_poly[['Male Population','geometry' ]]
tract_poly.to_file('data/tracts_bayarea2.geojson', driver="GeoJSON")
C:\ProgramData\Anaconda3\lib\site-packages\geopandas\io\file.py:108: FionaDeprecationWarning: Use fiona.Env() instead.
  with fiona.drivers():

我们使用ChoroplethViz类实现分级统计图的可视化，并且我们可以在Mapbox Studio上生成定制的底图作为此处的底图：

custom_style

1
2
3
4
5
6
7
8
9
10
11
12
13
vizClor = ChoroplethViz('http://127.0.0.1:8888/files/data/tracts_bayarea2.geojson',
                        access_token=public_token,
                        color_property='Male Population',
                        color_stops=create_color_stops([0, 2000, 3000,5000,7000, 15000],colors='YlOrRd'),
                        color_function_type='interpolate',
                        line_stroke='-',
                        line_color='rgb(128,0,38)',
                        line_width=1,
                        opacity=0.6,
                        center=(-122, 37.75),
                        zoom=9)
vizClor.style='mapbox://styles/signup-zhang/cjsrih04z0qr61fudnlsrfsgv'
vizClor.show()

ch_viz

保存地图为html

通过使用create_html()方法，我们可以把生成的地图保存在本地：

1
2
with open('male_pop.html', 'w') as f:
    f.write(vizClor.create_html())

创建热力图(Heat map)

使用HeatmapViz类创建热力图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
tract_poly = gpd.read_file('data/tract_points.geojson')
measure = 'Female Population'
heatmap_color_stops = create_color_stops([0.01, 0.25, 0.5, 0.75, 1], colors='YlOrRd')
heatmap_radius_stops = [[0, 3], [14, 100]]
color_breaks = [round(tract_poly[measure].quantile(q=x*0.1), 2) for x in range(2,10)]
color_stops = create_color_stops(color_breaks, colors='Spectral')
heatmap_weight_stops = create_weight_stops(color_breaks)
vizheat = HeatmapViz('http://127.0.0.1:8888/files/data/tract_points.geojson',
                     access_token=token,
                     weight_property = "Female Population",
                     weight_stops = heatmap_weight_stops,
                     color_stops = heatmap_color_stops,
                     radius_stops = heatmap_radius_stops,
                     opacity = 0.8,
                     center=(-122, 37.78),
                     zoom=7,
                     below_layer='waterway-label'
)
vizheat.show()
with open('fpop_heat.html', 'w') as f:
    f.write(vizheat.create_html())

heat_map

小结

本文介绍了Mapbox与Python的集成使用，具体包括MapboxGL-Jupyter库和Mapbox Python SDK，利用这两个工具实现了数据的可视化，可以定制化制作底图，绘制分级统计图、热力图等，以及如何进行数据的上传、转换等处理。

算是总结

从第一篇介绍Python的地理空间分析库开始，前前后后用了三个月的时间，基本算是完成了用Python进行地理空间分析的入门学习，甚至可以说，这也只是入门的入门。比如，用GPU数据库（例如MapD）或者Hadoop进行地理空间分析、用Flask或者Django开发地理空间分析的网络应用、其他的GIS服务如QGIS和CARTO……这些都没有介绍。

一味贪多容易消化不良，不如先把主要的一些工具掌握好，再在时机合适或者需要的时候再去学习也不迟。当然，我个人觉得八篇文章，也有了一个框架的雏形：GIS相关概念、Python代码包、GIS数据库、栅格和矢量数据处理、主流GIS软件的结合等等。一窍通百窍通，再要去接触新的东西也不会无所适从。

接下来做些什么呢？拭目以待吧~