% covEstimationS            - covariance estimation tutorial using spatial data

% Generate, if needed, the data files for the spatial fields 
% ----------------------------------------------------------
echo off;

if exist('covexp.txt')~=2 | exist('covgau.txt')~=2 | exist('covnuggaus.txt')~=2 | exist('covanis.txt')~=2,
  %
  % Generate the values for exponential field and save it to file
  %
  rand('state',2);
  randn('state',2);
  nSim=1;
  ch=20*rand(100,2);
  covmodel='exponentialC';
  covparam=[1.5 10];
  [zh,L]=simuchol(ch,covmodel,covparam,nSim);
  zh=zh-min(zh)+1;
  varNames={'x','y','zh'};
  writeGeoEAS([ch zh],varNames,'Hard data with exponential covariance','covexp.txt');
  %
  % Generate the values for gaussian field and save it to file
  %
  nSim=1;
  ch=20*rand(100,2);
  covmodel='gaussianC';
  covparam=[1.5 10];
  [zh,L]=simuchol(ch,covmodel,covparam,nSim);
  zh=zh-min(zh)+1;
  varNames={'x','y','zh'};
  writeGeoEAS([ch zh],varNames,'Hard data with gaussian covariance','covgau.txt');
  %
  % Generate the values for gaussian field with nugget effect and save it to file
  %
  ch=20*rand(100,2);
  covmodel={'nuggetC','gaussianC'};
  covparam={[0.5],[1.0 20]};
  [zh,L]=simuchol(ch,covmodel,covparam,nSim);
  zh=zh-min(zh)+1;
  varNames={'x','y','zh'};
  writeGeoEAS([ch zh],varNames,'Hard data with nugget/gaussian cov','covnuggau.txt');
  %
  % Generate the values for anisotropic field and save it to file
  %
  rand('state',2);
  randn('state',2);
  ch=20*rand(200,2);
  covmodel='exponentialC';
  covparam=[1.5 10];
  angle=45;           % counterclockwise angle in degrees between 
                      % east and principal direction
  ratio=4;            % range in major direction divided by the  
                      % range in principal dir (>1)
  chiso=aniso2iso(ch,angle,ratio);
  [zh,L]=simuchol(chiso,covmodel,covparam,nSim);
  zh=zh-min(zh)+1;
  writeGeoEAS([ch zh],varNames,'Hard with anisotric covariance','covanis.txt');
end;

%%% Clear memory content and set echo to on.

clear;
clc;
echo on;

%%% ESTIMATION OF AN EXPONENTIAL COVARIANCE
%%% ---------------------------------------

%%% Read the hard data from file, and display it
%%% using circles having a size which is proportional to 
%%% value of each hard data.

[val,valname,filetitle]=readGeoEAS('covexp.txt');
ch=val(:,1:2);
zh=val(:,3);

figure;hold on;
markerplot(ch,zh,[4 30]);
xlabel('x (Km)');
ylabel('y (Km)');
title('Markerplot of data');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in all
%%% directions, and plot the result on a new graphic
%%% window

figure;
hold on;
cl=[0 1.5 4 7 20];
[d,C,o]=covario(ch,zh,cl,'kron');
variance=var(zh);
h(1)=plot([0 d'],[variance C'],'.r:');
plot([0 15],[0 0],'k-');
title('Covariance C(r)');

%%% Type any key to continuing...
pause;
clc;

%%% A reasonable model is an exponential covariance
%%% function with a sill (variance) cc=1.00,
%%% and a range aa=4.5 Km. 

r=0:.2:15;
covmodel='exponentialC';
cc=1.00;
aa=4.5;
covparam=[cc aa];
modelplot(r,covmodel,covparam);

%%% Type any key to continuing...
pause;
clc;

%%% ESTIMATION OF A GAUSSIAN COVARIANCE
%%% -----------------------------------

%%% Read the hard data from file, and display it
%%% using circles having a size which is proportional to 
%%% value of each hard data.

[val,valname,filetitle]=readGeoEAS('covgau.txt');
ch=val(:,1:2);
zh=val(:,3);

figure;hold on;
markerplot(ch,zh,[4 30]);
xlabel('x (Km)');
ylabel('y (Km)');
title('Markerplot of data');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in all
%%% directions, and plot the result on a new graphic
%%% window

figure;
hold on;
cl=[0 2 4 6 8 10 12 14 16];
[d,C,o]=covario(ch,zh,cl,'kron');
variance=var(zh);
h(1)=plot([0 d'],[variance C'],'.r:');
plot([0 15],[0 0],'k-');
title('Covariance C(r)');

%%% Type any key to continuing...
pause;
clc;

%%% A reasonable model is an gaussian covariance
%%% function with a sill (variance) cc=1.7, 
%%% and a range aa=6.5 Km. 

r=0:.2:15;
covmodel='gaussianC';
covparam=[1.7 6.5];
modelplot(r,covmodel,covparam);
xlabel('spatial lag r (Km)');
ylabel('Covariance C(r)');

%%% Type any key to continuing...
pause;
clc;

%%% ESTIMATION OF A GAUSSIAN COVARIANCE WITH NUGGET EFFECT
%%% ------------------------------------------------------

%%% Read the hard data from file, and display it
%%% using circles having a size which is proportional to 
%%% value of each hard data.

[val,valname,filetitle]=readGeoEAS('covnuggau.txt');
ch=val(:,1:2);
zh=val(:,3);

figure;hold on;
markerplot(ch,zh,[4 30]);
xlabel('x (Km)');
ylabel('y (Km)');
title('Markerplot of data');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in all
%%% directions, and plot the result on a new graphic
%%% window

figure;
hold on;
cl=[0 1.0 2.2  4 6  11];
[d,C,o]=covario(ch,zh,cl,'kron');
variance=var(zh);
h(1)=plot([0 d'],[variance C'],'.r:');
plot([0 15],[0 0],'k-');
title('Covariance C(r)');

%%% Type any key to continuing...
pause;
clc;

%%% A reasonable model is an gaussian covariance
%%% function with a sill (variance) cc=0.22
%%% and a range aa=10 Km, and 
%%% a nugget effect with sill 0.43

r=0:.2:15;
covmodel={'nuggetC','gaussianC'};
covparam={[0.43],[0.22 10]};
modelplot(r,covmodel,covparam);

%%% Type any key to continuing...
pause;
clc;

%%% ESTIMATION OF AN ANISOTROPIC COVARIANCE
%%% ---------------------------------------

%%% Read the hard data from file, and display it
%%% using circles having a size which is proportional to 
%%% value of each hard data.

clear;
[val,valname,filetitle]=readGeoEAS('covanis.txt');
ch=val(:,1:2);
zh=val(:,3);

figure;hold on;
markerplot(ch,zh,[4 30]);
xlabel('x (Km)');
ylabel('y (Km)');
title('Markerplot of data');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in all
%%% directions, and plot the result on a new graphic
%%% window

figure;
hold on;
cl=[0 2 4 6 15];
[d,C,o]=covario(ch,zh,cl,'kron');
variance=var(zh);
h(1)=plot([0 d'],[variance C'],'.r-');
plot([0 15],[0 0],'k-');
title('Covariance C(r)');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in the 
%%% direction with an angle=45 degree clockwise
%%% for the horizontal direction (with a tolerance 
%%% of +/-45 degrees), and plot the result in a same
%%% graphic window

cl=[0 2 4 8 15];
[d,C,o]=covario(ch,zh,cl,'kron',[0 0 90]);
h(2)=plot([0 d'],[variance C'],'.r--');
legend(h,'all directions','45 degree');

%%% Type any key to continuing...
pause;
clc;

%%% Estimate the covariance by using pairs in the 
%%% direction with an angle=-45 degree clockwise
%%% for the horizontal direction (with a tolerance 
%%% of +/-45 degrees), and plot the result in a same
%%% graphic window

cl=[0 1 2 4 15];
[d,C,o]=covario(ch,zh,cl,'kron',[0 -90 0]);
h(3)=plot([0 d'],[variance C'],'.r:');
legend(h,'all directions','45 degree','-45 degrees');

%%% Type any key to continuing...
pause;
clc;

%%% A reasonable model is an exponential covariance
%%% function with a sill (variance) cc=1.55,
%%% and a range aa=7 Km is the major direction at 45 degree,
%%% and a range aa=3 Km in the minor direction at -45 degree,
%%% which corresponds to an anisotropy ratio of 7/3=2.3

r=0:.2:15;
covmodel='exponentialC';
cc=1.55;

aa=7;
covparam=[cc aa];
modelplot(r,covmodel,covparam);

aa=3;
covparam=[cc aa];
modelplot(r,covmodel,covparam);