interpolateNetCDF

 Function expecting an EMODnet Bathymetry-compliant NetCDF4 file
 Input:
   - method: the method to use. Available:
               - 'Nearest': Nearest Neighbor.
               - 'Delaunay': Delaunay triangulation (linear interpolation).
               - 'Natural': Natural neighbors.
               - 'IDW': Inverse Distance Weighted.
               - 'Kriging': Kriging interpolation.
               - 'MLS': Moving Least Squares.
               - 'RBF.<rbf_type>': Radial Basis Functions.
               - 'QTPURBF.<rbf_type>': QuadTree Partition of Unity BRF.
                   For the last two cases, you need to indicate in <rbf_type>
                   the type of the Radial Basis Function to use. Available:
                       - 'linear'
                       - 'cubic'
                       - 'quintic'
                       - 'multiquadric'
                       - 'thinplate'
                       - 'green'
                       - 'tensionspline'
                       - 'regularizedspline'
                       - 'gaussian'
                       - 'wendland'
               (For more information on each RBF, please check their
               individual documentation in their corresponding functions on
               the "rbf" folder of this toolbox)
               - 'inpainting.<type>': Inpainting method. Available types:
                   - 'sobolev'
                   - 'tv'
                   - 'amle'
                   - 'ccst'
                   - 'bertalmio'

0001 function interp = interpolateNetCDF(inputNetCDF, outputNetCDF, method, options)
0002 % Function expecting an EMODnet Bathymetry-compliant NetCDF4 file
0003 % Input:
0004 %   - method: the method to use. Available:
0005 %               - 'Nearest': Nearest Neighbor.
0006 %               - 'Delaunay': Delaunay triangulation (linear interpolation).
0007 %               - 'Natural': Natural neighbors.
0008 %               - 'IDW': Inverse Distance Weighted.
0009 %               - 'Kriging': Kriging interpolation.
0010 %               - 'MLS': Moving Least Squares.
0011 %               - 'RBF.<rbf_type>': Radial Basis Functions.
0012 %               - 'QTPURBF.<rbf_type>': QuadTree Partition of Unity BRF.
0013 %                   For the last two cases, you need to indicate in <rbf_type>
0014 %                   the type of the Radial Basis Function to use. Available:
0015 %                       - 'linear'
0016 %                       - 'cubic'
0017 %                       - 'quintic'
0018 %                       - 'multiquadric'
0019 %                       - 'thinplate'
0020 %                       - 'green'
0021 %                       - 'tensionspline'
0022 %                       - 'regularizedspline'
0023 %                       - 'gaussian'
0024 %                       - 'wendland'
0025 %               (For more information on each RBF, please check their
0026 %               individual documentation in their corresponding functions on
0027 %               the "rbf" folder of this toolbox)
0028 %               - 'inpainting.<type>': Inpainting method. Available types:
0029 %                   - 'sobolev'
0030 %                   - 'tv'
0031 %                   - 'amle'
0032 %                   - 'ccst'
0033 %                   - 'bertalmio'
0034 
0035 %% Script
0036 img = ncread(inputNetCDF , 'elevation');
0037 mask = ncread(inputNetCDF, 'interpolation_flag');
0038 mask = ~isnan(mask); % Convert to boolean
0039 img(~mask) = 0; % To avoid numerical issues because of NaNs
0040 
0041 C = strsplit(method, '.');
0042 methodStart = C{1};
0043 if numel(C) > 1
0044     methodSub = lower(C{2});
0045 end
0046 
0047 if strcmpi(methodStart, 'nearest') || ...
0048    strcmpi(methodStart, 'delaunay') || ...
0049    strcmpi(methodStart, 'natural') || ...
0050    strcmpi(methodStart, 'idw') || ...
0051    strcmpi(methodStart, 'kriging') || ...
0052    strcmpi(methodStart, 'mls') || ...
0053    strcmpi(methodStart, 'rbf') || ...
0054    strcmpi(methodStart, 'qtpurbf')
0055 
0056     if nargin < 4 && ~strcmpi(methodStart, 'Nearest') && ~strcmpi(methodStart, 'Delaunay') && ~strcmpi(methodStart, 'Natural')
0057         options = hmitScatteredDefaultOptions(x, y, z, xi, yi);
0058     end
0059 
0060     % Scattered data interpolation method
0061     error('Still not implemented!');
0062     
0063 elseif strcmpi(methodStart, 'inpainting')
0064         
0065     if nargin < 4
0066         options = hmitInpaintingDefaultOptions(methodSub);
0067     end
0068     
0069     % Compute the grid step from the data and override defaults or user input
0070 %     lats = ncread(inputNetCDF , 'lat');
0071 %     lons = ncread(inputNetCDF , 'lon');
0072 %     hx = abs((lons(end)-lons(1))/numel(lons));
0073 %     hy = abs((lats(end)-lats(1))/numel(lats));
0074 %     options.GridStepX = hx;
0075 %     options.GridStepY = hy;
0076     
0077     % Inpainting method
0078     param = structToVarargin(options);
0079     switch methodSub
0080         case 'sobolev'            
0081             inpainter = SobolevInpainter(param{:});            
0082         case 'tv'
0083             inpainter = TVInpainter(param{:});
0084 %             interp = inpaintTV(img, mask, 1);
0085         case 'amle'
0086             inpainter = AMLEInpainter(param{:});
0087 %             interp = inpaintAMLE(img, mask);
0088         case 'ccst'
0089             inpainter = CCSTInpainter(param{:});
0090 %             interp = inpaintCCST(img, mask, 0);
0091         case 'bertalmio'
0092             inpainter = BertalmioInpainter(param{:});
0093         otherwise
0094             error('Unknown inpainting type');
0095     end    
0096     % Inpaint!
0097     interp = inpainter.inpaint(img, mask);
0098 else
0099     error(sprintf('Unknown method %s', method));
0100 end
0101 
0102 % Write results
0103 copyfile(inputNetCDF, outputNetCDF);
0104 ncwrite(outputNetCDF, 'elevation', interp);