不同维数的python netCDF4变量乘法

Question

我有第一个变量h有维(111,141)还有一个变量cs_w有维(51，)。基本上，我的数据是ROMS历史记录输出数据。现在我想用cs_w乘以h，最后的结果应该是维数(51,111,141)。在这里，我被困住了，无法继续前进。下面是我的代码

import numpy as np

import matplotlib.pyplot as plt

import netCDF4 as nc

f_in = nc.Dataset('ocean_his_0010.nc', "r")



h = f_in.variables['h']

cs_w = f_in.variables['Cs_w']


z=[[],[],[]]
for i in range(len(h[0])):
    for j in range(len(h[1])):
        for k in range(len(cs_w)):
            z[i][j][k] = h[i][j]*cs_w[k]

下面是我想使用的两个变量的描述。

Out88: float64 Cs_w(s_w) long_name: s坐标拉伸曲线在W点valid_min：-1.0 valid_max: 0.0字段: Cs_w，标量无限维:当前形状= (51 )填充，默认9.969209968386869e+36使用_FillValue

H Out89: float64 h(eta_rho，xi_rho) long_name: RHO点单位测深:米网格:网格位置:面坐标: lon_rho lat_rho字段:沐浴，标量无限维:当前形状= (111,141)填充，默认_FillValue使用9.969209968386869e+36

Below is the ncdump ocean_his_0010.nc

netcdf ocean_his_0010 {
dimensions:
        xi_rho = 141 ;
        xi_u = 140 ;
        xi_v = 141 ;
        xi_psi = 140 ;
        eta_rho = 111 ;
        eta_u = 111 ;
        eta_v = 110 ;
        eta_psi = 110 ;
        N = 50 ;
        s_rho = 50 ;
        s_w = 51 ;
        tracer = 2 ;
        boundary = 4 ;
        ocean_time = UNLIMITED ; // (360 currently)

        double Cs_w(s_w) ;
                Cs_w:long_name = "S-coordinate stretching curves at W-points" ;
                Cs_w:valid_min = -1. ;
                Cs_w:valid_max = 0. ;
                Cs_w:field = "Cs_w, scalar" ;
        double h(eta_rho, xi_rho) ;
                h:long_name = "bathymetry at RHO-points" ;
                h:units = "meter" ;
                h:grid = "grid" ;
                h:location = "face" ;
                h:coordinates = "lon_rho lat_rho" ;
                h:field = "bath, scalar" ;

Answer 1

您不能仅仅将一个(多维)列表定义为z=[[],[],[]]并开始以您正在尝试的方式填充它，它首先需要适当的大小。参见这个问题/答案，它处理相同的问题。

Numpy通常更方便地处理nD-数组，您的z数组可以简单地定义为：

z = np.zeros((51,111,141))

并使用类似嵌套循环的内容或使用向量化指令填充，例如：

for k in range(51):
    z[k,:,:] = cs_w[k] * h[:,:]

甚至是完全自动的(甚至不必事先定义z )：

import numpy as np

h    = np.zeros((111,141))
cs_w = np.zeros(51)

z = cs_w[:,np.newaxis,np.newaxis] * h

使用这些矢量化操作通常比手动编写循环要快得多。