% Demo du KRIGING TOOLBOX
clc
echo on
%
% Demo du kriging toolbox qui utilise les fonctions 
% vario2di, fitvario, confint et cokri pour interpoler
% objectivement la salinité à 15 m sous la surface de 
% l'estuaire maritime du St-Laurent pendant COUPPB-90.
%
% La distribution de la salinité est considérée isotrope
% et l'interpolation d'un point utilise dépend de tous
% les points connus.

% Trace la carte des stations d'échantillonnage de G3

  figure('units','normal','Position',[0.5 0.5 0.5 1])
  plt90stn
  text(0.5,1.2,'Localisation des échantillons','Units',...
	'normalized','HorizontalAlignment','center','fontsize',14,...
	'color','b')
  pause
     
% Load les données de salinité à 15 m

  load sal15
 
% Calcule le variogramme de la salinité
 % 1- Arrangement des données

  x = sal15(:,1)./1000;
  y = sal15(:,2)./1000;
  z = sal15(:,3);

 % 2- Paramètres de gamv2

  nd     = length(x);				
  nlag   = 22;		
  xlag   = 3.5; 	
  xltol  = 1.75;
  ndir   = 1;
  azm    = 0; 
  atol   = 90;
  bandw  = 80;
  nvarg  = 2;
  ivtail = [1;1];
  ivhead = [1;1];
  ivtype = [1;4];
  nvar   = 1;

 % 3- Variogramme

  [np,gam,hm,tm,hv,tv] =... 
		vario2di (nd,x,y,z,nlag,xlag,xltol,ndir,azm,...
		       atol,bandw,nvarg,ivtail,ivhead,ivtype,nvar);

% Trouve les distances du variogramme pour lesquelles le nombre de paires est
% supérieur à 20

  I = find(np(:,2) >= 20);

% Fit le variogramme avec le modèle no21 de la fonction variogr

  fitvario(21,gam(I,1:2),200,0.05)
  refresh
  xlabel('Distance (km)')
  ylabel('Variance (psu²)')
  title('Variogramme','fontsize',14,'color','b')

% Calcul et trace les intervalles de confiance de 95%
 
  yesno = input('Désirez-vous voir les intervalles de confiance (''o''/''n'')? ');

  if yesno == 'o'
    for i = 1:length(I)
       [k1, k2] = confint(0.95,np(I(i),2),gam(I(i),2)); hold on
       plot(ones(1,2)*np(I(i),1),[k1,k2],'-g')
    end
  pause
  end  

% Krigeage
 % 1- paramètres de cokri.m
 
  mod = 11;		% modèle de variogramme 	
  a   = 79.1;		% paramètre a
  b   = 0.086*a;	% paramètre b (changement de variable de r à h)
  c   = 1.31;		% paramètre c

 % 2- Arrangement des données

  data = [sal15(:,1:2)./1000 sal15(:,3)];

 % 3- Position des points de la grille à interpoler

  [x0, nx, ny] = kregrid (150,2.5,222.5,5360,2.5,5430);

 % 4- Entrées de la fonction cokri

  model = [mod a];
  itype = 2;
  avg   = mean (data(:,3));
  block = [1 1];
  nd    = [1 1];
  ival  = 0;
  nk    = length(data);
  rad   = 80;
  ntok  = 1;

 % 5- Krigeage  

   yesno = 'n';

   if yesno == 'y'
     [x0s,s,sv,id,l] = cokri (data,x0,model,c,itype,avg,block,nd,ival,nk,rad,ntok,b);

     zi = deplie(x0s(:,3),nx,ny);
     szi = deplie(s(:,3),nx,ny); 
   else 
     load kri
   end

% Trace les résultats du krigeage
 % 1- Load la grille bathymétrique

  load bathy

 % 2- Masque les points de grille icorrect

  [l,n] = size(bathy);
  smin = 23; smax = 31; ds = 1; z = 15;

   for i = 1:l
	I = find (bathy(i,:) < z);
	zi(i,I)  = zi(i,I) .*nan;
        szi(i,I) = szi(i,I).*nan;
        bat(i,I) = bathy(i,I).*nan;
   end

 % 3- Trace le contour de la salinité

  xi = 150:2.5:222.5;
  yi = 5360:2.5:5430;
  
  clf  
  colormap (gray(8))
  pcolor (xi*1000,yi*1000,zi);
  shading interp
  caxis([smin,smax])
  text(0.5,1.2,'Salinité à 15 m','Units','normalized',...
	'HorizontalAlignment','center','fontsize',14,...
	'color','b')
  hold on

 % 4- Trace la côte

  plt90stn

 % 5- Trace l'échelle de couleur

  subplot('Position',[0.3 0.15 0.4 0.04])
  colormap (gray(8))
  map = pcolor([smin:ds:smax],[1,2],[smin:ds:smax;smin:ds:smax]);
  shading flat
  set (gca,'yticklabels',[],'xlim',[smin+1,smax-1],...
	'xtick',smin+1:ds:smax-1,'Fontsize',9)
  pause

 % 6- Trace le contour de l'erreur = sqrt(variance)

  clf
  cont = contour(xi.*1000,yi.*1000,sqrt(szi),(0.2:0.2:3));
  legend('0.2 ','0.4 ','0.6 ','0.8 ','1.0 ',-1)
  hold on
  plt90stn
  text(0.5,1.2,'Contour de l''erreur d''interpolation','color','b',...
	'Units','normalized','HorizontalAlignment','center','fontsize',14)
  pause

% Fin du demo 

  close