How to wrap around longitude?
Added by Laure Resplandy over 7 years ago
Hi,
I am working with nc ocean files that starts and ends in the middle of the Indian Ocean. I am trying to merge the 2 parts of the Indian Ocean to get a continuous basin.
I tried cdo sellonlatbox to change the boundaries of the file, which I think would be the easiest -> not working with these files.
I also tried some of the ideas on the forum with mergegrid or shiftx cyclic -> not working either.
Finally I tried to regrid with remapbil -> not working either.
I always get warnings about the grid and lon/lat definition.
The only thing I managed to do is to extract 2 subdomains of the Indian Ocean with
cdo selindexbox,3001,3600,251,1565 sample_file.nc file1.nc
cdo selindexbox,1,500,251,1565 sample_file.nc file2.nc
... but can't merge them into one file.
I can't even attach 1 time-step of the surface 2D field as it is 111 Mb (gridsize 9720000) but I attach the cdo sinfo and ncdump -h below.
Thank you very much for any piece of advise on this,
Best regards,
Laure
------------- cdo sinfo ------------------------------------------
Warning (cdf_set_var) : Inconsistent variable definition for geolat_c!
Warning (cdf_set_var) : Inconsistent variable definition for geolat_t!
Warning (cdf_set_var) : Inconsistent variable definition for geolon_c!
Warning (cdf_set_var) : Inconsistent variable definition for geolon_t!
Warning (find_time_vars) : Found more than one time variable, skipped variable average_T1!
Warning (find_time_vars) : Found more than one time variable, skipped variable average_T2!
File format : NetCDF2
-1 : Institut Source Steptype Levels Num Points Num Dtype : Parameter ID
1 : unknown unknown instant 50 1 9720000 1 F32 : -1
2 : unknown unknown instant 1 2 1 2 F64 : -2
Grid coordinates :
1 : curvilinear : points=9720000 (3600x2700)
geolon_t : -279.9712 to 1e+20 degrees_east
geolat_t : -81.10863 to 1e+20 degrees_east
mapping : undefined
xt_ocean : -279.95 to 79.95 by 0.1 degrees_E
yt_ocean : -81.10863 to 89.9789 degrees_N
2 : generic : points=1
Vertical coordinates :
1 : generic : levels=50
st_ocean : 5.03355 to 5395.023 meters
2 : surface : levels=1
Time coordinate : 1 step
RefTime = 0001-01-01 00:00:00 Units = days Calendar = 365_day Bounds = true
YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss
0195-07-02 12:00:00
--------- ncdump h ------------
dimensions:
time = UNLIMITED ; // (1 currently)
yu_ocean = 2700 ;
xu_ocean = 3600 ;
yt_ocean = 2700 ;
xt_ocean = 3600 ;
nv = 2 ;
st_ocean = 50 ;
st_edges_ocean = 51 ;
variables:
float geolat_c(yu_ocean, xu_ocean) ;
geolat_c:long_name = "uv latitude" ;
geolat_c:units = "degrees_N" ;
geolat_c:valid_range = -91.f, 91.f ;
geolat_c:missing_value = 1.e+20f ;
geolat_c:_FillValue = 1.e+20f ;
geolat_c:cell_methods = "time: point" ;
geolat_c:coordinates = "geolon_c geolat_c" ;
float geolat_t(yt_ocean, xt_ocean) ;
geolat_t:long_name = "tracer latitude" ;
geolat_t:units = "degrees_N" ;
geolat_t:valid_range = -91.f, 91.f ;
geolat_t:missing_value = 1.e+20f ;
geolat_t:cell_methods = "time: point" ;
geolat_t:coordinates = "geolon_t geolat_t" ;
float geolon_c(yu_ocean, xu_ocean) ;
geolon_c:long_name = "uv longitude" ;
geolon_c:units = "degrees_E" ;
geolon_c:valid_range = -281.f, 361.f ;
geolon_c:missing_value = 1.e+20f ;
geolon_c:_FillValue = 1.e+20f ;
geolon_c:cell_methods = "time: point" ;
geolon_c:coordinates = "geolon_c geolat_c" ;
float geolon_t(yt_ocean, xt_ocean) ;
geolon_t:long_name = "tracer longitude" ;
geolon_t:units = "degrees_E" ;
geolon_t:valid_range = -281.f, 361.f ;
geolon_t:missing_value = 1.e+20f ;
geolon_t:_FillValue = 1.e+20f ;
geolon_t:cell_methods = "time: point" ;
geolon_t:coordinates = "geolon_t geolat_t" ;
double nv(nv) ;
nv:long_name = "vertex number" ;
nv:units = "none" ;
nv:cartesian_axis = "N" ;
float o2(time, st_ocean, yt_ocean, xt_ocean) ;
o2:long_name = "Oxygen" ;
o2:units = "mol/kg" ;
o2:missing_value = -1.e+20f ;
o2:_FillValue = -1.e+20f ;
o2:coordinates = "geolon_t geolat_t" ;
o2:cell_methods = "time: mean" ;
o2:time_avg_info = "average_T1,average_T2,average_DT" ;
double st_edges_ocean(st_edges_ocean) ;
st_edges_ocean:long_name = "tcell zstar depth edges" ;
st_edges_ocean:units = "meters" ;
st_edges_ocean:cartesian_axis = "Z" ;
st_edges_ocean:positive = "down" ;
double st_ocean(st_ocean) ;
st_ocean:long_name = "tcell zstar depth" ;
st_ocean:units = "meters" ;
st_ocean:cartesian_axis = "Z" ;
st_ocean:positive = "down" ;
st_ocean:edges = "st_edges_ocean" ;
double time(time) ;
time:long_name = "time" ;
time:units = "days since 0001-01-01 00:00:00" ;
time:cartesian_axis = "T" ;
time:calendar_type = "NOLEAP" ;
time:calendar = "NOLEAP" ;
time:bounds = "time_bounds" ;
double time_bounds(time, nv) ;
time_bounds:long_name = "time axis boundaries" ;
time_bounds:units = "days" ;
time_bounds:missing_value = 1.e+20 ;
time_bounds:_FillValue = 1.e+20 ;
double xt_ocean(xt_ocean) ;
xt_ocean:long_name = "tcell longitude" ;
xt_ocean:units = "degrees_E" ;
xt_ocean:cartesian_axis = "X" ;
double xu_ocean(xu_ocean) ;
xu_ocean:long_name = "ucell longitude" ;
xu_ocean:units = "degrees_E" ;
xu_ocean:cartesian_axis = "X" ;
double yt_ocean(yt_ocean) ;
yt_ocean:long_name = "tcell latitude" ;
yt_ocean:units = "degrees_N" ;
yt_ocean:cartesian_axis = "Y" ;
double yu_ocean(yu_ocean) ;
yu_ocean:long_name = "ucell latitude" ;
yu_ocean:units = "degrees_N" ;
yu_ocean:cartesian_axis = "Y" ;
double average_T1(time) ;
average_T1:long_name = "Start time for average period" ;
average_T1:units = "days since 0001-01-01 00:00:00" ;
double average_T2(time) ;
average_T2:long_name = "End time for average period" ;
average_T2:units = "days since 0001-01-01 00:00:00" ;
double average_DT(time) ;
average_DT:long_name = "Length of average period" ;
average_DT:units = "days" ;
Replies (8)
RE: How to wrap around longitude? - Added by Ralf Mueller over 7 years ago
Hi Laure!
I got some questions:
- how did you try to merge the files?
- which CDO version do you use?
- can you upload some sample data to our ftp server?
- url: ftp://ftp.zmaw.de/incoming/
- howto: user anonymous, passwd email
- create a subdir for you and upload
cheers
ralf
btw. gzip work with nc files pretty good - you might get it under 70mb out of the box and skip the ftp server
RE: How to wrap around longitude? - Added by Laure Resplandy over 7 years ago
Hi,
Thanks for the quick reply,
how did you try to merge the files?
I tried cdo mergegrid after creating an enlarged file with the following:
#Make a file called "grid_indian.txt" and put in it the right info:
gridtype = curvilinear
gridsize = 1446500
xname = lon
xlongname = longitude
xunits = degrees east
yname = lat
ylongname = latitude
yunits = degrees north
xsize = 1100
ysize = 1315
xfirst = 20.05
xinc = 0.1
yfirst = 70.54879
yinc = 0.1
-------
- Extract the 2 subdomains
cdo selindexbox,3001,3600,251,1565 sample_file.nc file1.nc
cdo selindexbox,1,500,251,1565 sample_file.nc file2.nc
#First make an enlarged file with the info of file3. Its contents are rubbish but that doesnt matter:
cdo enlarge,grid_indian.txt file1.nc file4.nc
#Then merge the enlarged domain with the content of the first and consequently with the second file:
cdo mergegrid file4.nc file1.nc file5.nc
cdo mergegrid file5.nc file2.nc file6.nc
... However my problems is beyond just mergegrid, I can’t use none of the following: remapbil, shiftx or sellonlatbox
which CDO version do you use?
I installed cdo with macport
cdo @1.8.0_0
can you upload some sample data to our ftp server?
I followed your advice and gzip it. See attached sample_file.tar.gz, with only the 2D surface field.
Thanks for your help!
Laure
sample_file.tar.gz (26 MB) sample_file.tar.gz |
RE: How to wrap around longitude? - Added by Laure Resplandy over 7 years ago
Hi,
I was just wondering if you had any ideas based on the file that I sent?
Thanks again for your help,
Laure
RE: How to wrap around longitude? - Added by Ralf Mueller over 7 years ago
Sorry,Laure. I totally missed your message. currently I'm sitting in train, but i will have a look at it tomorrow.
Best wishes,
ralf
RE: How to wrap around longitude? - Added by Ralf Mueller over 7 years ago
Hi Laure!
here is my first analysis of the file:
% cdo sinfov sample_file.nc [Thu 2017-06-29|13:15:42] File format : NetCDF2 -1 : Institut Source Steptype Levels Num Points Num Dtype : Parameter name 1 : unknown unknown instant 1 1 9720000 1 F32 : o2 2 : unknown unknown instant 1 2 1 2 F64 : average_DT Grid coordinates : 1 : curvilinear : points=9720000 (3600x2700) geolon_t : -279.9712 to 1e+20 degrees_east geolat_t : -81.10863 to 1e+20 degrees_east mapping : undefined xt_ocean : -279.95 to 79.95 by 0.1 degrees_E yt_ocean : -81.10863 to 89.9789 degrees_N 2 : generic : points=1 Vertical coordinates : 1 : generic : levels=1 st_ocean : 5.03355 meters 2 : surface : levels=1 Time coordinate : 1 step RefTime = 0001-01-01 00:00:00 Units = days Calendar = 365_day Bounds = true YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss 0195-07-02 12:00:00
IMO there is something wrong with the coordinates. the original coordinates xt_ocean,yt_ocean
seem to be ok, but the maped ones (geolon_t,gelolat_t
) seem to hold wrong values (1e+20).
I checked geolon_t in detail:
- 2913273 values are 1.0000000200408773e20
- 6806727 are not, i.e. varying from -279.9711608886719 to 79.96676635742188
for CDO this looks like a mask, but in the coordinates! This is pretty much unusual, the _FillValue
attribute for the data variable should be used for this, but not the coordinates. Or did you intent to do any kind of projection with this second set of coordinates (geolon,geolat)?
With these values in the coordinates any kind on interpolation fails. I tried a nearest neighbour interpolation, but stoped it after a couple of minutes.
If you want, I can repair the file - just let me know
cheers
ralf
RE: How to wrap around longitude? - Added by Uwe Schulzweida over 7 years ago
The CDO operator mergegrid is only implemented for regular lon/lat grids and will not work on your curvilinear grid.
You can use selindexbox to select both parts in one step:
cdo selindexbox,3001,500,251,1565 sample_file.nc result.nc
RE: How to wrap around longitude? - Added by Laure Resplandy over 7 years ago
Thank you for the feedback.
This is the way the nc files are written by this climate model. I have used many models before and this is the first one I have this issue with.
In the new runs, I'll fix this but I already have hundreds of files with the same issue. I would be very very interested in being able to fix the problem in a systematic way and be able to use cdo to analyse them.
Would you know what the best options would be?
Thanks again!
-Laure
RE: How to wrap around longitude? - Added by Ralf Mueller over 7 years ago
We need to make sure, which coordinate is the correct one for o2? geolon_t,geolat_t or xt,yt. Could you find that out?
best
ralf