Version 5 (modified by nicolasmartin, 3 years ago) (diff)

Interpolation with CDO

Last edition: 04/12/17 23:07:16 by nicolasmartin

You can find here an example of bash script to interpolate forcing core2 files.

To interpolate files you need to do 3 steps:

  1. have good grid description file (if not you can create it, see instructions 1.)
  2. prepare files to be interpolated
  3. remapping on ORCA2 grid

NOTE: This script works with cdo version 1.4.0.1

1) Prepare grid description file for the ORCA curvilinear grid

You need to prepare a netCDF grid description file for the ORCA curvilinear grid.
This file must contain:

  • a longitude variable with exactly 2 dimensions and the attribute: units = "degrees_east"
  • a latitude variable with exactly 2 dimensions and the attribute: units = "degrees_north"
  • a dummy variable with the attribute: coordinates = "xxx yyy" where xxx and yyy are the name of the longitude and latitude variables mentioned just above

There is an example of a netCDF grid description file for the T grid of ORCA2:

netcdf grid_ORCA2_T {
dimensions:
        y = 149 ;
        x = 182 ;
variables:
        byte dummy(y, x) ;
                dummy:coordinates = "glamt gphit" ;
        double glamt(y, x) ;
                glamt:units = "degrees_east" ;
        double gphit(y, x) ;
                gphit:units = "degrees_north" ;
}

There is an example of nco commands to create this file from coordinates.nc (input file used by NEMO for its grid description)

# select glamt and gphit variables from the coordinates file
ncks -O -C -a -v glamt,gphit coordinates_ORCA_R2.nc grid_ORCA2_T.nc

make sure that coordinates variables contains only 2 dimension. Use ncwa -a to remove degenerated dimensions (with a size of 1) for example:

# remove degenerated dimention time (if existing)
ncwa -O -a time grid_ORCA2_T.nc  grid_ORCA2_T.nc
# remove degenerated dimention z (if existing)
ncwa -O -a z grid_ORCA2_T.nc  grid_ORCA2_T.nc
# add a dummy variable
ncap2 -O -s 'dummy[y,x]=1b' grid_ORCA2_T.nc grid_ORCA2_T.nc
# add needed attributes
ncatted -a coordinates,dummy,c,c,'glamt gphit' \
        -a units,glamt,c,c,'degrees_east'      \
        -a units,gphit,c,c,'degrees_north'  grid_ORCA2_T.nc

2) Prepare file to be interpolated

For a curvilinear grid:

  • you have to follow above spefication

For a regular lon/lat grid:

  • The coordinate variables must have the same name as the dimensions. You can chose any name. If you use the name x for the dimension then you have to use also x for the coordinate variable.
  • The coordinate variable must also have the attribute units = "degrees_east" or "degrees_north" to identify the lon/lat variables:
    netcdf etopo5 {
     dimensions:
             x = 4320 ;
             y = 2161 ;
     variables:
             float x(x) ;
                      x:units = "degree_east" ;
             float y(y) ;
                     y:units = "degree_north" ;
             float bath(y, x) ;
     }
    

3) Remapping on ORCA2

Generate weights with bilinear interpolation, and then remapping for a scalar field:

    cdo genbil,grid_ORCA2_T.nc $file bil_orca2_weights.nc
    cdo remap,grid_ORCA2_T.nc,bil_orca2_weights.nc $file_fill.nc $file_orca2.nc

where

  • $file_fill.nc is a file whose missing values is filled (see below interpolation_loop.sh)
  • $file_orca2.nc is the output
  • and bil_orca2_weights.nc is a weight file generated (see below)

Generate weights with bilinear interpolation, and then remapping for a vectorial field:

    cdo genbic,grid_ORCA2_T.nc $file bic_orca2_weights.nc
    cdo remap,grid_ORCA2_T.nc,bic_orca2_weights.nc $file_fill.nc $file_orca2.nc

where

  • $file_fill.nc is a file whose missing values is filled (see below interpolation_loop.sh)
  • $file_orca2.nc is the output
  • and bic_orca2_weights.nc is a weight file generated (see below)

Example of interpolation_loop.sh script

This script needs to have input directory in your execution's directory to copy *.nc files and it calls 4 others sub-scripts: create_orca2.sh, prepare_mask.sh, bilweights.sh and bicweights.sh

NOTE: Before running interpolation_loop.sh, you nedd to checks different things:

  1. control latitude of your mask file. If it goes from north to south and latitude in your file goes from south to north you have to invert it:
    cdo invertlat ${file_mask_land}.nc ${file_mask_land}.SN.nc
    
  1. control format of land variable, if it is short and your file has float format you will have problems with cdo created missing values, and filling missing values
    #convert in float (from short) land variable
    ncap -O -s "land=float(land)" ${file_mask_land}.SN.nc ${file_mask_land}.SN.float.nc
    
  1. control values of your land variable, if it is 1 on earth and 0 on ocean is not good if you are interpolating files for ocean forcing
    #convert land=0 on earth, 1 on ocean
    ncap -O -s 'land=1-land' mask.nc mask.nc
    
    #rename "land" in "mask"
    ncrename -v land,mask mask.nc
    
  1. control values of calendar's attribute in your file, if it is "NOLEAP" it is better to convert it in a recognized format by nccated (ex. 365_day)
    ncatted -a calendar,${variable},m,c,365_day ${file}.nc ${file}_365.nc
    

main programm: interpolation_loop.sh

#!/bin/bash
#script to interpolate files on ORCA2 grids

set -u
LC_ALL=C; export LC_ALL

type cdo
status_cdo=${?}
if [ ${status_cdo} != 0 ] ; then
        echo "ERROR: cdo not found in your PATH"
        exit 1
fi

file_orca2="grid_ORCA2_T.nc"
#build ORCA2_grid if necessary
if [ ! -f ${file_orca2} ] ; then
 /bin/sh create_orca2.sh
fi

#preparing mask file if necessary
if [ ! -f mask.nc ] ; then
 /bin/sh prepare_mask.sh
fi

dir=/PATH/FILES/TO/BE/INTERPOLATED

`cp ${dir}/*.nc input/.`
cd input

for file in $(echo `ls *.nc | sed -e "s/\.nc//"`) ; do
 
 # create missing values over land with e.g. 'ifthen':
 echo "create missing values ifthen mask"
 cdo ifthen ../mask.nc ${file}_365.nc ${file}_mask.nc
 status_ifthen=${?}
 if [ ${status_ifthen} != 0 ] ; then
         echo "ERROR ifthen mask"
         exit 1
 fi
 
 
 # fill missing value (here over land):
 echo "fill missing values"
 cdo fillmiss ${file}_mask.nc ${file}_fill.nc
 status_fillmiss=${?}
 if [ ${status_fillmiss} != 0 ] ; then
  echo "ERRORE cdo fillmiss"
         exit 1
 fi
 
 
 #remapping on ORCA2 grid
 echo "remapping "
 cd ..
# looking for first part of filename, to see if it is vectorial field, for interannual files
 name=$(echo $file | awk -F \_ {'print $1'} )
 if [ ${name} = "u" -o ${name} = "v" ] ; then
    echo "vectorial field"
    /bin/sh bicweights.sh input/${file}_fill.nc
    cdo remap,grid_ORCA2_T.nc,bic_orca2_weights.nc input/${file}_fill.nc input/${file}_orca2.nc
 else 
    echo "scalar field"
    /bin/sh bilweights.sh input/${file}_fill.nc
    cdo remap,grid_ORCA2_T.nc,bil_orca2_weights.nc input/${file}_fill.nc input/${file}_orca2.nc
 fi
 status_remap=${?}
 if [ ${status_remap} != 0 ] ; then
         echo "ERROR cdo remap"
         exit 1
 fi
 echo "status remap: ${status_remap}"
 cd input
done

first sub-programm: file create_orca2.sh

#!/bin/bash
#script to create ORCA2 grid

set -u
LC_ALL=C ; export LC_ALL

echo "create orca2 grid"
ncks -O -C -a -v glamt,gphit coordinates.nc grid_ORCA2_T.nc 

#remove not necessar axes:
ncwa -O -a time,z grid_ORCA2_T.nc grid_ORCA2_T.nc

#add dummy variable:
ncap -O -s 'dummy[y,x]=1b' grid_ORCA2_T.nc grid_ORCA2_T.nc

#add needed attributes
ncatted -a coordinates,dummy,c,c,'glamt gphit'  \
        -a units,glamt,c,c,'degrees_east'       \
        -a units,gphit,c,c,'degrees_north' grid_ORCA2_T.nc

Second sub-programm: prepare_mask.sh

#!/bin/bash
#preparing mask file

set -u
LC_ALL=C ; export LC_ALL

file_mask_land="land.sfc.gauss"
# invert latitude (from SOUTH to NORTH)
cdo invertlat ${file_mask_land}.nc ${file_mask_land}.SN.nc

#convert in float (from short) land variable
ncap -O -s "land=float(land)" ${file_mask_land}.SN.nc ${file_mask_land}.SN.float.nc 

# remove all attribute to land
ncatted -a ,land,d,, ${file_mask_land}.SN.float.nc ${file_mask_land}.SN.float.miss.nc

# extract land
#ncks -O -v land land.sfc.gauss.SN.float.nc mask.nc 
ncks -O -v land ${file_mask_land}.SN.float.miss.nc mask.nc 
ncwa -O -a time mask.nc mask.nc

#convert land=0 on earth, 1 on ocean
ncap -O -s 'land=1-land' mask.nc mask.nc

#rename "land" in "mask"
ncrename -v land,mask mask.nc
chmod 755 mask.nc

Third sub-programm: bicweights.sh

#!/bin/bash

#remapping on ORCA2 grid
echo "gen bicubic WEIGHTS"
export REMAP_EXTRAPOLATE=on
cdo genbic,grid_ORCA2_T.nc $1 bic_orca2_weights.nc

Forth sub-programm: bilweights.sh

#!/bin/bash

#remapping on ORCA2 grid
echo "gen bilinear WEIGHTS"
export REMAP_EXTRAPOLATE=on
cdo genbil,grid_ORCA2_T.nc $1 bil_orca2_weights.nc