Project

General

Profile

cdo remapcon adds mysterious unnecessary interpolation wh... » CDORemapcon_example.txt

Toby Marthews, 2022-07-04 12:27

 
## WARNING: these commands will overwrite any file called jkl, SixPointsOnly.nc, regridfile_tmp1.nc, regridfile_tmp2.nc, regridfile_lonbnds.nc and gridfile.txt in your current directory. ##


setup cdo

cat >jkl <<EOF
netcdf SixPointsOnly {
dimensions:
ncells = 6 ;
variables:
double lon(ncells) ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:units = "degrees_east" ;
double lat(ncells) ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:units = "degrees_north" ;
float elevation(ncells) ;
elevation:units = "m asl" ;
elevation:title = "elevation" ;
elevation:grid_type = "unstructured" ;
elevation:coordinates = "lon lat" ;
elevation:_FillValue = -1.e+20f ;

data:

lon = -179.625, -161.125, -152.875, -133.625, -133.625, -127.125 ;

lat = 71.125, 70.875, 68.625, 68.625, 67.375, 66.375 ;

elevation = 392, 116, 42, 30, 91, 221 ;
}
EOF
ncgen -o SixPointsOnly.nc jkl
rm jkl
ncview SixPointsOnly.nc &
#You should now see the image SixPointsOnly.png

cat >gridfile.txt <<EOF
gridtype = lonlat
gridsize = 1036800
xsize = 1440
ysize = 720
xname = lon
xlongname = longitude
xunits = degrees_east
yname = lat
ylongname = latitude
yunits = degrees_north
xfirst = -179.875
xinc = 0.25
yfirst = -89.875
yinc = 0.25
EOF
























#REGRIDDING UNSTRUCTURED -> STRUCTURED
ncap2 -Oh -s "latN=lat+0.125" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a long_name,latN,o,c,"Gridcell bound N" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,latN,o,c,"Gridcell bound N" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_min,latN,d,f, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_max,latN,d,f, SixPointsOnly.nc SixPointsOnly.nc

ncap2 -Oh -s "latS=lat-0.125" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a long_name,latS,o,c,"Gridcell bound S" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,latS,o,c,"Gridcell bound S" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_min,latS,d,f, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_max,latS,d,f, SixPointsOnly.nc SixPointsOnly.nc

ncap2 -Oh -s "lonE=lon+0.125" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a long_name,lonE,o,c,"Gridcell bound E" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,lonE,o,c,"Gridcell bound E" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_min,lonE,d,f, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_max,lonE,d,f, SixPointsOnly.nc SixPointsOnly.nc

ncap2 -Oh -s "lonW=lon-0.125" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a long_name,lonW,o,c,"Gridcell bound W" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,lonW,o,c,"Gridcell bound W" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_min,lonW,d,f, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a valid_max,lonW,d,f, SixPointsOnly.nc SixPointsOnly.nc


ncatted -Oh -a grid_type,elevation,o,c,"unstructured" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a coordinates,elevation,o,c,"lon lat" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a _FillValue,elevation,o,f,1e10 SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a missing_value,elevation,o,f,1e10 SixPointsOnly.nc SixPointsOnly.nc

ncatted -Oh -a long_name,lon,o,c,"longitude" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a long_name,lat,o,c,"latitude" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a standard_name,lon,o,c,"longitude" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a standard_name,lat,o,c,"latitude" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a units,lon,o,c,"degrees_east" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a units,lat,o,c,"degrees_north" SixPointsOnly.nc SixPointsOnly.nc
ncks -Oh -4 SixPointsOnly.nc SixPointsOnly.nc

#Longitude bounds
ncks -Oh -v lonE SixPointsOnly.nc regridfile_tmp1.nc
ncrename -Oh -v lonE,lon_bnds regridfile_tmp1.nc regridfile_tmp1.nc
ncks -Oh -v lonW SixPointsOnly.nc regridfile_tmp2.nc
ncrename -Oh -v lonW,lon_bnds regridfile_tmp2.nc regridfile_tmp2.nc

#Go round anti-clockwise starting at SE corner:
ncecat -Oh regridfile_tmp2.nc regridfile_tmp1.nc regridfile_tmp1.nc regridfile_tmp2.nc regridfile_lonbnds.nc

ncrename -Oh -d record,bnds regridfile_lonbnds.nc regridfile_lonbnds.nc
ncks -Oh --fix_rec_dmn bnds regridfile_lonbnds.nc regridfile_lonbnds.nc
rm regridfile_tmp1.nc regridfile_tmp2.nc
ncpdq -Oh -v lon_bnds -a land,bnds regridfile_lonbnds.nc regridfile_lonbnds.nc
ncks -Ah -v lon_bnds regridfile_lonbnds.nc SixPointsOnly.nc
rm regridfile_lonbnds.nc
ncatted -Oh -a long_name,lon_bnds,o,c,"Gridcell west-east bounds" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,lon_bnds,d,c, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a units,lon_bnds,o,c,"degrees_east" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a bounds,lon,o,c,"lon_bnds" SixPointsOnly.nc SixPointsOnly.nc

#Latitude bounds
ncks -Oh -v latN SixPointsOnly.nc regridfile_tmp1.nc
ncrename -Oh -v latN,lat_bnds regridfile_tmp1.nc regridfile_tmp1.nc
ncks -Oh -v latS SixPointsOnly.nc regridfile_tmp2.nc
ncrename -Oh -v latS,lat_bnds regridfile_tmp2.nc regridfile_tmp2.nc

#Go round anti-clockwise starting at SE corner:
ncecat -Oh regridfile_tmp2.nc regridfile_tmp2.nc regridfile_tmp1.nc regridfile_tmp1.nc regridfile_latbnds.nc

ncrename -Oh -d record,bnds regridfile_latbnds.nc regridfile_latbnds.nc
ncks -Oh --fix_rec_dmn bnds regridfile_latbnds.nc regridfile_latbnds.nc
rm regridfile_tmp1.nc regridfile_tmp2.nc
ncpdq -Oh -v lat_bnds -a land,bnds regridfile_latbnds.nc regridfile_latbnds.nc
ncks -Ah -v lat_bnds regridfile_latbnds.nc SixPointsOnly.nc
rm regridfile_latbnds.nc
ncatted -Oh -a long_name,lat_bnds,o,c,"Gridcell north-south bounds" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a title,lat_bnds,d,c, SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a units,lat_bnds,o,c,"degrees_north" SixPointsOnly.nc SixPointsOnly.nc
ncatted -Oh -a bounds,lat,o,c,"lat_bnds" SixPointsOnly.nc SixPointsOnly.nc

#Strangely, if I leave the lonW,lonE,latN,latS variables in the file, cdo remapcon aborts with the error "cdo remapcon (Abort): Unsupported grid type: generic"
ncks -Oh -x -v lonW,lonE,latN,latS SixPointsOnly.nc SixPointsOnly.nc


cdo -P 4 remapcon,gridfile.txt SixPointsOnly.nc out.nc
#Not sure why, but the cdo remapcon command here produces a lot of warnings (even though it does produce an output file):
#wllf003 ancils $ cdo -P 4 remapcon,gridfile.txt SixPointsOnly.nc out.nc
#cdo remapcon: SCRIP first order conservative remapping from unstructured (6) to lonlat (1440x720) grid
#cdo remapcon: 75%
cdo remapcon: Map weight > 1! grid1idx=0 grid2idx=925953 nlink=98 wts=1.45193
#cdo remapcon: Map weight > 1! grid1idx=0 grid2idx=925952 nlink=100 wts=557.318
#[... 197 'Map weight' errors in total ...]
#cdo remapcon: Map weight > 1! grid1idx=3 grid2idx=927367 nlink=1560 wts=2262.35
#cdo remapcon: Processed 6 values from 1 variable over 1 timestep ( 10.61s )
#wllf003 ancils $




rm SixPointsOnly.nc
rm gridfile.txt

ncview -extra -noautoflip -no_auto_overlay -minmax all out.nc &
(1-1/5)