Python气象信息可视化方法汇总

摘要

        气象数据可视化是一种将气象数据通过图形、图像等视觉化手段进行呈现，以便更好地理解、分析和预测天气状况的系统。它的发展历程可以追溯到20世纪70年代，随着计算机技术和大数据处理能力的不断提升，气象数据可视化系统也得到了迅速发展，成为气象领域中不可或缺的一部分。

        我第一次接触气象信息可视化是在大气物理这门课。在大气物理课上我第一次使用Python metpy库进行埃玛图的绘制。接下来就是在天气学原理这门课，在课后我们实现了包括平均温度场经向剖面图、北半球平均500hPa高度场、1000hPa平均风场等在内的十余幅天气学插图复现。随后我们专业开设了气象数据可视化课程，在这里我学到了更多可视化方法以及气象数据处理技巧。包括利用scipy库进行数据处理以及更多关于cartopy库的技巧。我将使用到的方法进行总结，形成了此篇文章。

关键词

可视化 Python 气象可视化 数据处理

正文

        可视化是一门涉及计算机科学、图形学、人机交互、认知科学等多个领域的交叉学科，其应用领域广泛，包括数据分析、科学研究、医学成像、工程设计等。而气象信息可视化是一种将气象数据以图形或图像的形式呈现出来的技术。它可以帮助人们更好地理解和分析气象数据，发现其中的规律和趋势，为气象预测和决策提供支持。气象信息可视化包括多种形式，如雷达图、气象卫星图、气象地图等。这些图形或图像可以显示各种气象数据，如降雨强度、风速、云层、气压等。通过这些可视化图像，人们可以更直观地了解气象情况，预测天气趋势，制定应对措施。此外，三维气象数据可视化系统也是气象信息可视化的一种重要形式。它基于地理信息绘制和大量的卫星影像数据和地形数据，对具体地理位置的数据进行实时、准确、直观地观察和统计，实现气象部门所需的气象监测与预警。三维气象数据可视化系统还可以根据行业特征及自身需求，灵活调配预警逻辑。在学习完气象信息可视化这门专业课后，我学到了更多可视化方法以及气象数据处理技巧。接下来我将对我最近在气象信息可视化中所学所写进行简要的汇总。

        此篇文章主要分为四个部分，我会首先从Python语言中库的使用入手，进行简单的介绍，再用一些实例进行绘图步骤的分析：包括数据的读取、处理以及图像的绘制过程。

一、库的介绍

• netCDF4

        NetCDF(network Common Data Form)网络通用数据格式是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作 NetCDF 数据集。netCDF4库是Python中专门对.nc文件进行处理的库。

• Matplotlib

        Matplotlib 是Python的一个综合性的库，可创建静态的、动画的和可交互的可视化图形图像。库中大部分组件都是类，可以创建并修改这些类的实例以创建你想要的图形。拥有可交互式接口，可以在图形显示区域直接对图形进行操作，例如点击，拖拽等。支持多种输出格式：包括但不限于PNG，JPG，SVG，PDF，PS等。

• NumPy

        NumPy(Numerical Python) 是 Python 语言的一个扩展程序库，支持大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。

        NumPy 的前身 Numeric 最早是由 Jim Hugunin 与其它协作者共同开发，2005 年，Travis Oliphant 在 Numeric 中结合了另一个同性质的程序库 Numarray 的特色，并加入了其它扩展而开发了 NumPy。NumPy 为开放源代码并且由许多协作者共同维护开发。

• Scipy

        Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。

• Cartopy

        Cartopy地图绘图包——“专为地理空间数据处理而设计，以生成地图和其他地理空间数据分析。”，是在PROJ、pyshp、shapely、GEOS等Python包的基础上编写的，在安装时，需要同时安装相关的依赖包。Cartopy包对Matplotlib包的功能进行了扩展，两者结合使用能绘制各种地图。它的特点是：①面向对象的投影定义 ②投影之间的点、线、多边形和图像转换 ③集成以在 Matplotlib 中公开高级映射，具有简单和直观的界面 ④通过将 shapefile 读取与 Shapely 集成，实现强大的矢量数据处理能力。

二、数据的读取

        在我进行气象数据可视化的过程中，主要使用Python对.nc文件或.txt文件进行处理。下面我们分成两类进行讨论。

（一）.nc文件

        在开始之前，我们首先应该知道.nc文件中的数据类型是什么。在此，我们使用xarray库编写一个简单的程序，来确定其内容。我们以数据“hgt.mon.ltm.1991-2020”为例：

import xarray as xr
filename=r'hgt.mon.ltm.1991-2020.nc'#地形文件储存地址
f=xr.open_dataset(filename)#读取文件
print(f)#打印其属性

 （▲图2.1 查找数据类型）

（▲图2.2 图2.1运行结果）

        在这里我们可以看到，数据“hgt.mon.ltm.1991-2020”有五个纬度，分别为：level、lat、lon、time、nbnds。以其中“level”为例，他的坐标为1000，925，850…30，20，10，数据类型为float32。在绘制高度场时，我们主要使用的变量为“hgt”，该变量有四个纬度，分别为level、lat、lon和time。其储存的数据类型是float32。

        在了解了数据类型后，我们就可以对数据进行读取了。以ERA5的数据为例，我们使用netCDF4库中Dataset函数进行文件的读取，并把数据以数组的形式储存：

import netCDF4 as nc
fu= nc.Dataset(r"era5.u_component_of_wind.20180816.nc",'r')
fv = nc.Dataset(r"era5.v_component_of_wind.20180816.nc",'r')
fm=nc.Dataset(r"era5.mslp.20180816.nc",'r')
levels = fu.variables['level'][:]
lat = fu.variables['latitude'][:]
lon0 = fu.variables['longitude'][:]
time = nc.num2date(fu.variables['time'][:],
                   fu.variables['time'].units).data
months = np.array([i.month for i in time])
u = fu.variables['u'][0,0,:,:]
v = fv.variables['v'][0,0,:,:]
mslp = fm.variables['msl'][:]

（▲图2.3 数据的读取）

        此番操作后，数据便可以被直接调用：

（▲图2.4 Latitude数据示例）

（二）.txt文件

        相对于.nc文件，.txt文件的读取会简单很多。我们只需要用到一些Python基本的方法即可。读取的方法如下图所示：

def cc(j):
    f=open('20180630/aws_20180630'+j+'00.txt','r')
    lines = f.readlines()
    lat = []
    lon = []
    V = []
    for i in lines:
        a=i.split()
        V.append(eval(a[3]))
        lat.append(float(a[1]))
        lon.append(float(a[2]))
    V=np.array(V)

（▲图2.5 .txt文件的读取）

三、数据的处理

        在读取完数据后，有些情况下我们可以直接进行绘图，而有些情况下我们需要首先对数据进行处理。在本篇文章中，我将例举三个数据处理在课程中使用到的处理方法。包括风场数据的处理、降雨量数据的处理、以及极小值点叠加涡度条件的选取。

（一）风场数据的处理

        在绘制1000hPa平均风场时，我们需要对风向进行判断。在我们手中只有uwnd和vwnd两个风速文件的情况下，我们需要对其进行叠加处理。其基本原理是：利用for循环将u、v数据的每一纬度都变成数组以便于直接运算，转换完成后利用公式计算风矢量的大小和方向。接下来判断方向是否在(0,180)范围内，若在，为负；不在，为正。接下来重新将其封装为列表并返回，即可得到风矢量数据。

u = fu.variables['uwnd'][:]
v = fv.variables['vwnd'][:]
def uv(u,v):
    i=len(u)
    j=len(u[0])
    k=len(u[0][0])
    l=len(u[0][0][0])
    u=np.array(u)
    v=np.array(v)
    for a in range(i):
        u[a]=np.array(u[a])
        v[a]=np.array(v[a])
        for b in range(j):
            u[a][b]=np.array(u[a][b])
            v[a][b]=np.array(v[a][b])
            for c in range(k):
                u[a][b][c]=np.array(u[a][b][c])
                v[a][b][c]=np.array(v[a][b][c])
    uv=(u**2+v**2)**0.5
    wdir = 180.0 + np.arctan2(u, v)*180/np.pi
    for a in range(i):
        for b in range(j):
            for c in range(k):
                for d in range(l):
                    if 0 < wdir[a][b][c][d] < 180:
                        uv[a][b][c][d]=-uv[a][b][c][d]
    uv.tolist()
    for a in range(i):
        uv[a].tolist()
        for b in range(j):
            uv[a][b].tolist()
            for c in range(k):
                uv[a][b][c].tolist()
    return uv

（▲图3.1 风场数据的处理）

（二）降雨量数据的处理

        在处理降雨量数据时，我们主要使用scipy库中的griddata函数进行插值处理。在处理时有多种方法可供选择，在griddata函数的参数中可以更改。我们主要使用了‘linear’和‘cubic’两种方法。具体的使用方法是：首先确定插值使用的网格，用np.meshgrid函数生成网格点坐标矩阵，接下来直接使用griddata函数进行插值，之后转换为自己想要的数据格式即可。代码如下图所示：

from scipy.interpolate import interp1d
from scipy.interpolate import LinearNDInterpolator
from scipy.interpolate import griddata
def cc(j):
    f=open('20180630/aws_20180630'+j+'00.txt','r')
    lines = f.readlines()
    lat = []
    lon = []
    V = []
    for i in lines:
        a=i.split()
        V.append(eval(a[3]))
        lat.append(float(a[1]))
        lon.append(float(a[2]))
    V=np.array(V)

    X = np.arange(min(lon), max(lon),0.03)
    Y = np.arange(min(lat), max(lat),0.03)
    X, Y = np.meshgrid(X, Y)
    rain1 = griddata((lon,lat), V,(X,Y), method='linear')
    rain2 = griddata((lon,lat), V,(X,Y), method='cubic')
    for p in range(len(rain2)):
        for q in range(len(rain2[0])):
            if rain2[p][q]<0:
                rain2[p][q]=0
    return X,Y,rain1,rain2
X1,Y1,Z1,Z2=cc(h[0])

（▲图3.2 降雨量数据的处理）

（三）极小值点叠加涡度条件的选取

        此项数据处理是在判断台风路径时所需要的。其要求是：①提取台风区域MSLP场的局部极小值点对应的坐标；②对于每个候选坐标，判断其周围278KM范围内的最大涡度是否大于5×10-5 。若大于，则认为其为台风中心。在了解完处理原理后，我们首先使用scipy.signal中的argrelextrema函数提取全部MSLP场的极小值，提取后生成行和列的索引以列表的形势储存。接下来使用np.gradient函数（它相当于梯度算子）对u、v数据进行处理，接下来就可进行涡度判断。我们以此提取行和列中的每一个索引，用列表储存以它为中心2.5°范围内的所有涡度。储存后用其最大值与设定好的涡度（5×10-5）进行判断 。若其涡度大于5×10-5 ，则在外层列表中储存该涡度值；若小于则舍去。当储存完成后，我们提取最小的涡度值并返回其索引，储存在列表中，为之后将其与经纬度坐标相联系做准备。具体代码如下图所示：

I=[]
J=[]
for y in range(24):
    a = argrelextrema(mslp[y], np.less)
    rows, cols = a[0], a[1]
    duy, dux = np.gradient(u[y])
    dvy, dvx = np.gradient(v[y])
    vorticity = dvx-duy
    t=1.1131955*10**5
    w=5*10**(-5)
    W=[]
    v1=0
    for i,j in zip(rows,cols):
        W1=[]
        for k in range(-10,11):
            for l in range(-10,11):
                if 190<j+k<280 and 100<i+l<150:
                    W1.append(vorticity[i+l][j+k])
        if W1!= [] and max(W1)>t*w/4:
            W.append(mslp[y][i][j])
        else:
            v1+=1
    for p in range(len(mslp[0])):
        for q in range(len(mslp[0][0])):
            if mslp[y][p][q]==min(W):
                I.append(q)
                J.append(p)
LO=[]
LA=[]
for y in range(len(I)):
    LO.append(lon0[I[y]])
    LA.append(lat[J[y]])

（▲图3.3 极小值点叠加涡度条件）

四、图像的绘制

        在这一部分，我们着重介绍Python中Cartopy库的使用。我画图的方法主要分为三步：①确定绘图角度 ②确定绘图使用的函数 ③绘制图例信息。在开始绘制之前，我会首先创建画布：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure(figsize=[8,8])

        一般来说，绘图角度需要根据所展示的内容确定。例如：

        在绘制北半球500hPa高度场时，我们取北极极射投影。创建的图像如下图所示：

proj =ccrs.NorthPolarStereo(central_longitude=180)
ax = plt.subplot(121,projection=proj)

（▲图4.1 北极极射投影）

        在绘制1000hPa平均风场时，我们选择的是圆柱投影。创建的图像如下图所示：

ax = plt.subplot(211,projection=ccrs.PlateCarree(central_longitude=180))

（▲图4.2 圆柱投影）

        在绘制588线时，我们选取朗伯等角圆锥投影。创建的图像如下图所示：

ax = plt.subplot(211,projection=ccrs.LambertConformal(central_longitude=125))

（▲图4.3 朗伯等角圆锥投影）

        在选择完绘图角度后，我们要根据需要呈现的数据选择合适的绘图方法。我们可选择的绘图方法有等值线、填色图、色码图、流线图、雷达图以及直接绘制等等。其函数分别为：

#等值线
ac=ax.contour(Lon1,Lat1, hgt1,range(5879, 5881,30),colors=K[i-8],linewidths=0.5,linestyles='-',transform=ccrs.PlateCarree())
#填色图
ec = ax.contourf(Lon2,Lat2, hgt2/10,cmap = 'jet',levels=levels,transform=ccrs.PlateCarree())
#色码图
ec=ax.pcolormesh(X, Y, Z2,norm=norm,cmap=cmap,transform=ccrs.PlateCarree())
#流线图
ax.streamplot(Lon1,Lat1,u1,v1,transform=ccrs.PlateCarree(),color='k',linewidth=0.75,arrowsize=0.75,density=2)
#雷达图
ax.scatter(361,91, 1,c='b',alpha=1,marker='o',linewidths=2,transform=ccrs.PlateCarree())
#直接绘制（台风移动路径）
ax.plot(LO,LA,color='green')
ax.plot(LO,LA,'o',color='green',markersize='2')

                                        （a）等值线                                              （b）填色图

                      

                    （c）色码图                                                             （d）流线图

      

                                （e）雷达图                                                    （f）台风路径

（▲图4.4 六种绘图方法）

        在完成数据可视化绘制后，我们还要添加响应的地理信息、经纬网并对图片进行修饰。地理信息的添加如下图所示，我们主要对陆地、海洋、海岸线、河流、湖泊进行描绘。在其函数中可以设施绘制地图的分辨率以及修改某些性质，如下所示：

import cartopy.feature as cfeature
ax.add_feature(cfeature.LAND,color='white')
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN,color='cyan')
ax.add_feature(cfeature.LAKES,color='cyan')
ax.add_feature(cfeature.RIVERS,color='k')
LAKES_border=cfeature.NaturalEarthFeature('physical','lakes','10m',edgecolor='black',facecolor='never')
ax.add_feature(LAKES_border,linewidth=0.6)

        其中值得注意的是，我使用了NaturalEarthFeature函数设置湖泊边界为黑色并添加进入绘制当中。所得效果如下图所示：

（▲图4.5 设置地理信息）

        接下来开始绘制经纬网。在这里我们使用gridlines函数绘制。代码如下图所示：

gl = ax.gridlines(draw_labels=True, color='gray',alpha=0.5,linestyle='--',xlocs=np.arange(-180,181,5),ylocs=np.arange(-90,91,5),x_inline=False,y_inline=False)

        其中，我设置经纬网的颜色为灰色，线状为虚线并设置经纬网的步长为5°。

        最后，我们可以添加图例信息进行图片的修饰。

        在添加标题时，我们可以使用title函数设置每幅图的标题，使用subtitle函数设置大标题：

ax.set_title("8-10月",fontsize=3,x=0.135,y=0.8,color='k')
plt.suptitle("图5-4-4   500hPa西太平洋副高脊的月平均位置",fontsize=4)

        若图中存在多个线条，我们可以单独设置标签。在这里，如果有些图无法通过函数直接添加标签，我的方法是在经纬度范围外画短线添加标签。这样既可以实现要求，又不会呈现在图片上。具体方法如下图所示：

labels=['0.1mm','10mm','25mm','50mm','100mm','250mm']
color2=['#FFFFFF00','#a6f28e','#3dba3d','#61b8ff','#0000e1','#fa00fa','#800040']
for i in range(len(labels)):
    plt.plot(x,y,color=color2[i+1],label=labels[i])
plt.legend()

（▲图4.6 设置标签）

        除了添加标签，我们也可以通过直接在图片上添加文字来实现线条的区分。代码如下所示：

ax.set_title("5-8月",fontsize=3,x=0.13,y=0.8,color='k')
ax.text(-2500000, 0, '5月(——)', color='red',fontsize=2,weight='bold')
ax.text(-2500000, -200000, '6月(——)', color='green',fontsize=2,weight='bold')
ax.text(-2500000, -400000, '7月(——)', color='blue',fontsize=2,weight='bold')
ax.text(-2500000, -600000, '8月(——)', color='purple',fontsize=2,weight='bold')

（▲图4.7 添加文字）

        我们还可以设置图片的范围和经纬度坐标。使用set_extent函数设置图片范围：

ax.set_extent([100,160,5,40],crs=ccrs.PlateCarree())

         在设置坐标时，使用LONGITUDE_FORMATTER, LATITUDE_FORMATTER函数直接进行转化：

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
gl.xformatter = LONGITUDE_FORMATTER  #使横坐标转化为经纬度格式            
gl.yformatter = LATITUDE_FORMATTER   #使纵坐标转化为经纬度格式

        也可以自行添加坐标：

ax.set_xticks([-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180])#指定要显示的经纬度 
ax.set_yticks([-60,-30,0,30,60])
ax.xaxis.set_major_formatter(LongitudeFormatter())#刻度格式转换为经纬度样式 
ax.yaxis.set_major_formatter(LatitudeFormatter())

        或者直接进行字符串形式的填写：

ax.set_yticks(np.linspace(100, 1000, 10),[1000,900,800,700,600,500,400,300,200,100][::-1])
ax.set_xlim(-30, 60)
ax.set_xticks(range(-30,61,10),['30°','20°','10°S', 'EQ','10°','20°', '30°','40°','50°','60°N'])
ax.set_xlim([-30,60])

        最后我们可以添加colorbar，并在设置colorbar时对cmap和norm进行约束。

        Colormap (cmap)是 Python 中一个用于可视化的重要工具。它是一种颜色映射方式，将数据映射到颜色空间中，使得数据的不同值能够以不同的颜色显示。在数据分析、科学可视化和机器学习等领域中，cmap 通常被用来表示不同的数据范用或者数据类型。

        下面举两个例子：

        ①自定义colorbar

color_list=['#FFFFFF00','#a6f28e','#3dba3d','#61b8ff','#0000e1','#fa00fa']
levels=[0.0,0.1,10.0,25.0,50.0,100.0,250.0]
cmap=mcolors.ListedColormap(color_list)
cmap.set_over('#800040')
norm = mcolors.BoundaryNorm(levels, cmap.N)
cb=fig.colorbar(ec,extend='max')
ax2 = cb.ax
ax2.set_ylim(0,250)
ax2.set_yticks([0,0.1,10,25,50,100,250])

        ②混合色彩的colorbar

cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=0, vmax=30)
ac=ax.streamplot(Lon,Lat,u3,v3,transform=ccrs.PlateCarree(),color=hs3,norm=norm,cmap=cmap,linewidth=1,arrowsize=0.75,density=4)
fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), label='(m/s)',shrink=0.75)

          

                 （a）自定义colorbar                                      （b）混合色彩的colorbar

        在所有工作做完后，我们可以设置文件名并保存： 

plt.savefig("200hPa水平风场流线图.png",dpi=800)
plt.show()

        一幅气象信息图就可呈现在各位眼前了。

参考文献

[1]利用Python(netCDF4库)读取.nc文件（NetCDF气象数据文件）的基本操作

利用Python(netCDF4库)读取.nc文件（NetCDF气象数据文件）的基本操作_<class 'netcdf4._netcdf4.variable'>-CSDN博客

[2]Matplotlib: Visualization with Python Matplotlib — Visualization with Python

[3]NumPy官网    NumPy

[4]Scipy官网    SciPy

[5]Cartopy官网   Cartopy · PyPI