Python遍历NetCDF时间步计算后合并为三维数组问题求助

我有一个存储了风场u、v分量的NetCDF文件，包含7305个时间步的日均值数据。需要逐时间步完成指定计算后，将每个时间步的输出结果合并为一个总数组，每个子数组对应一天的计算结果。原始NetCDF文件的纬度维度为101、经度维度为129，因此最终输出的涡度（vort）数组维度应为(7305,101,129)。目前编写的代码单独存储单日计算结果时可正常运行，但使用append方法合并数组时出现错误。

from matplotlib import pyplot
from matplotlib.cm import get_cmap
from __future__ import print_function
from netCDF4 import Dataset,num2date,date2num
from matplotlib.colors import from_levels_and_colors
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from metpy.units import units
from datetime import datetime
from metpy.plots import StationPlot
#
import metpy.calc as mpcalc
import xarray as xr
import cartopy.crs as ccrs
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import datetime
import cartopy.feature as cfeature
import mygrads as mg
#
#
root_dir = '/users/pr007/mkaryp/'
nc = Dataset(root_dir+'ERA5_ensmean_daymean.nc')
#
lon=nc.variables['longitude'][:]
lat=nc.variables['latitude'][:]
#
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
#
u = []
v = []
vort = np.array([])
#
for i in range(7305):
    v = np.array(nc.variables['v'][i,:,:])                    
    u = np.array(nc.variables['u'][i,:,:]) 
    #   
    v = units.Quantity(v, 'm/s')
    u = units.Quantity(u, 'm/s')
    #
    # Compute dx and dy spacing for use in divergence calculation
    vort1 = np.array(mpcalc.vorticity(u, v, dx=dx, dy=dy))
    vort = np.append([vort, vort1])
    #
    print(i)

提供了仅包含10个时间步的测试用NetCDF子样本文件用于验证。

核心错误点 #

1. np.append语法错误：该函数要求至少传入两个参数（原数组、待追加的内容），你将两个数组打包成一个列表作为唯一参数传入，不符合函数调用要求，直接触发报错。
2. 数组维度不匹配：初始化的vort是空一维数组，直接追加二维的单日涡度结果会将数据扁平化，最终无法得到(7305,101,129)的三维结构，且反复append会频繁重新申请内存，运行效率极低。

推荐修复方案（预分配数组，性能最优） #

直接提前创建符合目标维度的空数组，循环中直接对对应时间步赋值即可，无需使用append：

• 将vort = np.array([])替换为：

vort = np.empty((7305, len(lat), len(lon)))

• 将循环内的vort = np.append([vort, vort1])替换为：

vort[i,:,:] = vort1

修改后运行即可直接得到维度为(7305,101,129)的涡度数组。

可选优化方案（无循环实现） #

使用xarray直接读取数据配合metpy的广播计算能力，无需手动写循环处理时间步，代码更简洁：

ds = xr.open_dataset(root_dir+'ERA5_ensmean_daymean.nc')
u = ds['u'] * units.m/s
v = ds['v'] * units.m/s
vort = mpcalc.vorticity(u, v, dx=dx, dy=dy).to_numpy()

内容的提问来源于stack exchange，提问作者Maria Chara Karypidou