Plotly绘制地图——基本用法

流量预警:本篇文章需要加载数据,大概需要30MB!

前言

在plotly里自带Maps相关的模块,可以用来绘制分层统计图(Choropleth map)、散点图、气泡图等,但是里面主要是美国或者的地理数据,如果想做中国的可能就不行了,所以还是会有一些限制。简单想一下,如果我们只是要画个大概的地图,有些可视化的功能,其实也不见得那么麻烦,所以这篇文章就介绍如何用plotly绘制一些基本的地图。

1
2
3
4
5
6
7
8
9
10
11
12
13
import xarray
import plotly.offline as py_offline
from plotly.graph_objs import *
import numpy as np
import pandas as pd
import plotly.offline as py_offline
import plotly.graph_objs as go
from plotly.offline import init_notebook_mode
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['figure.figsize']=(16,9)
init_notebook_mode(connected=True)

我们在很早很早以前就用cartopy画过陆地的轮廓地图:

1
2
3
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
plt.show()

output_20_0.png

这个图一看就很基本,但是我在plotly里面并没有找到一个同样简便的方法。参考了plotly Chart Studio里的一些做法,要实现类似的效果,基本的思路是:用matplotlib的Basemap获得轮廓的经纬度数据,将轮廓转换为plotly的"trace",然后用plotly画出图来——这也是我能想到的比较好的办法了。

Basemap转换为plotly

这种做法在许多plotly的作图中已经提到过,我这里也主要是重复他们之前的工作:

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
from mpl_toolkits.basemap import Basemap

# Make shortcut to Basemap object,
# not specifying projection type for this example
m = Basemap()

# Make trace-generating function (return a Scatter object)
def make_scatter(x,y):
return Scatter(
x=x,
y=y,
mode='lines',
line=Line(color="black"),
name=' ' # no name on hover
)

# Functions converting coastline/country polygons to lon/lat traces
def polygons_to_traces(poly_paths, N_poly):
'''
pos arg 1. (poly_paths): paths to polygons
pos arg 2. (N_poly): number of polygon to convert
'''
# init. plotting list
data = dict(
x=[],
y=[],
mode='lines',
line=dict(color="black"),
name=' '
)

for i_poly in range(N_poly):
poly_path = poly_paths[i_poly]

# get the Basemap coordinates of each segment
coords_cc = np.array(
[(vertex[0],vertex[1])
for (vertex,code) in poly_path.iter_segments(simplify=False)]
)

# convert coordinates to lon/lat by 'inverting' the Basemap projection
lon_cc, lat_cc = m(coords_cc[:,0],coords_cc[:,1], inverse=True)

# add plot.ly plotting options
data['x'] = data['x'] + lon_cc.tolist() + [np.nan]
data['y'] = data['y'] + lat_cc.tolist() + [np.nan]

# traces.append(make_scatter(lon_cc,lat_cc))

return [data]

# Function generating coastline lon/lat traces
def get_coastline_traces():
poly_paths = m.drawcoastlines().get_paths() # coastline polygon paths
N_poly = 91 # use only the 91st biggest coastlines (i.e. no rivers)
return polygons_to_traces(poly_paths, N_poly)

# Function generating country lon/lat traces
def get_country_traces():
poly_paths = m.drawcountries().get_paths() # country polygon paths
N_poly = len(poly_paths) # use all countries
return polygons_to_traces(poly_paths, N_poly)

# Get list of of coastline, country, and state lon/lat traces
traces_cc = get_coastline_traces()
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
data = go.Data(traces_cc)

title = u"Coast lines based on Basemap"

axis_style = dict(
zeroline=False,
showline=False,
showgrid=False,
ticks='',
showticklabels=False,
)

layout = Layout(
title=title,
showlegend=False,
hovermode="closest", # highlight closest point on hover
yaxis=go.YAxis(
axis_style,
),
autosize=False,
width=1200,
height=600,
)

fig = Figure(data=data, layout=layout)
py_offline.iplot(fig)
py_offline.plot(fig, include_plotlyjs="cdn", output_type='file', filename='basemap_based.html')

Cartopy转换为plotly

由于Basemap已经不再维护,而cartopy将代替Basemap,所以我想看一下cartopy能不能完成类似的工作。此外,前面我们没有设置地图的坐标系,所以实际上用的就是“空投影”(Plate Carree),那么能不能选择其他投影呢?

答案——当然是肯定的啦!哈哈哈~

1
2
3
import cartopy
import itertools
from cartopy.mpl.patch import geos_to_path

形状提取

首先,我们要获取海岸线、国家边境和经纬度的轮廓线,cartopy提供了直接从"Natural Earth"下载shapefile的接口,并且用cartopy的shapereader模块从shapefile中提取geometries属性,在这里我们选择了低分辨率(110m)的数据(之前用的是中分辨率的,结果实在是有些慢)。

实际上,我们完全可以不用cartopy,直接自己下载shapefile,然后提取里面的几何形状。

在基于Basemap的时候,我们是将matplotlib.collections.PathCollection的转换为plotly里的traces,为了复用前面的代码,我们也将geometries转换为matplotlib.collections.PathCollection

1
2
3
coastline_110m = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m')
countries_110m = cartopy.feature.NaturalEarthFeature('cultural', 'admin_0_countries', '110m')
graticules_110m = cartopy.feature.NaturalEarthFeature('physical', 'graticules_30', '110m')

坐标转换

所以,完成提取和转换以后,接下来的工作就和前面Basemap的做法应该一样了。但是,如果考虑转换坐标系的事情呢?

这里面有一个小技巧,cartopy提供很多坐标转换方法,例如:cartopy.crs.Projection.project_geometry()cartopy.crs.CRS.transform_points()都可以完成坐标的转换,但是我们需要用前者:后者只改变了坐标,而前者改变的不只是点的坐标,还有点的连接关系(设想原来相邻的两个点经过了映射以后分到了地图的两侧)。

1
2
3
4
5
6
def get_traces(feature, target_projection):
geoms = feature.geometries()
geoms = [target_projection.project_geometry(geom, feature.crs)
for geom in geoms]
paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms))
return polygons_to_traces_cartopy(paths, len(paths), feature.crs,target_projection)

地图显示

坐标转换之后,实际上就可以显示了,但是如果这个时候我们看一下坐标轴,会发现地图中的坐标并不是经纬度,所以为了保证显示的效果,我关闭了坐标轴,但是还是不能查看正确的经纬度,所以首先用transform_points将traces的数据转换为经纬度坐标,然后将其设置为hovertext

此外,为了让图形更加美观,我们增加了经纬线,并且用浅灰色标识。

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
# Functions converting coastline/country polygons to lon/lat traces
def polygons_to_traces_cartopy(poly_paths, N_poly, target_crs, orig_crs):
'''
pos arg 1. (poly_paths): paths to polygons
pos arg 2. (N_poly): number of polygon to convert
'''
# init. plotting list
data = dict(
x=[],
y=[],
hovertext=[],
hoverinfo='text',
mode='lines',
line=Line(color="black"),
name=' '
)

for i_poly in range(N_poly):
poly_path = poly_paths[i_poly]

# get the Basemap coordinates of each segment
coords_cc = np.array(
[(vertex[0],vertex[1])
for (vertex,code) in poly_path.iter_segments(simplify=False)]
)

# convert coordinates to lon/lat by 'inverting' the Basemap projection
xy = target_crs.transform_points(orig_crs, coords_cc[:,0],coords_cc[:,1])[:,:2]
lon_cc = coords_cc[:,0]
lat_cc = coords_cc[:,1]
hovertext = xy.round(2)


# add plot.ly plotting options
data['x'] = data['x'] + lon_cc.tolist() + [np.nan]
data['y'] = data['y'] + lat_cc.tolist() + [np.nan]
data['hovertext'] = data['hovertext'] + hovertext.tolist()+ [np.nan]

# traces.append(make_scatter(lon_cc,lat_cc))

return [data]
1
2
3
4
5
6
7
8
9
10
target_projection = ccrs.Robinson(central_longitude=120)


traces_coastline = get_traces(coastline_110m, target_projection)
traces_countries = get_traces(countries_110m, target_projection)
traces_graticules = get_traces(graticules_110m, target_projection)
traces_graticules[0]['line']['color'] = 'rgb(192,192,192)'

traces_countries[0]['line']['dash']='dash'
traces_cc = traces_graticules + traces_countries + traces_coastline

设置Plotly绘图的相关参数:

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
data = Data(traces_cc)

title = u"Coast lines based on Cartopy(CRS:Robinson)"

axis_style = dict(
zeroline=False,
showline=False,
showgrid=False,
ticks='',
showticklabels=False,
)

layout = Layout(
title=title,
showlegend=False,
hovermode="closest",
yaxis=YAxis(
axis_style,
),
xaxis=XAxis(
axis_style,
),
autosize=False,
width=1000,
height=600,
)
fig = go.Figure(data=data, layout=layout)

py_offline.iplot(fig)