% edit entries below this line to set the mapping options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option = 'griddata';
option = 'loess2dg';
option = 'oi';

switch option
  
  case 'griddata'
    % method = 'linear';
    method = 'cubic';
    
  case 'loess2dg'
    % The 'number' method searches for R points nearest to each grid location
    % and applies the loess smoother
    %
    % The 'distance' method searches for all the points within R (in units
    % of the input data, in this case degrees lon/lat) of each grid
    % location and applies the loess smoother
    
    loess_method = 'distance';
    R = 0.4; % degrees lon/lat
    
    %loess_method = 'number';
    %R = 400; % nearest points
    
  case 'oi'
    R = 0.5; % length scale in deg lon/lat of Gaussian covariance 
             % function
    lambda = 0.1; % error squared to signal squared

end

% edit entries above this line to set the mapping options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load the data from RV Cape Hatteras and RV Oceanus 
load latte2005_surface_ships
% if you have trouble loading the data the the other file
% latte2005_surface_ships_v7.mat

disp('Data are in the time range:')
disp(datestr(time([1 end])))
disp(' ')

% nominal date of the map
% You may try changing this to see how the interpolation methods treat 
% time dependence 
d = datenum(2005,4,18);

disp(['Analysis date will be ' datestr(d)])
disp(' ')

% time data are in matlab datenum form
% subset +/- N days from nominal date of map
Ndays = 5;
s = find(abs(time-d)<=Ndays);
for varlist = {'lon','lat','sal','temp','time'}
  v = char(varlist);
  eval([v ' = ' v '(s);'])
end

% report the option being used
disp(['Using mapping option ' option])

% regular lon/lat target grid at whicyh to compute the map
x = -74.15:0.01:-73.55;
y = 40.01:0.01:40.58;
[X,Y] = meshgrid(x,y);

% map these data
data = sal; cax = [20 32]; % salinity

switch option
  
  case 'griddata'
    % map a map with griddata
    D = griddata(lon,lat,data,X,Y,method);

  case 'loess2dg'
    % Use John's 2-dimensional quadratic loess smoother (loess2dg.m)

    % the 'number' option searches for R points nearest to each grid location
    % and applies the loess smoother
    
    % the 'distance' option searches for all the points within R (in units
    % of the input data, in this case degrees lon/lat) of each grid
    % location and applies the loess smoother
    disp(['with method ' loess_method])
    disp(['and length or number scale R = ' num2str(R)])
    
    D = loess2dg(lon,lat,data,X,Y,R,loess_method); 

  case 'oi'
    
     disp(['length scale ' num2str(R) ' lon/lat'])
     disp(['and error to signal variance ratio lambda of ' num2str(lambda)])
    % optimal interpolation
    % R = length scale in degrees
    
    xr = repmat(lon(:),[1 length(lon)]);
    yr = repmat(lat(:)',[length(lon) 1]);
    Rdd = sqrt((xr-xr').^2+(yr-yr').^2);
    Cdd0 = exp(-Rdd/R^2);
    
    % Add error/signal variance to the diagonal
    Cdd = Cdd0 + lambda*eye(size(Cdd0));
    
    xr = repmat(lon(:),[1 length(X(:))]);
    yr = repmat(lat(:),[1 length(X(:))]);
    xg = repmat(X(:),[1 length(lon)]);
    yg = repmat(Y(:),[1 length(lon)]);
    
    Rmd = sqrt((xr-xg').^2+(yr-yg').^2);
    Cmd = exp(-Rmd/R^2);
    Cmd=Cmd';
    
    data=data(:);
    D = mean(data) + Cmd*(Cdd\(data-mean(data)));
    D = reshape(D,size(X));
    
  otherwise
    error('Oops: you have not set the OPTION for mapping method')
    
end

% plot the results
pcolor(X,Y,D)
shading flat
axis([-74.2 -73.5 40 40.65])
caxis(cax);
colorbar
hold on
han = scatter(lon(:),lat(:),30,data(:),'filled');
set(han,'markeredgecolor','k')

% plot the coastline so we recognize where we are
plotmabcoast
hold off


% interpolate solution to the data points to examine the residuals
data_m = interp2(X,Y,D,lon,lat);
residuals = data(:) - data_m(:);
rmse = std(residuals(find(isnan(residuals)~=1)));
disp(' ')
disp(['Mean squared error of mapping was ' num2str(rmse)])
disp(' ')