cartopy绘制中国降雨地图

news2025/7/23 5:41:20

常用的地图可视化的编程工具有 MATLAB、IDL、R、GMT、NCL 等。相比于ArcGIS、QGIS和ArcGISpro用鼠标点来点去，编程绘图也是有很大的优点的，方便，可批量，美观。

大气科学和气象的朋友们一直使用的应该是 NCL，易用性不错，画地图的效果也很好。然而 2019 年初，NCAR 宣布 NCL 将停止更新，并会在日后转为 Python 的绘图包。

此前 Python 最常用的地图包是 Basemap，然而它将于 2020 年被弃用，官方推荐使用 Cartopy 包作为替代。Cartopy 是英国气象局开发的地图绘图包，实现了 Basemap 的大部分功能，还可以通过 Matplotlib 的 API 实现丰富的自定义效果。

似乎Cartopy是目前比较不错的选择，但是实在是太不稳定了，很多功能还不完备，目前的版本是0.21。

首先读取中国降雨格点数据，读成矩阵：

import pandas as pd
import numpy as np
import os
dir = 'D:/Acdemic/xibei/data_1/grid_prcp/'
txtLists = os.listdir(dir)
files = list(filter(lambda x: x[-4:] in ['.txt'], txtLists))
df = pd.DataFrame()
i = 0;
prcp = np.zeros((72, 128, len(files)))
for file in files:
    data = pd.read_table(dir+file, sep='\s+', header=None, index_col=False, skiprows=6)
    print(str(i) + ': ' + file)
    prcp[:, :, i] = np.array(data)
    i += 1
prcp1 = np.mean(prcp, 2) * 12
prcp1[prcp1 < -1000] = None

数据和读取之前介绍过了，参考：

写个函数绘制省界矢量，这个函数也可以用于各类矢量绘制

import cartopy.io.shapereader as shpreader
def add_shp(ax, **kwargs):
    '''
    在地图上画出中国省界的shapefile.
    Parameters
    ----------
    ax : GeoAxes
        目标地图.
    **kwargs
        绘制shape时用到的参数.例如linewidth,edgecolor和facecolor等.
    '''
    proj = ccrs.PlateCarree()
    reader = shpreader.Reader('D:/OneDrive/data/china.shp')
    provinces = reader.geometries()
    ax.add_geometries(provinces, proj, **kwargs)
    reader.close()

绘图：

import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# 导入Cartopy专门提供的经纬度的Formatter
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.ticker as mticker
# 导入底图包
import cartopy.feature as cfeature

fig = plt.figure()
proj = ccrs.PlateCarree()
ax = fig.add_subplot(111, projection=proj)

# 设置经纬度范围,限定为中国
# 注意指定crs关键字,否则范围不一定完全准确
extents = [70, 140, 15, 55]
ax.set_extent(extents, crs=proj)
# 添加各种特征
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAND, edgecolor='black')
#ax.add_feature(cfeature.LAKES, edgecolor='black')
#ax.add_feature(cfeature.RIVERS)
#ax.add_feature(cfeature.BORDERS)
# 添加网格线
ax.gridlines(linestyle='--')

# 设置大刻度和小刻度
ax.set_xticks(np.arange(70, 140 + 20, 20), crs=proj)
ax.set_xticks(np.arange(70, 140 + 10, 10), minor=True, crs=proj)
ax.set_yticks(np.arange(15, 55 + 20, 20), crs=proj)
ax.set_yticks(np.arange(15, 55 + 10, 10), minor=True, crs=proj)

# 利用Formatter格式化刻度标签
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 添加底图
long = np.linspace(72, 136, 128); lat = np.linspace(18, 54, 72)
im = ax.contourf(
    long, lat[::-1], prcp1,
    levels=np.linspace(0, 2000, 11), cmap='RdYlBu_r',
    extend='both', alpha=0.8
)
cbar = fig.colorbar(
    im, ax=ax, shrink=0.9, pad=0.1, orientation='horizontal'
)

cbar.ax.tick_params(labelsize='small')
ax.set_title('Annul precipitaion (mm)')

# 添加矢量
add_shp(ax, lw=1, ec='k', fc='none')
    
plt.show()

extents = [70, 140, 15, 55]是图层的经纬度范围，简单来说是框的大小
ax.add_feature(cfeature.OCEAN)用于添加cartopy的各类内置要素，海洋大陆等等
ax.set_xticks(np.arange(70, 140 + 20, 20), crs=proj)和ax.set_xticks(np.arange(70, 140 + 10, 10), minor=True, crs=proj)分别是主次刻度线，70, 140是x的经度开始和截至，20是间隔
ax.contourf( long, lat[::-1], prcp1, levels=np.linspace(0, 2000, 11), cmap='RdYlBu_r', extend='both', alpha=0.8用于添加栅格，x和y对应经纬度，一般来说nc文件自带这个数据，prcp是栅格矩阵，levels是颜色映射的level。

这样就绘制好了，结果如下：