cdo remap rotated grid help
Added by Melchior van Wessem about 7 years ago
Dear all,
I am trying to reproject data from a rotated grid to a simple rectangular grid. In the past I succeeded in getting this to work, but with CDO 1.9.0 it doesn't seem to recognize my grid definitions anymore.
My griddef= https://www.dropbox.com/s/4lq9q7z6c12btlj/grid_ANT27.txt?dl=0
- gridID 0
#
gridtype = lonlat
gridsize = 62880
xname = rlon
xlongname = longitude in rotated grid
xunits = degrees
yname = rlat
ylongname = latitude in rotated grid
yunits = degrees
xsize = 262
ysize = 240
xfirst = -32.75
xinc = 0.25
yfirst = -30
yinc = 0.25
xnpole = -170.
ynpole = 180.
When I run my script, that includes:
cdo remapcon,$outgrid -setgrid,$ingrid $infile $outfile
it does not recognize the xnpole and ynpoles anymore, and fails in remapping: "Invalid parameter : >xnpole<". Has these projection parameters been changed in the cdo updates? Please help me in getting these interpolations done! Thanks a lot.
Melchior
Replies (7)
RE: cdo remap rotated grid help - Added by Karin Meier-Fleischer about 7 years ago
Hi Melchior,
using a rotated pole data file which has the following definition for the rotated_pole e.g.
char rotated_pole ; rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ; rotated_pole:grid_north_pole_latitude = 39.25 ; rotated_pole:grid_north_pole_longitude = -162. ;
the grid information among others depending on the rotated pole can be retrieved using
cdo griddes infile
The rotated pole grid information variables are not xnpole and ynpole, but you get
grid_mapping = rotated_pole grid_mapping_name = rotated_latitude_longitude grid_north_pole_latitude = 39.25 grid_north_pole_longitude = -162.
Hope this will help.
-Karin
RE: cdo remap rotated grid help - Added by Melchior van Wessem about 7 years ago
Hi karin,
Thanks for your reply. Unfortunately I can't seem to end up with a good interpolation that way.
I get something like this:
#
- gridID 6 #
gridtype = curvilinear
gridsize = 62880
xsize = 262
ysize = 240
xname = lon
xdimname = rlon
xlongname = "longitude"
xunits = "degrees_east"
yname = lat
ydimname = rlat
ylongname = "latitude"
yunits = "degrees_north"
gridtype = projection
gridsize = 62880
xsize = 262
ysize = 240
xname = rlon
xlongname = "longitude in rotated pole grid"
xunits = "degrees"
yname = rlat
ylongname = "latitude in rotated pole grid"
yunits = "degrees"
xfirst = -32.75
xinc = 0.25
yfirst = -30
yinc = 0.25
grid_mapping = rotated_pole
grid_mapping_name = rotated_latitude_longitude
grid_north_pole_latitude = -180.f
grid_north_pole_longitude = -170.f
proj4_params = "-m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lat_p=-180.0 +lon_0=10.0"
proj_parameters = "-m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lat_p=-180.0 +lon_0=10.0"
projection_name = "rotated_latitude_longitude"
> long_name = "projection details" I attached the file I wanted to reproject (on anything), and it doesn't seem to work. We are familiar with the fact that our rotated grid makes no sense and hardly any people can work with it, haha. It worked easily with the previous grid def files (the ones with xnpole and ynpole), but since CDO 1.8 it does not work anymore. Only way now seems to be to just point-wise interpolate. I hope you can help further! Melchior
RE: cdo remap rotated grid help - Added by Karin Meier-Fleischer about 7 years ago
Hi Melchior,
the rotated data can be remaped to a 0.25 degrees latlon grid using the following grid file new_grid.txt
gridtype = lonlat xsize = 1440 ysize = 720 xfirst = 0.0 xinc = 0.25 yfirst = -90.0 yinc = 0.25
The data is on the southern hemisphere, therefore you have to set the yfirst value to -90.
cdo -f nc -remapbil,new_grid.txt infile outfile
The plot shows the remaped data on a polar map:
plot_remapped.png (71.1 KB) plot_remapped.png | plot of remaped data |
RE: cdo remap rotated grid help - Added by Uwe Schulzweida about 7 years ago
Thanks for the example file! The problematic attribute is proj4_params. This is a special attribute for CDO. If this attribute is found the proj4 library is used together with all parameters found in the attribute proj4_params to convert the projected coordinates to geographical coordinates.
It seems that there is something wrong or missing in the proj4_params.
If the proj4_params attribute is missing then CDO is using it's own routines to calculate the geographic coordinates for rotated lon/lat grids. That means if you remove this attribute the results will be correct.
Cheers,
Uwe
RE: cdo remap rotated grid help - Added by Melchior van Wessem about 7 years ago
Karin Meier-Fleischer wrote:
Hi Melchior,
the rotated data can be remaped to a 0.25 degrees latlon grid using the following grid file new_grid.txt
[...]
The data is on the southern hemisphere, therefore you have to set the yfirst value to -90.
[...]
The plot shows the remaped data on a polar map:
So I must have a problem with my libraries surely, as Uwe mentions below. Because regridding it like this does not work at all. How can I fix that?
I have tried reinstalling cdo from my MacPorts, but to no avail.
RE: cdo remap rotated grid help - Added by Melchior van Wessem about 7 years ago
I get it to work with an older version of cdo (1.6.0) now, making it active again. Let's leave it at that....
RE: cdo remap rotated grid help - Added by Leo van Kampenhout about 7 years ago
Hello,
I bumped into the same issue as Melchior described, about 'xnpole' and 'ynpole' not being supported parameters (not surprisingly as my data is from the same RCM).
The problem lies partly with the data as the netCDF headers in these files are not complete. Therefore, when doing remapping, we use external grid description files (as .txt) as a workaround.
This worked fine until CDO 1.8 , when the the formatting of 'griddes' changed. Now, the poles are no longer described by xnpole and ynpole but by the parameter names that Karin mentioned.
Long story short, what I did was insert the "correct" grid definition into the NetCDF using CDO 1.7.2, and create a new grid description file (.txt) from that using CDO 1.8.2.
The differences are as follows:
Old¶
gridtype = lonlat gridsize = 95472 xname = rlon xlongname = longitude in rotated grid xunits = degrees yname = rlat ylongname = latitude in rotated grid yunits = degrees xsize = 306 ysize = 312 xfirst = 17.35 xinc = -0.1 yfirst = 15.8 yinc = -0.1 xnpole = -37.5 ynpole = -18.0
New¶
# # gridID 1 # gridtype = projection gridsize = 95472 xsize = 306 ysize = 312 xname = rlon xlongname = "longitude in rotated grid" xunits = "degrees" yname = rlat ylongname = "latitude in rotated grid" yunits = "degrees" xfirst = 17.35 xinc = 0.1 yfirst = 15.8 yinc = -0.0999999999999996 grid_mapping = rotated_pole grid_mapping_name = rotated_latitude_longitude grid_north_pole_latitude = -18. grid_north_pole_longitude = -37.5
As you can see, the x-increment changed from -0.1 to 0.1, which gives a wrong longitude range, as revealed by operator SINFO:
rlon : 17.35 to 47.85 by 0.09999999 degrees rlat : 15.8 to -15.3 by -0.1000004 degrees
whereas this should be
rlon : 17.35 to -13.15 by 0.1 degrees rlat : 15.8 to -15.3 by -0.1 degrees
My current understanding is that this is a bug. (Also, the y-increment is rounded differently, but probably this is negligible.)
Steps to repeat bug¶
Use CDO 1.7.2 to remap some random file to the first grid definition file (called grid_ZGRN11.txt below).
Use CDO 1.8.2 to generate the second grid description.
~/bin/cdo-1.7.2/src/cdo remapbil,grid_ZGRN11.txt tmp.nc old.nc cdo griddes old.nc > grid_ZGRN11_new.txt
Perhaps related, perhaps not, operator CONST suffers from a similar issues related to 'xinc'...
cdo const,1,grid_ZGRN11_new_fixed.txt tmp2.nc cdo griddes tmp2.nc .... xinc = 1.080328
Cheers, Leo