Calculating grid-centred values when using data on a staggered grid

Question

I am working with MITgcm to make some simulations, specifically to work with internal waves models; I got .nc files with my results but some variables are not in exactly the same coordinates. I'll explain myself: I want to work out the components of velocity but, for some numerical reasons that I don't fully understand, horizontal velocity coordinates are in the left side of the cells and vertical coordinates on the bot side of the cells. To operate with velocity data I need to unify the reference of all coordinates.

I thought about doing something like this

 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)

So I'll have my coordinates all in the center of the cells and in the same reference.

I don't know how to do it, using python, editing a NetCDF file. I could be happy just extracting all u and v data, editing like I said and creating a new NetCDF file with just this two variables.

Is that possible? How can I do that?

Edit: Added ncdump information

   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 ;
}

N1B4 · Accepted Answer · 2016-06-02T21:37:30.340

3

That model is using a staggered grid where u-velocity is solved on the west/east grid faces, while v-velocity is solved on the north/south faces.

You are right that now you need to post-process the components so that u- and v- are each placed in the center of each grid cell.

Let's define nx to be the number of grid cells in the x-dimension (i.e. where the u-component is solved on) and ny to be the number of grid cells in the y-dimension (i.e. where the v-component is solved on). nz is the number of vertical model layers.

Then u has dimensions nx+1 x ny x nz and v has dimensions nx x ny+1 x nz. It's a simple average then to get u and v into the centers of each cell:

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()

edited Jun 02 '16 at 21:37

answered Jun 02 '16 at 14:58

N1B4

3,377
1
21
24

It went well until `u_center =0.5*...`, Then it say this error: `ValueError: operands could not be broadcast together with shapes (9999,20,1,61) (10000,20,1,61) ` – Alejandro Jiménez Rico Jun 02 '16 at 18:59
I've edited the answer to correct the indexing error (`nx-1` to `nx`, and `ny-1` to `ny`). I forgot that `nx` and `ny` were previously defined with the minus 1 attached to it, so we don't need that in the centering lines. – N1B4 Jun 02 '16 at 19:22
Great! Now it works. With this, my netcdf file is edited now? Or I can export my new relocated variables? If not, How can I do it? – Alejandro Jiménez Rico Jun 02 '16 at 20:57
I have now included a sample write-out of a netCDF file. I recommend reading the python-netCDF4 documentation for a better understanding of how to do this for your case. If you believe I have answered your question, you may also opt to mark this answer as solved. – N1B4 Jun 02 '16 at 21:39
It tells an error at line 22. `Value error: total size of new array must be unchanged` – Alejandro Jiménez Rico Jun 03 '16 at 06:35
Your data may be (lat x lon x level) instead of (lon x lat x level)...or maybe it's 4D with time in there. You'll need to provide more details on what `u` and `v` look like and/or adjust the code above to fit your data appropriately. – N1B4 Jun 03 '16 at 11:15
Yeah, my system is time-dependent. I'm sorry to forget mentioning it. How can I complete what you did with time dependence? – Alejandro Jiménez Rico Jun 03 '16 at 18:22
I edited the question providing more details about 'u' and 'v'. With this, could you help me to finish the code? – Alejandro Jiménez Rico Jun 06 '16 at 21:42

ClimateUnboxed · Answer 2 · 2017-10-18T15:20:43.840

How about doing this from the command line with cdo using a 2-point average using the shift function combined with the ensemble mean function?

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: if your fields are global, then you want to do a cyclic shift in the longitude direction using:

cdo ensmean -shiftx,1,cyclic u.nc u.nc ucen.nc

Calculating grid-centred values when using data on a staggered grid

2 Answers2