PARAMOUR_TOOLS/interp

CDO horizontal interpolation.

1. Read the CDO documentation.

2. Prepare CDO grid description files (see get_cdo_griddes.sh), which are kind of "blueprints" for CDO - it describes your grid.

3. Your input data should be properly masked. This means that the masked cells should be _FillValue, and appear in white when you look at them with ncview. This way, CDO will know that the fields are masked there, and that the values should not be used. If there's one "dummy" value that is not the specified _FillValue, then CDO will use these values as if they were good. Which you don't want to do. To properly mask a field, try to do something like that:

   ncatted -a _FillValue,<some_var>,o,f,1.e20 file.nc # in <some_var> is a SIMPLE precision FLOAT
   ncatted -a _FillValue,<some_var>,o,d,1.e20 file.nc # in <some_var> is a DOUBLE precision float
   

   And then:

   ncap2 -O -s 'where(mask==0) <some_var>=<some_var>@_FillValue;' file.nc output.nc
   

   This assumes that file.nc has a field called mask, and that cells where mask is zero ought to be masked out. <some_var> can have more dimensions than mask: for example, if mask(y,x) and <some_var>(t,y,x), things will be broadcasted along the t dimension (all t indices will be masked). Think about what you do: if you're masking a 3D (z,y,x) ocean field, use a 3D mask; the mask is not the same at the surface and at 5000m depth. Other things can be used in the where statement, e.g. where(abs(<some_var>) > 1.e10) will mask anywhere where <some_var> exceeds 1.e10.

4. You can "drown" the _FillValue (i.e., set masked cells to near values) by:

   • Using my nn_interp.py function: it's quite fast, and will fill a masked array with a the nearest neighbor in terms of discrete coordinates (i.e., array indices, not physical distance).
   • Using the undocument setmisstonn operator in CDO: cdo setmisstonn -setgrid,<cdo_griddes> masked.nc drowned. In terms of code, it's straightforward and reliable, but this actually computes loads of distances (orthodromy) to get the physically nearest neighbor. So it's slow. But reliable. Just do it overnight!
   • If your grid is regular (same lat/lon at each x/y), then better use fillmiss which fill masked values with linear reinterpolation. Warning: these drowning operations only work for array slices which have at least one unmasked cell in the (x,y) plane. If, on an array slice, all cells are masked, I think my nn_interp.py will just loop forever, and the CDO operators won't do anything. You probably need to do something else (like, replicate one vertical level to the other).

5. If you want to perform the same interpolation a lot of times, it may be a good idea to generate a weight file, and then use it. The weight file is kind of the DNA of the interpolation. If CDO has it, it just uses the weight, and does a linear combination, which is much faster than interpolation without having the weights in advance. To do so, for bilinear interpolation:

   cdo genbil,<dest_cdo_griddes> -setgrid,<src_cdo_griddes> input.nc weights.nc
   cdo remap,<dest_cdo_griddes>,weights.nc input.nc output.nc
   

   Warning: weight files depend on the mask of the input data in input.nc. If the input mask changes, e.g. if the data has a depth dimension along which the mask changes, this won't be useful.

6. Most times, conservative interpolation is "cleaner" but it's much harder to implement. In particular, you will need grid corners. Look at the CDO documentation and at the comments at the end of get_cdo_griddes.sh.

7. Conservative interpolation from a masked input field is doable by:

   • Putting zero's instead of _FillValue and removing the _FillValue attribute (this is an exception to point 3. above):

     ncatted -a _FillValue,<some_var>,o,f,0. input.nc # store zero's on masked cells, but the netCDF will still consider them masked
     ncatted -a _FillValue,<some_var>,d,, input.nc # delete _FillValue attribute: the previous _FillValue will be there, hence, zero's
     

   • Preparing a float or double precision input mask field, which is 0. on masked cells and 1. on unmasked cells.

   • Perform conservative interpolation of the input data and of the mask.

   • The result, on the destination grid and taking into account the mask, is data/mask. If mask == 0, this will look bad, but it's normal: you're looking at a destination masked cell (i.e. which has no overlap with any unmasked source cell).

CDO vertical interpolation

Just the same as horizontal interpolation, CDO uses description files for describing a vertical axis. They can be generated from netCDF files like this:

cdo zaxisdes input.nc > description.zaxis

I never quite cornered what input.nc must exactly look like, but it should contain a "depth" field whose standard name is something like "depth", with informations on how it's stored. I know that raw 3D NEMO output fields work OK with this type of depth description:

	float deptht(deptht) ;
		deptht:long_name = "Vertical T levels" ;
		deptht:units = "m" ;
		deptht:axis = "Z" ;
		deptht:positive = "down" ;
		deptht:bounds = "deptht_bounds" ;
	float deptht_bounds(deptht, axis_nbounds) ;