format of input files for intlevel3d
Added by Fred W almost 5 years ago
Hello all,
I am trying to use cdo intlevel3d/intlevelx3d to regrid data in the vertical from a hybrid s-coordinate system (NEMO ocean model) to standard depth-levels (i.e. z-coordinates).
I have prepared all the input files as best I can, but it's not working as intended. Here's where I got to so far:
cdo --verbose intlevel3d,MODEL_gdept.nc model_output_original.nc MODEL_std_vertgrid.nc model_output_regrid.nc
The file MODEL_gdept.nc contains the gdept variable (30 levels) which is a 3D description of the depth at each grid point.
netcdf MODEL_gdept {
dimensions:
depth = 30 ;
lat = 172 ;
lon = 244 ;
variables:
float depth(depth) ;
depth:standard_name = "model_level_number" ;
depth:units = "m" ;
depth:positive = "down" ;
depth:long_name = "Vertical T levels" ;
float gdept(depth, lat, lon) ;
gdept:_FillValue = 9.96921e+36f ;
gdept:coordinates = "depth lat lon" ;
gdept:long_name = "Ocean depth (masked)" ;
gdept:units = "m" ;
float lat(lat, lon) ;
lat:standard_name = "latitude" ;
lat:units = "degrees_north" ;
float lon(lat, lon) ;
lon:standard_name = "longitude" ;
lon:units = "degrees_east" ;
// global attributes:
...
}
The gdept variable is masked, and contains NaNf in dry grid cells. Values on the last level are all NaNf (no data values are ever present on the last level either).
The model_output_original.nc file is sample model output with 24 records of temperature (244x172x30), salinity (244x172x30) and sea surface height (244x172). The horizontal and vertical grid matches that in MODEL_gdept.nc.
Note the coordinate variables are 2D fields, as it's a curvilinear grid. The horizontal regridding is something I have yet to tackle in step 2 after getting the vertical regridding right.
netcdf model_output_original {
dimensions:
depth = 30 ;
time = UNLIMITED ; // (24 currently)
lat = 172 ;
lon = 244 ;
variables:
float depth(depth) ;
depth:axis = "Z" ;
depth:standard_name = "model_level_number" ;
depth:units = "m" ;
depth:positive = "down" ;
short sossheig(time, lat, lon) ;
sossheig:_FillValue = -32768s ;
sossheig:missing_value = -32768s ;
sossheig:standard_name = "sea_surface_height_above_geoid" ;
sossheig:units = "m" ;
sossheig:coordinates = "lon lat time" ;
double time(time) ;
time:standard_name = "time" ;
time:units = "seconds since 2019-04-20 00:00:00" ;
time:calendar = "gregorian" ;
short vosaline(time, depth, lat, lon) ;
vosaline:_FillValue = -32768s ;
vosaline:missing_value = -32768s ;
vosaline:standard_name = "sea_water_salinity" ;
vosaline:units = "1e-3" ;
vosaline:coordinates = "lon lat depth time" ;
short votemper(time, depth, lat, lon) ;
votemper:_FillValue = -32768s ;
votemper:missing_value = -32768s ;
votemper:standard_name = "sea_water_potential_temperature" ;
votemper:units = "C" ;
votemper:coordinates = "lon lat depth time" ;
float lat(lat, lon) ;
lat:standard_name = "latitude" ;
lat:units = "degrees_north" ;
lat:long_name = "Latitude" ;
float lon(lat, lon) ;
lon:standard_name = "longitude" ;
lon:units = "degrees_east" ;
lon:long_name = "Longitude" ;
// global attributes:
...
}
To generate MODEL_std_vertgrid.nc I wrote a python program that creates a file with the same horizontal grid (244x172, but 1D coordinate variables as vectors), but a new vertical grid containing the WOA13 standard depth levels.
netcdf MODEL_std_vertgrid {
dimensions:
depth = 102 ;
time = UNLIMITED ; // (1 currently)
lat = 172 ;
lon = 244 ;
variables:
int depth(depth) ;
depth:long_name = "depth" ;
depth:standard_name = "depth" ;
depth:units = "meters" ;
depth:positive = "down" ;
int gdept(time, depth, lat, lon) ;
gdept:long_name = "gdept" ;
gdept:standard_name = "gdept" ;
gdept:units = "meters" ;
gdept:positive = "down" ;
float lat(lat) ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:units = "degrees_north" ;
float lon(lon) ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:units = "degrees_east" ;
double time(time) ;
time:long_name = "time" ;
time:standard_name = "time" ;
time:units = "hours since 1970-01-01 00:00:00.0" ;
time:calendar = "gregorian" ;
// global attributes:
...
}
Now here is the output from the cdo command line mentioned above:
OpenMP: num_procs=40 max_threads=1 cdo intlevel3d: 30 records input 3d vertical height cdo intlevel3d: 102 records target 3d vertical height and gridsize 41968 cdo intlevel3d: Input vertical coordinate has no missing values. cdo intlevel3d: Output vertical coordinate has no missing values. cdo intlevel3d (Warning): Non monotonic zaxis! cdo intlevel3d: lev1 0: nan cdo intlevel3d: lev1 1: nan cdo intlevel3d: lev1 2: nan cdo intlevel3d: lev1 3: nan cdo intlevel3d: lev1 4: nan cdo intlevel3d: lev1 5: nan cdo intlevel3d: lev1 6: nan cdo intlevel3d: lev1 7: nan cdo intlevel3d: lev1 8: nan cdo intlevel3d: lev1 9: nan cdo intlevel3d: lev1 10: nan cdo intlevel3d: lev1 11: nan cdo intlevel3d: lev1 12: nan cdo intlevel3d: lev1 13: nan cdo intlevel3d: lev1 14: nan cdo intlevel3d: lev1 15: nan cdo intlevel3d: lev1 16: nan cdo intlevel3d: lev1 17: nan cdo intlevel3d: lev1 18: nan cdo intlevel3d: lev1 19: nan cdo intlevel3d: lev1 20: nan cdo intlevel3d: lev1 21: nan cdo intlevel3d: lev1 22: nan cdo intlevel3d: lev1 23: nan cdo intlevel3d: lev1 24: nan cdo intlevel3d: lev1 25: nan cdo intlevel3d: lev1 26: nan cdo intlevel3d: lev1 27: nan cdo intlevel3d: lev1 28: nan cdo intlevel3d: lev1 29: nan cdo intlevel3d: lev1 30: 0 cdo intlevel3d: lev1 31: 0 cdo intlevel3d (Abort): Level 0 not found!
Several things strike me as odd:
The first thing is that cdo reports "Input vertical coordinate has no missing values." which is clearly untrue. The gdept variable in MODEL_gdept.nc has missing values shown as NaNf by ncdump.
The second thing which is weird is the warning "Non monotonic zaxis!". The warning isn't clear which z-axis it is referring to (from the input s-lvl grid with 30 levels, or the output z-lvl grid with 102 levels)?
But both z-axes are perfectly monotonically increasing, but if and only if the NaNs are correctly ignored.
The third thing is the abort error "Level 0 not found!" - I cannot make sense of that. Does cdo not recognise the missing values in the input grid? Where is it trying to find level 0?
I tried replacing all NaNs in the input grid MODEL_gdept.nc with 0, but that didn't work either.
Can anyone help??
Thanks,
Fred
Replies (6)
RE: format of input files for intlevel3d - Added by Ralf Mueller almost 5 years ago
hi Fred!
in order to have a reproducer in my hands I'd need all the input files. Can you please upload them?
I will have a look into this with the current 1.9.6 release and the development version on CDO - which version did you use for your tests?
cheers
ralf
RE: format of input files for intlevel3d - Added by Fred W almost 5 years ago
Yes, I was using the current 1.9.6 release (where I corrected one bit of source by hand to make it compile as instructed here in the forum).
$ cdo -V
Climate Data Operators version 1.9.6 (http://mpimet.mpg.de/cdo)
System: x86_64-pc-linux-gnu
CXX Compiler: g++ -std=gnu++11 -g -O2 -fopenmp
CXX version : g++ (GCC) 4.8.5 20150623 (Red Hat 4.8.5-11)
C Compiler: gcc -std=gnu99 -fPIC -fopenmp
C version : gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-11)
F77 Compiler: gfortran -g -O2
F77 version : GNU Fortran (GCC) 4.8.5 20150623 (Red Hat 4.8.5-11)
Features: 62GB 40threads C++11 Fortran DATA PTHREADS OpenMP3 HDF5 NC4/HDF5 OPeNDAP SSE2
Libraries: HDF5/1.10.5
Filetypes: srv ext ieg grb1 grb2 nc1 nc2 nc4 nc4c nc5
CDI library version : 1.9.6
cgribex library version : 1.9.2
ecCodes library version : 1.28.0
NetCDF library version : 4.6.3 of Apr 23 2019 11:44:45 $
hdf5 library version : 1.10.5
exse library version : 1.4.1
FILE library version : 1.8.3
| MODEL_gdept.nc (5.13 MB) MODEL_gdept.nc | original 3D depths of the model in variable gdept | ||
| MODEL_std_vertgrid.nc (16.3 MB) MODEL_std_vertgrid.nc | |||
| MODEL_gdept_with_NaNs.nc (5.13 MB) MODEL_gdept_with_NaNs.nc | like MODEL_gdept.nc, but zero depths replaced with NaN |
RE: format of input files for intlevel3d - Added by Fred W almost 5 years ago
I also attach an example of an original model output file.
The data is given on depth levels specified by gdept in MODEL_gdept.nc
what I am trying achieve is to regrid the same data onto the standard depth levels specified in MODEL_std_vertgrid.nc
So I was expecting model_output_regrid.nc to have the same horizontal grid, but standard depth levels.
My plan is to use cdo genbil and cdo remap to change the 2D horizontal grid into a regular grid with 1D coordinate vectors in each horizontal dimension.
(Note, that for this model the 2D coordinates are already regular - the 2D coordinates were simply generated by meshgrid() from 1D vectors, but for other data I have that isn't the case and the horz grid is irregular).
Thanks
Fred
| model_output_original.nc (24.6 MB) model_output_original.nc | example of an original model output file |
RE: format of input files for intlevel3d - Added by Ralf Mueller almost 5 years ago
hi Fred!
- I think what you want to do goes beyond the scope of
intlevel3d. This operator was designed to interpolate one 3D vertical coordinate to another one. Both have to have the same unit sinceintlevel3ddoes not do any unit conversion (possibly needed when switching between s- and z-coordinate - I am not an expert of vertical coordinates in ocean models). - the vertical coordinated are quite different on your uploads: The original model output has only a 1D-vertical axis (same list if heights for all horizontal locations).
MODEL_gdepth.nctruly holds a 3d-depth variable, but the gdepth variable inMODEL_std_vertgrid.ncis a time-dependent one. Such fields are not recognized as coordinates in CDO.
Sorry, but I don't think this is solvable with CDO atm.
btw: for handling missing values you need to use the value given in the _FillValue attribute of a variable. NaN is not a placeholder for this
RE: format of input files for intlevel3d - Added by Fred W almost 5 years ago
Hello Ralph,
Thanks for the explanation.
(1) The content in MODEL_std_vertgrid.nc is fully under my control, so I could change it into any format required.
(2) The 1D vertical axis in MODEL_gdept.nc is fake - the metres in there don't refer to anything. I think it's the wrong axis for this sort of file - it would be better if it contained an axis variable of just model level numbers. I can change that file if that would help.
(3) The 3D depth variable gdept in MODEL_gdept.nc describes the depth of each grid point (to be precise it's a "g"lobal matrix of the "dep"th of each "T"-point - hence g-dep-t). I want to interpolate from those s-levels to z-levels.
(4) when I said NaN I was referring to the value printed by ncdump. In the Netcdf file itself the masking is done using the _FillValue attribute as usual.
I'll probably have to accept that this is currently impossible with cdo.
I have a half-working solution using scrip, but was hoping this could be done simpler with cdo (as already works for horizontal remapping).
The fact that cdo doesn't precompute vertical weights is a give-away.
Can we have a cdo operator cdo genvertgrid ... and cdo vertregrid? Or (like my half-working scrip solution) one that does horizontal and vertical regridding in a single step once all the weights files have been computed.
Thanks,
Fred
RE: format of input files for intlevel3d - Added by Ralf Mueller almost 5 years ago
Hi Fred!
I had a deeper look into your uploaded files. Here is my summary:
MODEL_gdept.ncthis file is the model-internal has no missing values but zeros% cdo infov MODEL_gdept.nc -1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter name 1 : 0000-00-00 00:00:00 3.14595 41968 0 : 0.0000 0.24182 3.3863 : gdept 2 : 0000-00-00 00:00:00 10.4168 41968 0 : 0.0000 0.72847 10.646 : gdept 3 : 0000-00-00 00:00:00 19.2562 41968 0 : 0.0000 1.2208 18.830 : gdept 4 : 0000-00-00 00:00:00 30.0452 41968 0 : 0.0000 1.7218 28.403 : gdept 5 : 0000-00-00 00:00:00 43.2552 41968 0 : 0.0000 2.2352 40.044 : gdept 6 : 0000-00-00 00:00:00 59.467 41968 0 : 0.0000 2.7677 54.739 : gdept 7 : 0000-00-00 00:00:00 79.3944 41968 0 : 0.0000 3.3275 73.889 : gdept 8 : 0000-00-00 00:00:00 103.911 41968 0 : 0.0000 3.9268 99.431 : gdept 9 : 0000-00-00 00:00:00 134.08 41968 0 : 0.0000 4.5813 133.92 : gdept 10 : 0000-00-00 00:00:00 171.189 41968 0 : 0.0000 5.3088 180.47 : gdept 11 : 0000-00-00 00:00:00 216.786 41968 0 : 0.0000 6.1312 242.43 : gdept 12 : 0000-00-00 00:00:00 272.709 41968 0 : 0.0000 7.0632 322.54 : gdept 13 : 0000-00-00 00:00:00 341.123 41968 0 : 0.0000 8.1131 421.61 : gdept 14 : 0000-00-00 00:00:00 424.536 41968 0 : 0.0000 9.2619 537.09 : gdept 15 : 0000-00-00 00:00:00 525.805 41968 0 : 0.0000 10.470 662.45 : gdept 16 : 0000-00-00 00:00:00 648.113 41968 0 : 0.0000 11.682 788.43 : gdept 17 : 0000-00-00 00:00:00 794.909 41968 0 : 0.0000 12.834 905.76 : gdept 18 : 0000-00-00 00:00:00 969.81 41968 0 : 0.0000 13.891 1008.0 : gdept 19 : 0000-00-00 00:00:00 1176.45 41968 0 : 0.0000 14.837 1092.8 : gdept 20 : 0000-00-00 00:00:00 1418.28 41968 0 : 0.0000 15.682 1161.1 : gdept 21 : 0000-00-00 00:00:00 1698.36 41968 0 : 0.0000 16.435 1215.9 : gdept 22 : 0000-00-00 00:00:00 2019.09 41968 0 : 0.0000 17.126 1260.9 : gdept 23 : 0000-00-00 00:00:00 2382.06 41968 0 : 0.0000 17.779 1299.8 : gdept 24 : 0000-00-00 00:00:00 2787.84 41968 0 : 0.0000 18.414 1335.5 : gdept 25 : 0000-00-00 00:00:00 3235.98 41968 0 : 0.0000 19.040 1370.9 : gdept 26 : 0000-00-00 00:00:00 3725.03 41968 0 : 0.0000 19.676 1408.1 : gdept 27 : 0000-00-00 00:00:00 4252.67 41968 0 : 0.0000 20.343 1449.2 : gdept 28 : 0000-00-00 00:00:00 4815.91 41968 0 : 0.0000 21.054 1496.4 : gdept 29 : 0000-00-00 00:00:00 5411.34 41968 0 : 0.0000 21.818 1551.6 : gdept 30 : 0000-00-00 00:00:00 6035.34 41968 0 : 0.0000 0.0000 0.0000 : gdeptObviously the levels are non-monotonic because of the zero-layer at the bottom.MODEL_gdept_with_NaNs.ncUsing Nans does not lead to correct missing-values% cdo infov MODEL_gdept_with_NaNs.nc -1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter name 1 : 0000-00-00 00:00:00 3.14595 41968 0 : 0.17241 nan 3.3863 : gdept 2 : 0000-00-00 00:00:00 10.4168 41968 0 : 0.51724 nan 10.646 : gdept 3 : 0000-00-00 00:00:00 19.2562 41968 0 : 0.86207 nan 18.830 : gdept 4 : 0000-00-00 00:00:00 30.0452 41968 0 : 1.2069 nan 28.403 : gdept 5 : 0000-00-00 00:00:00 43.2552 41968 0 : 1.5517 nan 40.044 : gdept 6 : 0000-00-00 00:00:00 59.467 41968 0 : 1.8966 nan 54.739 : gdept 7 : 0000-00-00 00:00:00 79.3944 41968 0 : 2.2414 nan 73.889 : gdept 8 : 0000-00-00 00:00:00 103.911 41968 0 : 2.5862 nan 99.431 : gdept 9 : 0000-00-00 00:00:00 134.08 41968 0 : 2.9310 nan 133.92 : gdept 10 : 0000-00-00 00:00:00 171.189 41968 0 : 3.2759 nan 180.47 : gdept 11 : 0000-00-00 00:00:00 216.786 41968 0 : 3.6207 nan 242.43 : gdept 12 : 0000-00-00 00:00:00 272.709 41968 0 : 3.9655 nan 322.54 : gdept 13 : 0000-00-00 00:00:00 341.123 41968 0 : 4.3103 nan 421.61 : gdept 14 : 0000-00-00 00:00:00 424.536 41968 0 : 4.6552 nan 537.09 : gdept 15 : 0000-00-00 00:00:00 525.805 41968 0 : 5.0000 nan 662.45 : gdept 16 : 0000-00-00 00:00:00 648.113 41968 0 : 5.3448 nan 788.43 : gdept 17 : 0000-00-00 00:00:00 794.909 41968 0 : 5.6897 nan 905.76 : gdept 18 : 0000-00-00 00:00:00 969.81 41968 0 : 6.0345 nan 1008.0 : gdept 19 : 0000-00-00 00:00:00 1176.45 41968 0 : 6.3793 nan 1092.8 : gdept 20 : 0000-00-00 00:00:00 1418.28 41968 0 : 6.7241 nan 1161.1 : gdept 21 : 0000-00-00 00:00:00 1698.36 41968 0 : 7.0690 nan 1215.9 : gdept 22 : 0000-00-00 00:00:00 2019.09 41968 0 : 7.4138 nan 1260.9 : gdept 23 : 0000-00-00 00:00:00 2382.06 41968 0 : 7.7586 nan 1299.8 : gdept 24 : 0000-00-00 00:00:00 2787.84 41968 0 : 8.1034 nan 1335.5 : gdept 25 : 0000-00-00 00:00:00 3235.98 41968 0 : 8.4483 nan 1370.9 : gdept 26 : 0000-00-00 00:00:00 3725.03 41968 0 : 8.7931 nan 1408.1 : gdept 27 : 0000-00-00 00:00:00 4252.67 41968 0 : 9.1379 nan 1449.2 : gdept 28 : 0000-00-00 00:00:00 4815.91 41968 0 : 9.4828 nan 1496.4 : gdept 29 : 0000-00-00 00:00:00 5411.34 41968 0 : 9.8276 nan 1551.6 : gdept 30 : 0000-00-00 00:00:00 6035.34 41968 0 : 2000.0 nan 2000.0 : gdeptsince it is not the value given at the_FillValueattribute.
I removed the lowest level with cdo sellevidx,1/29 from the model outout (model_output_original.nc) and MODEL_gdept.nc and reran the vertical interpolation. The error with the non-monoton axis disappeared, but still there are points where all levels are 0. Hence there cannot be a reasonable vertical interpolation.
I designed the intlevel3d operator for atmosphere hight-fields that are used for the ICON-atmosphere model. These do not have any missing values, so I don't expect it to work with 3d ocean data. Regarding you suggestion for new operators (genvertgrid/vertregrid) the only availbe operators is setzaxis. But this is only for a simple 1D vertical coordinate. In general the vertical coordinates of models are very very different. Modellers create them just like that ;-) IMO there cannot be a single operator to generate them all - only some basic methods can be offered to change things.
hth
ralf