Cartopy绘图系列 | 绘制全球温度场+风矢量场

世上总会有“猝不及防的再见和毫无留恋的散场”， 但也会有 “突如其来的遇见和始料未及的欢喜”， 无论何时，记得不要失了好心情！

空间数据的可视化展示是空间数据可视化制图的常见需求。目前常见的地图绘图工具和软件有:

• (1) 支持 多种操作系统的命令行工具 GMT（Generic Mapping Tools）
• (2) ArcGIS、GrDAS等。
• (3) NCAR Command Language): 常用于气象数据读取和可视化分析的高级语言（目前也已经迁移到PyNGL上）
• (4)Python 绘图工具 matplotlib 的扩展包 Basemap（作者在前面的文章中有简单介绍）

除此之外，小编想给大家推荐 Python 的一种制图工具包 Cartopy,（内容比Basemap更加丰富和实用）

cartopy.png

Cartopy简介与安装

Cartopy 是一个开源免费的第三方 Python 扩展包，由英国气象办公室的科学家们开发，支持 Python 2.7 和 Python 3，致力于使用最简单直观的方式生成地图，并提供对 matplotlib 的接口，两者可以完美的协作。
1、Cartopy常用于用于地理空间数据处理，以便生成地图和其他地理空间数据分析。Cartopy利用了强大的PROJ.4、NumPy和Shapely库，并在Matplotlib之上提供了一个编程接口，用于创建出版高质量的地图。
2、Cartopy的关键特性是它面向对象的投影定义，以及在这些投影之间转换点、线、向量、多边形和图像的能力。

为什么要选用Cartopy ?

• Cartopy 的工作流非常简单：设置地图投影，添加地图要素，最后显示地图。
• Basemap绘图包配置相对较繁琐，自定义性不强而且Basemap 将在 2020 年版本的更新。因此，如果你在寻找一个新的Python制图工具包，Cartopy是Basemap官网钦定的继承人，无疑是最佳选择（https://matplotlib.org/basemap/users/intro.html#cartopy-new-management-and-eol-announcement）。

Cartopy安装：

官网教程：https://scitools.org.uk/cartopy/docs/latest/installing.html#installing
1）Anaconda环境
如果你正在使用 Python 的科学计算发行版 Anaconda，安装 Cartopy 非常容易。
命令行输入: conda install -c conda-forge cartopy
2）Windows环境
命令行输入：pip install cartopy
3）可以切换国内像源，安装速度更快
命令行输入：pip install -i https://pypi.tuna.tsinghua.edu.cn/simple cartopy

测试是否安装成功：
启动python命令行工具，输入 import cartopy 如果没有报错，则安装成功！

Cartopy 制图简介

• Cartopy官方网站中列举了许多绘图案例，并且有完整的代码演示。链接https://scitools.org.uk/cartopy/docs/latest/gallery/index.html
  

  image.png

• Cartopy地图绘制的总体流程就是：(1) 选择合适的投影； (2) 添加地图要素（自定义shp\海岸线\边界线等）；(3) 叠加空间数据；(4) 设置地图要素（经纬度标签、比例尺、文本等）

下面介绍几种作者自己总结的常见绘图案例

1、全球地图显示

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : cartopy_script.py
# Author : zengsk in NanJing
# Created: 2020/6/15 0:54

"""
Details: Cartopy package 绘图示例
"""

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt

def main():
    fig = plt.figure(figsize=(8, 8))

    # Label axes of a Plate Carree projection with a central longitude of 180:
    ax1 = fig.add_subplot(2, 1, 1,
                          projection=ccrs.PlateCarree(central_longitude=180))
    ax1.set_global()
    ax1.coastlines() # 添加海岸线

    ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax1.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax1.xaxis.set_major_formatter(lon_formatter)
    ax1.yaxis.set_major_formatter(lat_formatter)

    ax2 = fig.add_subplot(2, 1, 2,
                          projection=ccrs.Robinson(central_longitude=0))

    # make the map global rather than have it zoom in to
    # the extents of any plotted data
    ax2.set_global()
    ax2.stock_img()
    ax2.coastlines()
    ax2.add_feature(cfeature.BORDERS, linestyle='-') # 添加国家边界
    ax2.plot(-0.08, 51.53, 'o', transform=ccrs.PlateCarree())

    plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
    plt.show()

if __name__ == '__main__':
    main()

结果展示

test.jpg

2、显示自定义 shapefile文件**

往往我们需要绘制自定义范围的研究区域，需要绘制指定shapefile文件的边界！

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : cartopy_test2.py

"""
Created by s.k zeng in Chengdu on 2020/07/07
"""

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as mticker
import matplotlib.pyplot as plt


def main():

    # 设置投影
    proj = ccrs.PlateCarree() # 圆柱投影， 默认WGS1984
    extent = [70, 140, 0, 60] # 显示范围
    shpfn = r'E:\GisData\SHP\China_shp\bou2_4m\bou2_4l.shp'
    tpshp = r'E:\GisData\SHP\Tibetan\TP_boundary\TP_Clip.shp'
    glshp = r'E:\GisData\SHP\Global\country.shp'
    reader = shpreader.Reader(shpfn)
    states_provinces = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
    reader = shpreader.Reader(tpshp)
    tpfeat = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
    reader = shpreader.Reader(glshp)
    glfeat = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')

    fig = plt.figure(figsize=(8, 8), frameon=True)
    ax1 = fig.add_subplot(1, 1, 1, projection=proj)
    # Label axes of a Mercator projection without degree symbols in the labels
    # and formatting labels to include 1 decimal place:
    ax1.set_extent(extent, proj)
    ax1.add_feature(glfeat, linewidth=1.5, edgecolor='grey')
    ax1.add_feature(states_provinces, linewidth=1.5, edgecolor='blue')
    ax1.add_feature(tpfeat, linewidth=2.0, edgecolor='red')

    # 设置经纬网以及标签
    ax1.gridlines(proj, draw_labels=False, linewidth=1.2, color='k', alpha=0.5, linestyle='--')
    ax1.set_xticks(np.arange(extent[0], extent[1] + 10, 10), crs=proj)
    ax1.set_yticks(np.arange(extent[2], extent[3] + 10, 10), crs=proj)
    ax1.xaxis.set_major_formatter(mticker.LongitudeFormatter())
    ax1.yaxis.set_major_formatter(mticker.LatitudeFormatter())

    plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
    plt.show()

if __name__ == '__main__':
    main()

结果展示

test.jpg

3、综合案例：Cartopy绘制全球温度场+风矢量场数据**

NCEP/NCAR Reanalysis 1 再分析资料的地面10m风场数据（u v）和地面温度数据绘图示例。（以2020年01月01日为例）
数据可在官网下载：https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : uvwind_plot.py

"""
Details: NCEP/NCAR Reanalysis 1: Surface 10m风速数据绘制风矢量场图
Created by s.k zeng in Chengdu on 2020/07/07
"""

import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.mpl.ticker as mticker


def read_uv():
    '''
    :return: lon, lat, uwind, vwind
    '''
    uwind = r"E:\Scripts\Test\uvdata\uwnd.10m.gauss.2020.nc"
    vwind = r"E:\Scripts\Test\uvdata\vwnd.10m.gauss.2020.nc"
    ufid = nc.Dataset(uwind)
    vfid = nc.Dataset(vwind)
    print(ufid.variables)
    lon = ufid.variables['lon'][:]
    lat = ufid.variables['lat'][:]
    u = ufid.variables['uwnd'][0, :, :] # unit: m/s
    v = vfid.variables['vwnd'][0, :, :]
    return lon, lat, u, v


def read_airtemper():
    temperFn = r"E:\Scripts\Test\uvdata\air.sig995.2020.nc"
    tfid = nc.Dataset(temperFn)
    print(tfid.variables)
    lon = tfid.variables['lon'][:]
    lat = tfid.variables['lat'][:]
    air = tfid.variables['air'][0, :, :] # air temperature unit: degK
    return lon, lat, air


def main():
    lon, lat, u, v = read_uv()
    lon2, lat2, air = read_airtemper()
    proj = ccrs.PlateCarree(central_longitude=180)

    fig = plt.figure(figsize=(9, 6), dpi=300)
    ax = fig.add_subplot(1, 1, 1, projection=proj)
    cs = ax.contourf(lon2, lat2, air, transform=proj, cmap='RdBu_r') # RdBu_r nipy_spectral
    ax.quiver(lon, lat, u, v, transform=ccrs.PlateCarree(), regrid_shape=35, width=0.002)
    ax.coastlines(color='dimgray')
    ax.set_global()

    cbar = fig.colorbar(cs, orientation='vertical', pad=0.02, aspect=20, shrink=0.6)
    cbar.set_label('degK')

    xticks = [-180, -120, -60, 0, 60, 120, 180]
    ax.set_xticks(xticks, crs=proj)
    ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=proj)
    lon_formatter = mticker.LongitudeFormatter(zero_direction_label=True)
    lat_formatter = mticker.LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)

    plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
    plt.show()


if __name__ == '__main__':
    main()

结果展示

test.jpg

感谢支持，本人能力有限，如文章存在错误和高见欢迎留言！