地理信息与地图¶

地理数据可视化是数据科学中一种十分常见的可视化类型。 Mat plot lab做此类可视化的主要工具是。base map工具箱。它是。 Matplotlab中专门mpl_toolkits命名空间里众多工具箱之一。坦白地说。 Basemap用起来有点笨重。在处理比较复杂的。地图可视化任务时。更现代化的解决方案。如。谷歌的地图API。由于中国大陆无法接入。谷歌地图。同时basemap。符合Python用户的使用习惯。我们在这里演示一些利用base map工具箱绘制地图的可视化示例。

Basemap安装起来非常简单。如果使用conda命令,只需敲入以下代码。

$ conda install basemap

同样在导出的过程中也需要只增加一行代码。

In [2]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

%matplotlib inline import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap

一旦装入basemap，显示地图的过程将会变得十分简单。

In [5]:

#help(Basemap)

#help(Basemap)

In [81]:

help(plt.gci)

help(plt.gci)

Help on function gci in module matplotlib.pyplot:

gci() -> 'ScalarMappable | None'
    Get the current colorable artist.

    Specifically, returns the current `.ScalarMappable` instance (`.Image`
    created by `imshow` or `figimage`, `.Collection` created by `pcolor` or
    `scatter`, etc.), or *None* if no such instance has been defined.

    The current image is an attribute of the current Axes, or the nearest
    earlier Axes in the current figure that contains an image.

    Notes
    -----
    Historically, the only colorable artists were images; hence the name
    ``gci`` (get current image).

In [9]:

plt.figure(figsize=(8, 8))
m = Basemap(projection='ortho', resolution=None, lat_0=50, lon_0=-100)
m.bluemarble(scale=0.8);

plt.figure(figsize=(8, 8)) m = Basemap(projection='ortho', resolution=None, lat_0=50, lon_0=-100) m.bluemarble(scale=0.8);

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Basemap的相关参数简单介绍如下：

Basemap具有完整的，连续的图形体系。所有的数据调用都可以通过相应的参数来实现。如投影类型、图形宽度和高度，以及经纬度。还可以在图形上，设置相应的坐标点。

In [94]:

fig = plt.figure(figsize=(8, 8))
m = Basemap(projection='lcc', resolution=None,
            width=8E6, height=8E6, 
            lat_0=45, lon_0=100,)
m.etopo(scale=0.5, alpha=0.5)

# Map (long, lat) to (x, y) for plotting
x, y = m(116, 40)
plt.plot(x, y, 'ok', markersize=5)
plt.text(x, y, ' Beijing', fontsize=12);

fig = plt.figure(figsize=(8, 8)) m = Basemap(projection='lcc', resolution=None, width=8E6, height=8E6, lat_0=45, lon_0=100,) m.etopo(scale=0.5, alpha=0.5) # Map (long, lat) to (x, y) for plotting x, y = m(116, 40) plt.plot(x, y, 'ok', markersize=5) plt.text(x, y, ' Beijing', fontsize=12);

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

这个例子让我们发现，只需要几行简单的Python代码就可以画出地理可视化图。后面我们将深入的介绍bitmap的主要特征。并通过一些可视化的地图示例进行演示。

地图投影¶

当你使用地图时。首先要做的就是确定地图的投影类型。你可能已经知道。像地球这样的球体。可以通过球面透视法将三维球面投影成一个二维平面，不会造成变形，也不会破坏其连续性。这些投影类型随着人类历史的进程逐步发展起来，现在已经有了许多选择。根据地图投影类型的不同用途，有一些地图特征。如方向，面积距离形状或者其他要素。都需要通过设置实现。

Basemap程序包里包含了几10种投影类型。所有的投影都有一个简便的格式吗？下面对一些常用的投影类型进行简单演示。

In [98]:

from itertools import chain

def draw_map(m, scale=0.2):
    # draw a shaded-relief image
    m.shadedrelief(scale=scale)
    
    # lats and longs are returned as a dictionary
    lats = m.drawparallels(np.linspace(-90, 90, 13))
    lons = m.drawmeridians(np.linspace(-180, 180, 13))

    # keys contain the plt.Line2D instances
    lat_lines = chain(*(tup[1][0] for tup in lats.items()))
    lon_lines = chain(*(tup[1][0] for tup in lons.items()))
    all_lines = chain(lat_lines, lon_lines)
    
    # cycle through these lines and set the desired style
    for line in all_lines:
        line.set(linestyle='-', alpha=0.3, color='w')

from itertools import chain def draw_map(m, scale=0.2): # draw a shaded-relief image m.shadedrelief(scale=scale) # lats and longs are returned as a dictionary lats = m.drawparallels(np.linspace(-90, 90, 13)) lons = m.drawmeridians(np.linspace(-180, 180, 13)) # keys contain the plt.Line2D instances lat_lines = chain(*(tup[1][0] for tup in lats.items())) lon_lines = chain(*(tup[1][0] for tup in lons.items())) all_lines = chain(lat_lines, lon_lines) # cycle through these lines and set the desired style for line in all_lines: line.set(linestyle='-', alpha=0.3, color='w')

圆柱投影（cylindrical projections）¶

圆柱投影是最简单的地图投影类型。经度线和纬度线分别映射成竖直线和水平线。采用这种投影类型的话，赤道区域的显示效果非常好。但是南北极附近的区域就会产生严重变形。由于纬度线的间距会因圆柱形投影的不同而不同。在下面这个等距圆柱投影示例中，不同纬度在子午线方向的间距保持不变。

另外两种常用的投影类型为墨卡托投影和圆柱等积投影，分别为 projection='merc' 和projection='cea'。

In [103]:

fig = plt.figure(figsize=(8, 6), edgecolor='w')
m = Basemap(projection='cyl', resolution=None,
            llcrnrlat=-90, urcrnrlat=90,
            llcrnrlon=-180, urcrnrlon=180, )
draw_map(m)

fig = plt.figure(figsize=(8, 6), edgecolor='w') m = Basemap(projection='cyl', resolution=None, llcrnrlat=-90, urcrnrlat=90, llcrnrlon=-180, urcrnrlon=180, ) draw_map(m)

The additional arguments to Basemap for this view specify the latitude (lat) and longitude (lon) of the lower-left corner (llcrnr) and upper-right corner (urcrnr) for the desired map, in units of degrees.

伪圆柱投影（Pseudo-cylindrical projections）¶

Pseudo-cylindrical projections relax the requirement that meridians (lines of constant longitude) remain vertical; this can give better properties near the poles of the projection. 摩尔维德投影是一种常见的伪圆柱投影 (projection='moll') 。 It is constructed so as to preserve area across the map: though there are distortions near the poles, the area of small patches reflects the true area. Other pseudo-cylindrical projections are the sinusoidal (projection='sinu') and Robinson (projection='robin') projections.

In [40]:

fig = plt.figure(figsize=(8, 6), edgecolor='w')
m = Basemap(projection='moll', resolution=None,
            lat_0=0, lon_0=0)
draw_map(m)

fig = plt.figure(figsize=(8, 6), edgecolor='w') m = Basemap(projection='moll', resolution=None, lat_0=0, lon_0=0) draw_map(m)

The extra arguments to Basemap here refer to the central latitude (lat_0) and longitude (lon_0) for the desired map.

透视投影方法（Perspective projections）¶

Perspective projections are constructed using a particular choice of perspective point, similar to if you photographed the Earth from a particular point in space (a point which, for some projections, technically lies within the Earth!). One common example is the orthographic projection (projection='ortho'), which shows one side of the globe as seen from a viewer at a very long distance. As such, it can show only half the globe at a time. Other perspective-based projections include the gnomonic projection (projection='gnom') and stereographic projection (projection='stere'). These are often the most useful for showing small portions of the map.

Here is an example of the orthographic projection:

In [106]:

fig = plt.figure(figsize=(8, 8))
m = Basemap(projection='ortho', resolution=None,
            lat_0=50, lon_0=0)
draw_map(m);

fig = plt.figure(figsize=(8, 8)) m = Basemap(projection='ortho', resolution=None, lat_0=50, lon_0=0) draw_map(m);

圆锥投影（Conic projections）¶

先将地图投影成一个圆锥体，然后将其展开。这样子会使远离圆锥点的地方变形很大。兰伯特投影 (projection='lcc')就是圆锥型投影。

In [121]:

fig = plt.figure(figsize=(8, 8))
m = Basemap(projection='lcc', resolution=None,
            lon_0=0, lat_0=50, lat_1=45, lat_2=55,
            width=1.6E7, height=1.2E7)
draw_map(m)

fig = plt.figure(figsize=(8, 8)) m = Basemap(projection='lcc', resolution=None, lon_0=0, lat_0=50, lat_1=45, lat_2=55, width=1.6E7, height=1.2E7) draw_map(m)

其他投影¶

参看Basemap程序包的说明文档： Basemap package.

Drawing a Map Background¶

bluemarble() 和 shadedrelief() 方法可以画出地球投影， drawparallels() 和 drawmeridians() 方法可以画出经线与纬线。basemap有一些其他独特功能。

• 物理边界与水域

  • drawcoastlines():绘制大陆海岸线
  • drawlsmask(): 为陆地和海洋设置重填色，这样可以在陆地和海洋投影其他图像。
  • drawmapboundary(): 绘制地图边界，包括为海洋充填颜色。
  • drawrivers(): 绘制河流。
  • fillcontinents(): 用一种颜色充填大陆，另一种颜色充填湖泊（可选）。

• 政治边界

  • drawcountries(): 绘制国界线
  • drawstates(): 绘制州边界图
  • drawcounties(): 绘制美国县域图

• 地图功能

  • drawgreatcircle(): 两点 之间绘制一个大圆
  • drawparallels(): 绘制经纬线
  • drawmeridians(): 绘制经线
  • drawmapscale(): 在地图上绘制一个线性比例尺

• 地球投影

  • bluemarble(): 使用 NASA的蓝色地球图像
  • shadedrelief(): 在地图上绘制地貌渲染图
  • etopo(): 在地图上绘制地形渲染图
  • warpimage():将用户提供的图像投影到地图上

如果使用边界特征，就必须在创建图像时候设置分辨率。Basemap使用 resolution 参数来设定，参数包括边界的细致程度, 原始 'c' (crude), 低的'l' (low), 中等'i' (intermediate), 高的'h' (high), 全的'f' (full), 或者使用 None表示不出现边界。 这个选择很重要，高分辨率在运行全球地图时会非常非常慢。

下面实例演示陆地与海洋的边界与resolution参数设置的关系。这是苏格兰美丽的天空岛，用一张9万*12万公里的区域可以展现出来。

In [137]:

fig, ax = plt.subplots(1, 2, figsize=(12, 8))

for i, res in enumerate(['l', 'f']):
    m = Basemap(projection='gnom', lat_0=57.3, lon_0=-6.2,
                width=90000, height=120000, resolution=res, ax=ax[i])
    m.fillcontinents(color="#FFDDCC", lake_color='#DDEEFF')
    m.drawmapboundary(fill_color="#DDEEFF")
    m.drawcoastlines()
    ax[i].set_title("resolution='{0}'".format(res));

fig, ax = plt.subplots(1, 2, figsize=(12, 8)) for i, res in enumerate(['l', 'f']): m = Basemap(projection='gnom', lat_0=57.3, lon_0=-6.2, width=90000, height=120000, resolution=res, ax=ax[i]) m.fillcontinents(color="#FFDDCC", lake_color='#DDEEFF') m.drawmapboundary(fill_color="#DDEEFF") m.drawcoastlines() ax[i].set_title("resolution='{0}'".format(res));

图上可以看出，低分辨率的海岸线不适合，但高分辨率的还不错。但是，低分辨率呈现的速度很快，因此，我们需要在速度与效果之间做好权衡。

在地图上画出数据¶

Basemap工具箱最实用的功能就是以地图为背景画出各种数据。使用任意plt函数就能在地图上绘制基本图形和文字，你可以将经纬度变成笛卡尔坐标中（x,y）值。

在basemap的实例中，大部分功能与matplotlib中类似，但增加了一个布尔型参数 latlon, 如果设置为 True，则允许直接接受经纬度数值，不需要投影到 (x, y) 坐标轴。

与地图有关的方法如下：

• contour()/contourf() : 绘制登高线，或者充填等高线
• imshow(): 绘制一个图像
• pcolor()/pcolormesh() : 对规则或者不规则的网格线绘制伪彩色图
• plot(): 绘制线图或者点图
• scatter(): 绘制点或者标记
• quiver(): 绘制向量图（箭头图形）
• barbs(): 绘制风羽图（wind barbs）.
• drawgreatcircle(): 绘制一个大圆

示例：加州城市¶

之前介绍过加州城市位置的图形，我们仍然使用它，绘制在地图之上。

In [114]:

import pandas as pd
cities = pd.read_csv('data/california_cities.csv')

# Extract the data we're interested in
lat = cities['latd'].values
lon = cities['longd'].values
population = cities['population_total'].values
area = cities['area_total_km2'].values

import pandas as pd cities = pd.read_csv('data/california_cities.csv') # Extract the data we're interested in lat = cities['latd'].values lon = cities['longd'].values population = cities['population_total'].values area = cities['area_total_km2'].values

In [116]:

cities

cities

Out[116]:

	Unnamed: 0	city	latd	longd	elevation_m	elevation_ft	population_total	area_total_sq_mi	area_land_sq_mi	area_water_sq_mi	area_total_km2	area_land_km2	area_water_km2	area_water_percent
0	0	Adelanto	34.576111	-117.432778	875.0	2871.0	31765	56.027	56.009	0.018	145.107	145.062	0.046	0.03
1	1	AgouraHills	34.153333	-118.761667	281.0	922.0	20330	7.822	7.793	0.029	20.260	20.184	0.076	0.37
2	2	Alameda	37.756111	-122.274444	NaN	33.0	75467	22.960	10.611	12.349	59.465	27.482	31.983	53.79
3	3	Albany	37.886944	-122.297778	NaN	43.0	18969	5.465	1.788	3.677	14.155	4.632	9.524	67.28
4	4	Alhambra	34.081944	-118.135000	150.0	492.0	83089	7.632	7.631	0.001	19.766	19.763	0.003	0.01
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
477	477	Yountville	38.403056	-122.362222	30.0	98.0	2933	1.531	1.531	0.000	3.966	3.966	0.000	0.00
478	478	Yreka	41.726667	-122.637500	787.0	2582.0	7765	10.053	9.980	0.073	26.036	25.847	0.188	0.72
479	479	YubaCity	39.134722	-121.626111	18.0	59.0	64925	14.656	14.578	0.078	37.959	37.758	0.201	0.53
480	480	Yucaipa	34.030278	-117.048611	798.0	2618.0	51367	27.893	27.888	0.005	72.244	72.231	0.013	0.02
481	481	YuccaValley	34.133333	-116.416667	1027.0	3369.0	20700	40.015	40.015	0.000	103.639	103.639	0.000	0.00

482 rows × 14 columns

接下来，我们设定投影，绘制图形，画出图例。

In [135]:

# 1. 画出地图背景
fig = plt.figure(figsize=(8, 8))
m = Basemap(projection='lcc', resolution='h', 
            lat_0=37.5, lon_0=-119,
            width=1E6, height=1.2E6)
m.shadedrelief()
m.drawcoastlines(color='gray')
m.drawcountries(color='gray')
m.drawstates(color='gray')

# 2. 绘制城市散点以及人口数量的反应面积
m.scatter(lon, lat, latlon=True,
          c=np.log10(population), s=area,
          cmap='Reds', alpha=0.5)

# 3. create colorbar and legend
plt.colorbar(label=r'$\log_{10}({\rm population})$')
plt.clim(3, 7)

# 使用虚拟点绘制图例 
for a in [100, 300, 500]:
    plt.scatter([], [], c='k', alpha=0.5, s=a,
                label=str(a) + ' km$^2$')
plt.legend(scatterpoints=1, frameon=False,
           labelspacing=1, loc='lower left');

# 1. 画出地图背景 fig = plt.figure(figsize=(8, 8)) m = Basemap(projection='lcc', resolution='h', lat_0=37.5, lon_0=-119, width=1E6, height=1.2E6) m.shadedrelief() m.drawcoastlines(color='gray') m.drawcountries(color='gray') m.drawstates(color='gray') # 2. 绘制城市散点以及人口数量的反应面积 m.scatter(lon, lat, latlon=True, c=np.log10(population), s=area, cmap='Reds', alpha=0.5) # 3. create colorbar and legend plt.colorbar(label=r'$\log_{10}({\rm population})$') plt.clim(3, 7) # 使用虚拟点绘制图例 for a in [100, 300, 500]: plt.scatter([], [], c='k', alpha=0.5, s=a, label=str(a) + ' km$^2$') plt.legend(scatterpoints=1, frameon=False, labelspacing=1, loc='lower left');

GeoPandas¶

本部分内容不在详细讲述，有需要的同学请查阅官方文档。