用 Python 读取气象环境数据并绘图

NCEP/NCAR再分析资料由美国国家环境预报中心（NCEP）和美国国家大气研究中心（NCAR）联合制作，数据涵盖全球范围内多个高度层次上的多种气象要素，通常被作为分析与研究天气气候异常的诊断资料，是气象工作者最常使用的气象数据之一。

该数据集以netCDF(netware Common Data Form)的形式存储，使用GrADS，NCL等常用气象软件可读取文件，对其中相关数据进行加工与处理，制作天气图。本文将以代码形式介绍如何使用python来读取此类数据，并对数据进行处理，实现自定义时段内平均环流场的可视化。

1.安装支持库

首先确保运行环境已安装netCDF4库文件，可用如下命令查询/安装。

查询：pip list (或 conda list)安装：pip install netCDF4 (或 conda install netCDF4)

其次需手动下载Basemap库文件并安装，下载地址：https://www.lfd.uci.edu/~gohlke/pythonlibs/。(需要先安装pyproj，地址相同)

安装：pip 路径 安装包文件名

2.数据及环境

本文示例中将使用2018年逐日全球高度场文件，包含高低空共17个等压面上的位势高度数据，水平分辨率2.5*2.5（文件可在NOAA官网中下载）。运行环境为python-3.6.5。

3.代码示例

from netCDF4 import Dataset,num2dateimport datetimefrom mpl_toolkits.basemap import Basemapimport matplotlib.pyplot as pltfrom matplotlib import colorsimport numpy as np# 自定义时间段month_B=7;day_B=1;month_E=8;day_E=31  #这里选取盛夏 # 自定义纬度范围lat_B=10;lat_E=90;lon_B=30;lon_E=190  #这里选取欧亚# 自定义等压面层次level=500 # 打开文件filename='e:/hgt.2018.nc'fh=Dataset(filename,mode='r')# 映射时间索引位置times_transformed=num2date(fh.variables['time'][:],units=fh.variables['time'].units,calendar='gregorian')times_transformed_date=list(map(lambda x:x.date(),times_transformed))datetime_B=datetime.date(times_transformed[0].year,month_B,day_B)datetime_E=datetime.date(times_transformed[0].year,month_E,day_E)time_B_index=times_transformed_date.index(datetime_B)time_E_index=times_transformed_date.index(datetime_E)# 映射经纬度索引位置lat_index_B=fh.variables['lat'][:].tolist().index(lat_B)lat_index_E=fh.variables['lat'][:].tolist().index(lat_E)lon_index_B=fh.variables['lon'][:].tolist().index(lon_B)lon_index_E=fh.variables['lon'][:].tolist().index(lon_E)if lat_index_B>lat_index_E:    lat_index_B,lat_index_E=lat_index_E,lat_index_Bif lon_index_B>lon_index_E:    lon_index_B,lon_index_E=lon_index_E,lon_index_B# 映射高度索引位置level_index=fh.variables['level'][:].tolist().index(level)# 对绘图变量进行赋值data_plot=fh.variables['hgt'][time_B_index:time_E_index,level_index,lat_index_B:lat_index_E+1,lon_index_B:lon_index_E+1]lons = fh.variables['lon'][lon_index_B:lon_index_E+1]lats = fh.variables['lat'][lat_index_B:lat_index_E+1]# 求要素的时间平均并保存要素单位data_plot=np.average(data_plot,0) Met_Var_units = fh.variables['hgt'].units+' ('+fh.variables['hgt'].var_desc+')'# 新建地图并设定经纬度范围m = Basemap(llcrnrlat = lat_B, urcrnrlat = lat_E, llcrnrlon = lon_B, urcrnrlon = lon_E) # 网格化经纬度并形成坐标矩阵lon, lat = np.meshgrid(lons, lats)xi, yi = m(lon, lat)# 设置等值线范围及间隔levels=range(5320,5961,40) # 创建填色图层，alpha为透明度0-1，Vmax表示大于该值使用同样的填色以突出副热带高压主体cs = m.contourf(xi, yi, np.squeeze(data_plot),levels=levels,vmin=5400,vmax=5880,alpha=1,cmap="RdBu_r") # 创建等值线图层cs2 = m.contour(xi, yi, np.squeeze(data_plot),levels=[5880],colors='black',linewidths=0.5) # 绘制经纬线背景网格 label为经纬度标签位置[左,右,上,下]m.drawparallels(np.arange(-90., 91., 20.), labels=[1,0,0,0], fontsize=10)m.drawmeridians(np.arange(-180., 181., 40.), labels=[0,0,0,1], fontsize=10)# 添加海岸线国界等m.drawcoastlines()m.drawstates()m.drawcountries()# 添加中国省界及台湾省地理信息(shape文件需自行下载)m.readshapefile('C:/matplotlib/gadm36_CHN_shp/gadm36_CHN_1','states',drawbounds=True)m.readshapefile('C:/matplotlib/gadm36_TWN_shp/gadm36_TWN_1','taiwan',drawbounds=True)# 添加图例及单位 pad调整相对位置cbar = m.colorbar(cs, location='bottom', pad="10%")cbar.set_label(Met_Var_units)# 添加标题 plt.title('%shPa   %s-%s'%(level,datetime_B.strftime('%Y%m%d'),datetime_E.strftime('%Y%m%d')))# 执行可视化(弹出窗口)plt.show()# 执行可视化(不弹出窗口，直接保存)# plt.savefig('D:/fig.png')plt.close()#--其他常用命令--print(fh)                 #查看文件描述print(fh.variables)       #查看所有变量描述print(fh.variables['hgt'])#查看特定变量描述print(fh.variables.keys())#快速查看文件变量

代码生成图片

NOAA在线绘图对比