# CDO horizontal interpolation.

1. Read the [CDO documentation](https://code.mpimet.mpg.de/projects/cdo/).

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) ;
```