Python地图绘制模块cartopy入门

Python中有很多用于地图绘制的模块，比如老牌的pyGMT，曾经主流的basemap。相比之下，cartopy是英国气象局开发的一款用于地理空间数据处理的python库，支持创建并发布高质量地图，随着basemap的停止维护，cartopy已成Python中地理绘图的首选模块。

海岸线

在cartopy中，最简单的地图即为海岸线的线条图，如下图所示

代码为

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.subplot(projection=ccrs.PlateCarree())
type(ax)    # <class 'cartopy.mpl.geoaxes.GeoAxes'>
ax.coastlines()
plt.show()

【PlateCarree】是可利投影，圆柱投影的原理是将地物投影到圆柱面上再展开，将这种投影作为一种坐标映射提供给坐标轴，便可得到一个cartopy定义的坐标轴对象GeoAxes。

【GeoAxes】显然继承自Axes，是绘图窗口中的坐标映射，并且在Axes的基础之上，封装了一些用于地图绘制的函数。

更多投影，可见projections.

地形图

【coastlines】即海岸线，由此得到的地图是极致的线条风格，光秃秃得什么都没有，并不好看。相比之下，下面这张图看上去十分不错

其绘图代码如下，【stock_img】为其生成函数，其功能是为地图添加低分辨率的地形图背景，但美中不足的是这个图就只有这一张。

ax = plt.subplot(projection=ccrs.PlateCarree())
ax.stock_img()
plt.show()

重载绘图函数

得益于cartopy中的图窗继承了matplotlib，从而后者封装的诸多函数也可以为其所用，例如在地形图上绘制风场，如下图所示。

其中风场数据通过下面的代码创建

import numpy as np
def sample_data(shape=(20, 30)):
    crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
    x = np.linspace(311.9, 391.1, shape[1])
    y = np.linspace(-23.6, 24.8, shape[0])
    X, Y = np.deg2rad(np.meshgrid(x, y))
    u = 10 * (2 * np.cos(2 * X + 3 * (Y + 30)) ** 2)
    v = 20 * np.cos(6 * np.deg2rad(X))
    return x, y, u, v, crs

绘图代码如下

ax = plt.subplot(projection=ccrs.PlateCarree())
ax.set_extent([-90, 80, 10, 85], crs=ccrs.PlateCarree())
ax.stock_img()
ax.coastlines()

x, y, u, v, vector_crs = sample_data(shape=(10, 14))
ax.barbs(x, y, u, v, length=5,
    sizes=dict(emptybarb=0.25, spacing=0.2, height=0.5),
    lw=0.95, transform=vector_crs)
plt.show()

首先，在窗口中新建一个坐标映射，然后通过【set_extent】设置一个元组$(x_0,x_1,y_0,y_1)$，以限制地图的显示范围。

【barbs】是plt中用于绘制风场的函数，但在cartopy中进行了重载，其中transform参数就是重载之后新增的，其目的自然是为了适应新的坐标映射。类似的重载图像如下表所示。

绘图特征

陆地和海洋的边缘，是最清晰的地理特征，所以海岸线图也是cartopy中最基础的地图。但除此之外，地图中还有其他地理特征，cartopy提供了相应的add_feature，均被封装在cartopy.feature中，如下表所示

这六种地图的区别如下

绘图代码如下

from cartopy.feature import *

proj = ccrs.PlateCarree()

feaDct = {
    "BORDERS":BORDERS,
    "COASTLINE":COASTLINE,
    "LAKES":LAKES,
    "LAND":LAND,
    "OCEAN":OCEAN,
    "RIVERS":RIVERS
}

fig = plt.figure()
for i,key in enumerate(feaDct, 1):
    ax = fig.add_subplot(2,3,i,projection=proj)
    ax.add_feature(feaDct[key])
    plt.title(key)

plt.show()

这些不同的特征，如果单挑出来，其实并不好看，其组合在一起之后，效果如下。这是一张非洲地图，途中给出了海岸线、国界线以及河流湖泊等地图元素。

绘图代码如下。

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-20, 60, -40, 40], crs=ccrs.PlateCarree())

ax.add_feature(LAND)
ax.add_feature(OCEAN)
ax.add_feature(COASTLINE)
ax.add_feature(BORDERS, linestyle=':')
ax.add_feature(LAKES, alpha=0.5)
ax.add_feature(RIVERS, color='blue')

plt.show()