python绘制极地投影/cartopy20.0＋解决python极地投影问题/python nc可视化——以北极10m风场为例

在之前的博客。python极地极地投影绘制中，我曾经提到过python cartopy在极地投影中难以添加坐标标签的问题，当时解决方法是使用文本的方式添加，但这种方式的麻烦且并不统一适用。而在之后的可视化过程，我发现python在极地投影时仍会出现一些奇怪的问题：比如，等值线扭曲，出现不规则多边形，风场分布不均匀等等。
然而，这些问题都在最新版的cartopy21中解决，因此被cartopy极地投影折磨了数月而不得不转用m_map的我得到了解脱，下面将以极地2018极地春季的10m风场为例，给大家看看效果。

安装

最新版的cartopy似乎无法通过conda upgrde与conda install的方式直接安装，但是没关系，我们直接去官网下载安装包：cartopy21.0，根据自己的python版本与系统选择即可，这里我选的是win64＋py39。
下载了安装包后，进入anaconda Prompt:

conda uninstall cartopy --force#只删除cartipy，不删除cartopy依赖
conda install cartopy-0.21.0-py39h4915f10_0.tar.bz2

即可。

使用

绘制风场的代码与其他并没有区别，只是由于bug的修复，原来用于普通投影的代码可以正常使用。

import os
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.path as mpath
import cmaps
import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from netCDF4 import Dataset
from wrf import getvar, interplevel, vertcross,vinterp, ALL_TIMES, CoordPair, xy_to_ll, ll_to_xy, to_np, get_cartopy, latlon_coords, cartopy_xlim, cartopy_ylim
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
d=nc.Dataset('F:/ERA5/met_em/use_1.nc')
u=np.array(d.variables['UU'])
v=np.array(d.variables['VV'])
os.chdir('F:/wrfout/wrfout2018/')
filepath='F:/wrfout/wrfout2018/'
file_list = os.listdir(filepath)
#数据导入
files = [os.path.join(filepath,x) for x in file_list]
wrflist=[Dataset(d) for d in file_list]
times = to_np(getvar(wrflist, "times", timeidx=ALL_TIMES))  #
lat2D = to_np(getvar(wrflist[10], "lat"  ))  # units: decimal degrees
lon2D = to_np(getvar(wrflist[10], "lon"  )) 
windspeed=np.sqrt(u**2 + v**2)
#绘图
mpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12   #字体大小
mpl.rcParams["axes.linewidth"] = 1 
proj =ccrs.NorthPolarStereo(central_longitude=0)#设置地图投影
#在圆柱投影中proj = ccrs.PlateCarree(central_longitude=xx)
leftlon, rightlon, lowerlat, upperlat = (-180,180,60,90)#经纬度范围

img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig1 = plt.figure(figsize=(12,10))#设置画布大小
f1_ax1 = fig1.add_axes([0.2, 0.3, 0.5, 0.5],projection = ccrs.NorthPolarStereo(central_longitude=0))#绘制地图位置
#注意此处添加了projection = ccrs.NorthPolarStereo()，指明该axes为北半球极地投影
f1_ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=1, color='grey',linestyle='--')
f1_ax1.set_extent(img_extent, ccrs.PlateCarree())
f1_ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
f1_ax1.set_boundary(circle, transform=f1_ax1.transAxes)
quiver = f1_ax1.quiver(lon2D, lat2D, musp, mvsp,windspeed, pivot='tail',width=0.002, scale=50, cmap=plt.cm.jet, headwidth=4,
                       regrid_shape=35,alpha=1,transform=ccrs.PlateCarree())
f1_ax1.quiverkey(quiver, 0.91, 1.03, 3, "3m/s",labelpos='E', coordinates='axes', fontproperties={'size': 10,'family':'Times New Roman'})
plt.show()

可以看出并无不同，唯一值得一提的是regrid_shape选项，它可以帮助我们绘制风场时不那么密集。
效果：

可以看出比之前的bug都有所修复。