2016-06-06 42 views
0

我正在使用MITgcm进行一些模拟,特别是使用内部波浪模型;我用我的结果得到了.nc文件,但是一些变量不在完全相同的坐标中。我会解释一下自己:我想研究速度的组成部分,但由于某些我不完全了解的数值原因,水平速度坐标位于单元格的左侧和单元格的机器人侧的垂直坐标。要使用速度数据进行操作,我需要统一所有坐标的参考。编辑netcdf变量并使用python写出一个新的netcdf文件

要做到这一点,我有这样的代码:

import netCDF4 
import numpy as np 

ncfile = netCDF4.Dataset('state.global.nc', 'r') 
u = ncfile.variables['U'][:,:,:] # nx+1 x ny x nz 
v = ncfile.variables['V'][:,:,:] # nx x ny+1 x nz 

nx = np.shape(u)[0] - 1 
ny = np.shape(v)[1] - 1 
nz = np.shape(u)[2] 

u_center = 0.5 * (u[0:nx,:,:] + u[1:nx+1,:,:]) 
v_center = 0.5 * (v[:,0:ny,:] + v[:,1:ny+1,:]) 

# Write out u_center and v_center into a new netCDF file 
ncfile_out = netCDF4.Dataset('./output.nc', 'w') 
ncfile_out.createDimension('longitude', nx) 
ncfile_out.createDimension('latitude', ny) 
ncfile_out.createDimension('level', nz) 
ncfile_out.createDimension('time', None) 
u_out = ncfile_out.createVariable('u_center', 'f4', ('time', 'longitude', 'latitude', 'level')) 
v_out = ncfile_out.createVariable('v_center', 'f4', ('time', 'longitude', 'latitude', 'level')) 
time = ncfile_out.createVariable('Time', 'i4', 'time') 
u_out[:,:,:] = u_center[:,:,:] 
v_out[:,:,:] = v_center[:,:,:] 
ncfile_out.close() 

不过,这并不编译,它告诉我一个错误的线24和25(U_OUT [:,:,:] ...);正好说'IndexError:数据数组的大小不符合切片'。我尝试通过u_out [:,:,:,]更改u_out [:,:,:]等。我不知道我的错误是什么。

欲了解更多信息,我在这里贴了我原来的netCDF文件信息:

netcdf state.global { 
dimensions: 
    T = UNLIMITED ; // (10001 currently) 
    Xp1 = 61 ; 
    Y = 1 ; 
    Z = 20 ; 
    X = 60 ; 
    Yp1 = 2 ; 
    Zl = 20 ; 
variables: 
    double Xp1(Xp1) ; 
     Xp1:long_name = "X-Coordinate of cell corner" ; 
     Xp1:units = "meters" ; 
    double Y(Y) ; 
     Y:long_name = "Y-Coordinate of cell center" ; 
     Y:units = "meters" ; 
    double Z(Z) ; 
     Z:long_name = "vertical coordinate of cell center" ; 
     Z:units = "meters" ; 
     Z:positive = "up" ; 
    double X(X) ; 
     X:long_name = "X-coordinate of cell center" ; 
     X:units = "meters" ; 
    double Yp1(Yp1) ; 
     Yp1:long_name = "Y-Coordinate of cell corner" ; 
     Yp1:units = "meters" ; 
    double Zl(Zl) ; 
     Zl:long_name = "vertical coordinate of upper cell interface" ; 
     Zl:units = "meters" ; 
     Zl:positive = "up" ; 
    double T(T) ; 
     T:long_name = "model_time" ; 
     T:units = "s" ; 
    int iter(T) ; 
     iter:long_name = "iteration_count" ; 
    double U(T, Z, Y, Xp1) ; 
     U:units = "m/s" ; 
     U:coordinates = "XU YU RC iter" ; 
    double V(T, Z, Yp1, X) ; 
     V:units = "m/s" ; 
     V:coordinates = "XV YV RC iter" ; 
    double Temp(T, Z, Y, X) ; 
     Temp:units = "degC" ; 
     Temp:long_name = "potential_temperature" ; 
     Temp:coordinates = "XC YC RC iter" ; 
    double S(T, Z, Y, X) ; 
     S:long_name = "salinity" ; 
     S:coordinates = "XC YC RC iter" ; 
    double Eta(T, Y, X) ; 
     Eta:long_name = "free-surface_r-anomaly" ; 
     Eta:units = "m" ; 
     Eta:coordinates = "XC YC iter" ; 
    double W(T, Zl, Y, X) ; 
     W:units = "m/s" ; 
     W:coordinates = "XC YC RC iter" ; 

// global attributes: 
     :MITgcm_version = "****************" ; 
     :build_user = "************" ; 
     :build_host = "**************" ; 
     :build_date = "*******************" ; 
     :MITgcm_URL = "***************" ; 
     :MITgcm_tag_id = "*******************" ; 
     :MITgcm_mnc_ver = 0.9 ; 
     :sNx = 30 ; 
     :sNy = 1 ; 
     :OLx = 2 ; 
     :OLy = 2 ; 
     :nSx = 2 ; 
     :nSy = 1 ; 
     :nPx = 1 ; 
     :nPy = 1 ; 
     :Nx = 60 ; 
     :Ny = 1 ; 
     :Nr = 20 ; 
} 

我在速度变量只是感兴趣。

非常感谢。

回答

1

您将u_out创建为具有四个维度('time','longitude','latitude','level')的netCDF变量,但您使用3D切片([:,:,])进行赋值数据。