用 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,num2date  
import datetime  
from mpl_toolkits.basemap import Basemap  
import matplotlib.pyplot as plt  
from matplotlib import colors  
import 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_B  
if 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())#快速查看文件变量

代码生成图片

picture.image

NOAA在线绘图对比

picture.image