在交错网格上使用数据时计算以网格为中心的值
Calculating grid-centred values when using data on a staggered grid
我正在与 MITgcm 合作进行一些模拟，特别是使用内波模型；我得到了包含结果的 .nc 文件，但一些变量的坐标并不完全相同。我会解释一下：我想计算出速度的分量，但是由于一些我不完全理解的数字原因，水平速度坐标在单元格的左侧，垂直坐标在单元格的底部.要使用速度数据进行操作，我需要统一所有坐标的参考。
我想过做这样的事情
u (i,j,k) = u(i,j,k) + u(i+1,j,k)
v (i,j,k) = v(i,j,k) + v(i,j+1,k)
所以我的坐标都在单元格的中心和同一个参考中。
我不知道该怎么做，使用 python，编辑 NetCDF 文件。我可能会很高兴提取所有 u 和 v 数据，像我说的那样编辑并用这两个变量创建一个新的 NetCDF 文件。
这可能吗？我怎样才能做到这一点？
编辑：添加了 ncdump 信息
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 ;
}
该模型使用 staggered grid，其中 u 速度在 west/east 网格面上求解，而 v 速度在 north/south 面上求解。
你是对的，现在你需要 post- 处理组件，以便 u- 和 v- 分别位于每个网格单元的中心。
让我们将 nx 定义为 x 维度（即求解 u 分量的位置）中的网格单元数，将 ny 定义为 x 维度中的网格单元数y 维度（即求解 v 分量的位置）。 nz为垂直模型层数。
那么 u 的尺寸为 nx+1 x ny x nz 并且 v 的尺寸为 nx x ny+1 x nz。然后将 u 和 v 放入每个单元格的中心是一个简单的平均值：
u_center = 0.5 * (u[0:nx,:,:] + u[1:nx+1,:,:]) # now has dimensions [nx,ny,nz])
v_center = 0.5 * (v[:,0:ny,:] + v[:,1:ny+1,:]) # now has dimensions [nx,ny,nz])
import netCDF4
import numpy as np
ncfile = netCDF4.Dataset('/path/to/file/foo.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)
u_out = ncfile_out.createVariable('u_center', 'f4', ('longitude', 'latitude', 'level')
v_out = ncfile_out.createVariable('v_center', 'f4', ('longitude', 'latitude', 'level')
u_out[:,:,:] = u_center[:,:,:]
v_out[:,:,:] = v_center[:,:,:]
ncfile_out.close()
在命令行中使用 cdo 使用 2 点平均并结合使用移位函数和整体均值函数来执行此操作怎么样？
cdo selvar,U mitfile.nc u.nc # select the velocity fields
cdo selvar.v mitfile.nc v.nc
# shift to the right/up and average both memberss
cdo ensmean -shiftx,1 u.nc u.nc ucen.nc
cdo ensmean -shifty,1 v.nc v.nc vcen.nc
# cat the files and calculate the wind:
cdo cat ucen.nc vcen.nc uv.nc
cdo expr,"wind=sqrt(U*U+V*V)" uv.nc wind.nc
ps：如果你的字段是全局的，那么你想在经度方向上做一个循环移位：
cdo ensmean -shiftx,1,cyclic u.nc u.nc ucen.nc
我正在与 MITgcm 合作进行一些模拟，特别是使用内波模型；我得到了包含结果的 .nc 文件，但一些变量的坐标并不完全相同。我会解释一下：我想计算出速度的分量，但是由于一些我不完全理解的数字原因，水平速度坐标在单元格的左侧，垂直坐标在单元格的底部.要使用速度数据进行操作，我需要统一所有坐标的参考。
我想过做这样的事情
u (i,j,k) = u(i,j,k) + u(i+1,j,k)
v (i,j,k) = v(i,j,k) + v(i,j+1,k)
所以我的坐标都在单元格的中心和同一个参考中。
我不知道该怎么做，使用 python，编辑 NetCDF 文件。我可能会很高兴提取所有 u 和 v 数据，像我说的那样编辑并用这两个变量创建一个新的 NetCDF 文件。
这可能吗？我怎样才能做到这一点？
编辑：添加了 ncdump 信息
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 ;
}
该模型使用 staggered grid，其中 u 速度在 west/east 网格面上求解，而 v 速度在 north/south 面上求解。
你是对的，现在你需要 post- 处理组件，以便 u- 和 v- 分别位于每个网格单元的中心。
让我们将 nx 定义为 x 维度（即求解 u 分量的位置）中的网格单元数，将 ny 定义为 x 维度中的网格单元数y 维度（即求解 v 分量的位置）。 nz为垂直模型层数。
那么 u 的尺寸为 nx+1 x ny x nz 并且 v 的尺寸为 nx x ny+1 x nz。然后将 u 和 v 放入每个单元格的中心是一个简单的平均值：
u_center = 0.5 * (u[0:nx,:,:] + u[1:nx+1,:,:]) # now has dimensions [nx,ny,nz])
v_center = 0.5 * (v[:,0:ny,:] + v[:,1:ny+1,:]) # now has dimensions [nx,ny,nz])
import netCDF4
import numpy as np
ncfile = netCDF4.Dataset('/path/to/file/foo.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)
u_out = ncfile_out.createVariable('u_center', 'f4', ('longitude', 'latitude', 'level')
v_out = ncfile_out.createVariable('v_center', 'f4', ('longitude', 'latitude', 'level')
u_out[:,:,:] = u_center[:,:,:]
v_out[:,:,:] = v_center[:,:,:]
ncfile_out.close()
在命令行中使用 cdo 使用 2 点平均并结合使用移位函数和整体均值函数来执行此操作怎么样？
cdo selvar,U mitfile.nc u.nc # select the velocity fields
cdo selvar.v mitfile.nc v.nc
# shift to the right/up and average both memberss
cdo ensmean -shiftx,1 u.nc u.nc ucen.nc
cdo ensmean -shifty,1 v.nc v.nc vcen.nc
# cat the files and calculate the wind:
cdo cat ucen.nc vcen.nc uv.nc
cdo expr,"wind=sqrt(U*U+V*V)" uv.nc wind.nc
ps：如果你的字段是全局的，那么你想在经度方向上做一个循环移位：
cdo ensmean -shiftx,1,cyclic u.nc u.nc ucen.nc