Cartopy绘图系列 | 绘制全球温度场+风矢量场
2020-07-09 本文已影响0人
zengsk
世上总会有“猝不及防的再见和毫无留恋的散场”， 但也会有 “突如其来的遇见和始料未及的欢喜”， 无论何时，记得不要失了好心情！
空间数据的可视化展示是空间数据可视化制图的常见需求。目前常见的地图绘图工具和软件有:
- (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
