Berechnen von Gitter-zentrierten Werten bei Verwendung von Daten auf einem gestaffelten Gitter

Question

Ich arbeite mit MITgcm, um einige Simulationen zu machen, speziell um mit internen Wellenmodellen zu arbeiten; Ich habe .NC-Dateien mit meinen Ergebnissen, aber einige Variablen sind nicht genau in den gleichen Koordinaten. Ich erkläre mich selbst: Ich möchte die Komponenten der Geschwindigkeit berechnen, aber aus einigen numerischen Gründen, die ich nicht vollständig verstehe, sind horizontale Geschwindigkeitskoordinaten auf der linken Seite der Zellen und vertikale Koordinaten auf der Bot-Seite der Zellen . Um mit Geschwindigkeitsdaten arbeiten zu können, muss ich die Referenz aller Koordinaten vereinheitlichen.

Ich dachte an so etwas wie diese

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) 

Damit werde ich meine Koordinaten alle in der Mitte der Zellen und in der gleichen Referenz.

Ich weiß nicht, wie es geht, mit Python, Bearbeiten einer NetCDF-Datei. Ich könnte glücklich sein, nur alle u und v Daten zu extrahieren, redigierend, wie ich sagte und eine neue NetCDF-Datei mit nur diesen zwei Variablen schaffend.

Ist das möglich? Wie kann ich das machen?

Edit: Hinzugefügt ncdump Informationen

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

Answer 1

Das Modell wird mit einem staggered grid wo u-Geschwindigkeit auf der West/Ost Gitter gelöst zugewandt ist, während V-Geschwindigkeit auf die Nord/Süd-Gesichter gelöst.

Sie haben recht, dass Sie nun die Komponenten nachbearbeiten müssen, so dass sich u- und v- jeweils in der Mitte jeder Gitterzelle befinden.

Definieren wir nx als die Anzahl der Gitterzellen in der x-Dimension (dh wo die u-Komponente gelöst ist) und als Anzahl der Gitterzellen in der y-Dimension (dh wo die v- Komponente ist gelöst). nz ist die Anzahl der vertikalen Modellschichten.

Dann hat u Abmessungen nx+1 x ny x nz und v Dimensionen hat nx x ny+1 x nz. Es ist ein einfacher Durchschnitt dann u und v in die Zentren jeder Zelle zu erhalten:

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

Answer 2

Wie wäre dies von der Kommandozeile mit cdo macht einen 2-Punkt-Durchschnitt mit Verwenden Sie die Shift-Funktion in Kombination mit der Ensemble-Mean-Funktion?

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: zu verwenden, wenn Sie Ihre Felder global sind, dann wollen Sie eine zyklische Verschiebung in der Längenrichtung tun:

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