基于Python的风矢量图绘制

好的伴侣，是精神上的导师，事业上的帮手，开心时的玩伴，失败后的后盾，其他的都只能算是“搭伙过日子”。最怕到后来除了在外面见识了一点世面，其余就是一穷二白的生活。。。。。

风场矢量图可以在ArcGIS中绘制，但是在Python中也是也可以的。主要使用到的方法是quiver()函数，该方法详细的使用可以看看matplotlib的介绍。但有时候会觉得ArcGIS绘制的也不错。。。。。。。。。。

那么在Python中实现起来的也不是什么苦难的事，前提是得有数据

使用的数据来自于AOOS的数据，地址是：[https://thredds.aoos.org/thredds/ncss/WRFSFC.nc/dataset.html]。如果感兴趣的话可以自己去下载数据来做测试

前面import 包那些都是老套路了，读取nc数据，可以选择使用netCDF4或者Xarray，我觉得后者还是比较方便的，特别是在直观呈现数据方面，很有层次。

import timeimport numpy as npfrom numpy import maimport matplotlib as mplfrom matplotlib import cmimport matplotlib.pyplot as pltfrom matplotlib import tickerimport netCDF4import xarrayimport osimport glob

读取数据之后可以切片一个时间节点来可视化看看。

ny=np.array(netcdf.variables['NY'][:])nx=np.array(netcdf.variables['NX'][:])xlon = np.array(netcdf.variables['XLON'][:,:,])xlat = np.array(netcdf.variables['XLAT'][:,:,])time_var = netcdf.variables['TIME']times=netCDF4.num2date(np.array(time_var), time_var.units)
i=-1 # 切最后一个时间点u10 = np.array(netcdf.variables['U10'][i,:,:,])v10 = np.array(netcdf.variables['V10'][i,:,:,])direction = np.array(netcdf.variables['DIR'][i,:,:,])spd = np.array(netcdf.variables['SPD'][i,:,:,])

def create_plot(lon, lat, varname):    i=743    fig = plt.figure(num=None, figsize=(234/15,135/15))    cm.bwr.set_bad('black', 0.8)        plt.title(netcdf.variables['SPD'].long_name + ' on ' + str(times[i]))
    data = np.array(netcdf.variables[varname][i,:,:,])    vmin = np.min(data)    vmax = np.max(data)    img = plt.pcolormesh(lon, lat, data, vmin=vmin, vmax=vmax, cmap = 'bwr')        #加上colorbar    cax=fig.add_axes([0.25, 0.02, 0.5, 0.03]) #[左, 底部, 宽度, 高度]    cb = plt.colorbar(img, cax=cax, orientation='horizontal', extend='both')    cb.locator = ticker.MaxNLocator(nbins=6)    cb.update_ticks()    plt.text(0.5, -2, netcdf.variables[varname].units)            plt.show()    return img

调用自定义函数就可以绘制出切最后一个时间点的图形了。有必要这么复杂吗，Xarray直接就是可以.plot出来的哈哈哈哈哈哈

接着绘制风场矢量图，就是使用quiver()函数，也没有什么其他的了

fig = plt.figure(num=None, figsize=(234/25,135/25))
plt.title(netcdf.variables['SPD'].long_name + ' from CBHAR')plt.xlabel('Longitude'), plt.ylabel('Latitude')
#嫌太密集的话，那就设置一下步长step呗plt.quiver(xlon[::10, ::10], xlat[::10, ::10],           u10[::10


    
, ::10], v10[::10, ::10],           spd[::10, ::10],           cmap='jet')
plt.show()

投影之后再看看

总的来说，要绘制风场矢量图的话还是比较easy的