Octave script to calculate the spatial semivariogram.

semivariogram.m

function [s lags ncp model]=semivariogram(pos,z,delta,varargin)
% Calculates an empirical spatial semivariogram
%
% SYNTAX:
%	[s lags ncp]=semivariogram(pos,z,delta,varargin)
%
% INPUTS:
%	Mandatory:
%	- pos: (NxD) matrix containing the D dimensional position of the control points.
%	- z: (Nx1) vector containing the observations at the control points.
%	- delta: (scalar) the width of the lag-intervals where the semi-variance will be calculated
%
%	Optional:
%	- MaxLag: maximum lag for which the semivariogram will be calculated. This will be adjusted
%		depending on the value of delta so that MaxLag is included in the plot.
%		Typically for large relative lags (getting close to the maximum) there will be less
%		control point pairs available, and therefore the semi-variance for those lags will be poorly
%		estimated. If not specified, a maxim lag equal to half of the maximum distance among control
%		points will be used. This is a "rule of thumb" distance of reliability (gstat uses 1/3 of the
%		diagonal of the bounding box containing all the control points).
%		If you want to use all the possible lags, use a large number (e.g. Inf)
%	- Method: string that can be either (where lower and upper case are admitted)
%		- Matheron: for the classic empirical semi-variogram using moments (Matheron, see  [1])
%		- Dowd: median absolute deviation.(Rousseeuw, P.J. and Croux, see  [4])
%		- CH: for the robust semi-variogram (Cressie and Hawkins, see [2])
%		- QN: Genton's highly robust semi-variogram (Genton, see [3])
%
%	Param-value pairs:
%	- 'FitModel': either a string or a cell of strings specifying the type of models to fit to the data.
%		Options:
%		-'Exponential'
%			         |0 if h=0
%			 gamma= < c_n + var_0*(3h/(2r) - 1/2*(h/r)^3)
%			         |c_0 if h>r
%		-'Gaussian'
%		       		  |0 if h=0
%			 gamma= < c_n + var_0*(3h/(2r) - 1/2*(h/r)^3)
%			         |c_0 if h>r
%		-'Spherical'
%			         |0 if h=0
%			 gamma= < c_n + var_0*(3h/(2r) - 1/2*(h/r)^3)
%			         |c_0 if h>r
%		where
%		 -c_n is the nugget
%		 -var_0 is the partial sill, so that c_0=c_n+var_0
%		 -r is the range
%		 -the sill c_0 = c_n + var_0
%	- 'FixParams': Can be used to fix the variogram model parameters to a certain value. If specified and if only one model
%		is to be fitted (FitModel is a string), it must be a vector of size 3 whose elements describe the fixed value for
%		c_n, var_0 and r in that order. Setting the value to NA means that the value is not fixed.
%		If multiple models are to be fitted (FitModel is a cell array), FixParams must be a cell array containing one
%		vector for each model specified in FitModel.
%	- 'FigNum': if specified, a plot of the semi-variogram will be produced, in the window
%		number specified by FigNum. If the user has requested to fit a model, the model will be also plotted.
%
% OUTPUTS
%	- s: (Nx1) vector containing values of the semi-variance for the lags
%	- lags:(Nx1) vector containing the lags for which the semi-variance was calculated. These are the center
%		of the intervals defined by 'delta'.
%	- ncp: the number of control points used to calculate the average semi-variance in each lag.
% 		The precission of the semi-variance estimates for each lag will depend on the number of
% 		samples used by the estimator. These estimators might be used later for other estimation,
%		and so the different precissions will have to be taken into account (heteroskedasticity).
%	- model: structure (or cell array of structures) containing the parameters of the model(s) fitted to the data (see 'FitModel' param-value pair).
%		The fields are:
%		- type: string containing the type of model fitted.
%		- c_n: nugget
%		- var_0: partial sill 
%		- r: the range
%		- g: semi-variance calculated with the fitted model for the lags given in the output
%
% EXAMPLES
%	[s lags np]=semivariogram(pos,z, 6,105);
%	[s lags np]=semivariogram(pos,z, 6,105,'CH');
%	[s lags np]=semivariogram(pos,z, 6,105,'Dowd','FigNum',1);
%	[s lags np model]=semivariogram(pos,z, 6,105,'QN','FigNum',1,'FitModel','gaussian');
%	[s lags np model]=semivariogram(pos,z, 6,105,'QN','FigNum',1,'FitModel',{'gaussian','exponential'});
%	[s lags np model]=semivariogram(pos,z, 6,105,'QN','FigNum',1,'FitModel',{'gaussian','exponential'},'FixParams',{[NA NA NA],[0 NA NA]}); %fix the nugget to zero in the exponential
%
% REFERENCES
% [1] - Matheron, G. (1963), “Principles of Geostatistics,” Economic Geology, 58, 1246–1266.
% [2] - Cressie, N. and Hawkins, D. M. (1980), “Robust Estimation of the Variogram: I,”
%	Mathematical Geology, 12(2), 115–125.
% [3] - Genton, M.G. (1998), "Highly Robust Variogram Estimation", Mathematical Geology,v. 30, n.2, p.213-221, 
% [4] - Rousseeuw, P.J. and Croux, C. (1993). ”Alternatives to the Median Absolute Deviation,”
%	Journal of the American Statistical Association, Vol. 88, 1273-1283.
%


%
% NOTE the generalized use of broadcasting operators!
%

if nargin < 3
  error("Expecting minimum 3 arguments.");
endif

opts=validate_input_args(pos,z,delta,varargin{:});

if strcmpi(opts.Method,"Matheron")
  Y=triu((z-z').^2); %(Upper triangular) matrix of squared differences
elseif strcmpi(opts.Method,"CH") || strcmpi(opts.Method,"Dowd")
  Y=triu(sqrt(abs(z-z'))); %(Upper triangular) matrix of squared root differences
elseif strcmpi(opts.Method,"QN")
  Y=triu(z-z'); %(Upper triangular) matrix of differences
endif

%(Upper triangular) matrix of distances among control points. D(n,k) = distance from n-th to k-th points
D=triu(permute(sqrt(sum((pos-permute(pos,[3 2 1])).^2,2)),[1 3 2]));

%Edges of bins and lags for each bin
if isna(opts.MaxLag)
  %opts.MaxLag=sqrt(sum((max(pos)-min(pos)).^2))/3; % Distance of reliability used in gstat
  opts.MaxLag=max(D(:))*1/2; 
else
  opts.MaxLag=min([opts.MaxLag max(D(:))]);
endif
bin_edges = (0:delta:opts.MaxLag+delta)';
lags=bin_edges(1:end-1)+delta/2;

s=zeros(size(lags));
ncp=zeros(size(lags));

%Calculate the average semi-covariance for each lag
for l_i=1:numel(lags) %lag index
  %Note that the > instead of >= in the first comparison discards the terms corresponding
  %to the diagonal and lower traingular elements of D
  idx=find(D>lags(l_i)-delta/2 & D<=lags(l_i)+delta/2);
  num_pairs=numel(idx);
  if num_pairs==0
    continue
  endif
  if strcmpi(opts.Method,"Matheron")
    s(l_i)=mean(Y(idx))/2;
  elseif strcmpi(opts.Method,"Dowd")
    Phi=median(Y(idx));
    s(l_i)= Phi^4 / 0.457 / 2;
  elseif strcmpi(opts.Method,"CH")
    Phi=mean(Y(idx));
    s(l_i)= Phi^4 / (0.457 + 0.494/num_pairs + 0.045/num_pairs^2) / 2;
  elseif strcmpi(opts.Method,"QN");
    if num_pairs==1 %There are no pairs of distances to compare
      continue;
    endif
    V=abs(triu(Y(idx)-Y(idx)')); %Set of all absolute diferences for i<j. This can be a very LARGE matrix!
    idx2=find(triu(ones(size(V)),1));
    V=sort(V(idx2));
    k=bincoeff( floor(num_pairs/2)+1 , 2 );
    s(l_i)=( (2.2191*V(k))^2 ) / 2; %The k-th order statistic of the interpoint distances, that is, the k-th smallest value.
    lags(l_i)=mean(D(idx));
  endif
  ncp(l_i)=num_pairs;
  samples_i{l_i}=idx;
endfor
e_i=find(ncp==0); %erase index
s(e_i)=[];
lags(e_i)=[];
ncp(e_i)=[];
%If the last bin contains no data, samples_i will have one element less than s or lags, as adding an empty vector to a cell
%does not effectively increase the cell's size. So we avoid the size reduction of samples_i in that case.
e_i(find(e_i==l_i))=[];
samples_i(e_i)=[];

%Fit models
model={};
if !iscell(opts.FitModel)
  opts.FitModel={opts.FitModel};
  opts.FixParams={opts.FixParams};
endif
for m_i=1:numel(opts.FitModel)
  if strcmpi(opts.FitModel{m_i},"None")
    continue;
  elseif strcmpi(opts.FitModel{m_i},"Spherical")
    model{m_i}=fit_spherical_model(lags,s,opts.FixParams{m_i});
  elseif strcmpi(opts.FitModel{m_i},"Exponential")
    model{m_i}=fit_exponential_model(lags,s,opts.FixParams{m_i});
  elseif strcmpi(opts.FitModel,"Gaussian")
    model{m_i}=fit_gaussian_model(lags,s,opts.FixParams{m_i});
  endif
endfor

%Plots
if !isna(opts.FigNum)
  figure(opts.FigNum);
  clf;
  hold on;
  %Plot the models
  if numel(model)!=0
    %col=hsv(numel(model)); %Is this correct in octave?
    N=numel(model);col=hsv2rgb ([linspace(0,1,N)'/N, ones(N,2)]);
    legnd={};
    for m_i=1:numel(model)
      if !strcmpi(model{m_i},"None")
	nugget=model{m_i}.c_n;
	plot([0;lags],[nugget;model{m_i}.g],'Color',col(m_i,:));
	legnd{end+1}=model{m_i}.type;
      endif
    endfor
    legend(legnd);
  endif
  %Plot the estimated semi-variances
  plot(lags,s,'o');
  grid on;
  xlabel('Lag');
  ylabel('Semi-variance');
  %{
  %Plot the squared differences
  if !strcmpi(opts.Method,"Matheron")
    Y=triu((z-z').^2); %(Upper triangular) matrix of squared differences
  endif
  for l_i=1:numel(lags)
    idx=samples_i{l_i};
    plot(D(idx),Y(idx),'.');
  endfor
  %}

endif

if numel(model)==1
  model=model{1};
endif

%%%%%%%%%%%%%%%%%%%%%%%
%
% Helper functions to fit the semi-variogram models.
%
% The form of the equations have been taken from 
% http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_variogram_a0000000579.htm
%
%%%%%%%%%%%%%%%%%%%%%%%


%Fit a spherical model, defined as
%         |0 if h=0
% gamma= < c_n + var_0*(3h/(2r) - 1/2*(h/r)^3)
%         |c_0 if h>r
% where
% 	-c_n is the nugget
%	-var_0 is the partial sill, so that c_0=c_n+var_0
% 	-r is the range
%
% The sill c_0 = c_n + var_0 
% 
function model=fit_spherical_model(lags,s,FixParams)
  model.type="Spherical";
  %Initial guess of the parameters
  stats_s=statistics(s);
  theta_0=[stats_s(1) stats_s(4) mean(lags)];
  %gamma function handle. Carefully ensure that gamma=c_0 for h>r
  gf_h=@(theta) [(theta(1) + theta(2)*(3*lags/(2*theta(3)) - (lags/theta(3)).^3/2))(find(lags<=theta(3))); repmat(theta(1)+theta(2),[rows(lags)-numel(find(lags<=theta(3))),1]) ];
  cost=@(theta) sum((s-gf_h(theta)).^2);
  %Do we have fixed parameters?
  fp_i=find(!isna(FixParams)); %index of fixed parameters
  fp_val=FixParams(fp_i); %value of the fixed parameters
  theta_min=[0 0 0];
  theta_min(fp_i)=fp_val;
  theta_max=[Inf Inf Inf];
  theta_max(fp_i)=fp_val;
  %Constrained optimization
  [theta,sse,exitflag] = sqp(theta_0,cost,[],[],theta_min,theta_max);
  if exitflag!=101 && exitflag!=103
    warning(sprintf("semivariogram: BFGS update failed when fitting %s model.",model.type));
  endif
  model.c_n=theta(1);
  model.var_0=theta(2);
  model.r=theta(3);
  model.g=gf_h(theta);
  model.sse=sse;
  model.exitflag=exitflag;
  

  
%Fit an exponential model, defined as
%         |0 if h=0
% gamma= < c_n + var_0(1-e^(-h/r))
%         \
% where
% 	-c_n is the nugget
%	-var_0 is the partial sill, so that c_0=c_n+var_0
% 	-r is the range
%
% The sill c_0 = c_n + var_0 
%
function model=fit_exponential_model(lags,s,FixParams)
%{
  %We use an optimizer to find a solution using the initial values of the parameters
  %Find an initial guess usin LSE and a linearized gamma
  % gamma = c_n  + var_0/r h - var_0/r^2 h^2 + ...
  % gamma = x(1) +   x(2)  h -   x(3)    h^2
  A=[ones(size(lags)) lags -lags.^2 lags.^3];
  x=A\s;
  c_n=x(1);
  r=x(2)/x(3);
  var_0=x(2)^2/x(3);
%}

  %Use the spherical model parameters as an initial guess. This seems more robust than the
  %previous (commented) method of using LSE on a linear approximation of the model.
  model=fit_spherical_model(lags,s,FixParams);
  theta_0=[model.c_n model.var_0 model.r];

  model.type="Exponential";
  gf_h=@(theta) theta(1) + theta(2)*(1-exp(-lags/theta(3))); %gamma function handle
  cost=@(theta) sum((s-gf_h(theta)).^2);
  %Do we have fixed parameters?
  fp_i=find(!isna(FixParams)); %index of fixed parameters
  fp_val=FixParams(fp_i); %value of the fixed parameters
  theta_min=[0 0 0];
  theta_min(fp_i)=fp_val;
  theta_max=[Inf Inf Inf];
  theta_max(fp_i)=fp_val;
  %Constrained optimization
  [theta,sse,exitflag] = sqp(theta_0,cost,[],[],theta_min,theta_max);
  if exitflag!=101 && exitflag!=103
        warning(sprintf("semivariogram: BFGS update failed when fitting %s model.",model.type));
  endif

  model.c_n=theta(1);
  model.var_0=theta(2);
  model.r=theta(3);
  model.g=gf_h(theta);
  model.sse=sse;
  model.exitflag=exitflag;

%Fit a gaussian model, defined as
%         |0 if h=0
% gamma= < c_n + var_0(1-e^(-(h/r)^2)
%         \
% where
% 	-c_n is the nugget
%	-var_0 is the partial sill, so that c_0=c_n+var_0
% 	-r is the range
%
% The sill c_0 = c_n + var_0
%
function model=fit_gaussian_model(lags,s,FixParams)
  model.type="Gaussian";
  
  %We use an optimizer to find a solution using the initial values of the parameters
  %Find an initial guess usin LSE and a linearized gamma
  % gamma = c_n  + 2*var_0/r h^2 + 2*var_0/r^2 h^3 + ...
  % gamma = x(1) +    x(2)   h^2 +     x(3)    h^2
  A=[ones(size(lags)) lags -lags.^2 lags.^3];
  %x=A\s;
  x=lsqnonneg(A,s);
  c_n=x(1);
  r=x(2)/x(3);
  var_0=x(2)^2/x(3)/2;
  theta_0=[c_n var_0 r];
  
  gf_h=@(theta) theta(1) + theta(2)*(1-exp(-(lags/theta(3)).^2)); %gamma function handle
  cost=@(theta) sum((s-gf_h(theta)).^2);
  %Do we have fixed parameters?
  fp_i=find(!isna(FixParams)); %index of fixed parameters
  fp_val=FixParams(fp_i); %value of the fixed parameters
  theta_min=[0 0 0];
  theta_min(fp_i)=fp_val;
  theta_max=[Inf Inf Inf];
  theta_max(fp_i)=fp_val;
  %Constrained optimization
  [theta,sse,exitflag] = sqp(theta_0,cost,[],[],theta_min,theta_max);
  if exitflag!=101 && exitflag!=103
        warning(sprintf("semivariogram: BFGS update failed when fitting %s model.",model.type));
  endif

  model.c_n=theta(1);
  model.var_0=theta(2);
  model.r=theta(3);
  model.g=gf_h(theta);
  model.sse=sse;
  model.exitflag=exitflag;
  
% Validation of input arguments
function opts=validate_input_args(pos,z,delta,varargin);
  p = inputParser;
  p.FunctionName = "semivariogram";

  %REQUIRED parameters in order (must be before optionals)
  p = p.addRequired ("pos", @(x)(ismatrix(x)&&(columns(x)==2)));
  p = p.addRequired("z", @(x)iscolumn(x));  
  if(rows(z)!=rows(pos))
    error("Rows of z must be the same as rows of pos");
  endif
  p = p.addRequired("delta", @(x)isscalar(x)&&isreal(x));

  %%% OPTIONAL parameters
  p = addOptional(p,"MaxLag", NA, @(x)isscalar(x)&&isreal(x));
  methods={
    "Matheron",...
    "Dowd",...
    "CH",... 
    "QN"};
  p = addOptional(p,"Method", 'Matheron', @(x) ischar(x) && @any(strcmpi(x,methods)));
  
  %%% PARAM VALUE pairs
  p = addParamValue(p,"FigNum",NA ,@(x)isscalar(x)&&isreal(x));
  p = addParamValue(p,"FitModel","None", @validate_models);
  p = addParamValue(p,"FixParams",[NA NA NA], @validate_fixparams);
  %Parse the arguments
  p = p.parse (pos,z,delta,varargin{:});
  opts = p.Results;
  
  if (ischar(opts.FitModel) && iscell(opts.FixParams))  || ...
      iscell(opts.FitModel) && iscell(opts.FixParams) && numel(opts.FitModel)!=numel(opts.FixParams)
    error("semivariogram: FitModel and FixParams must both be either vectors or cells with the same number of elements.")
  endif
  if iscell(opts.FitModel) && !iscell(opts.FixParams) %FixParams is not specified and is the default
    opts.FixParams={opts.FixParams};
    opts.FixParams=repmat(opts.FixParams,[numel(opts.FitModel) 1]);
  endif

function isvalid=validate_models(x)
  models={
  "None",...
  "Spherical",...
  "Exponential",...
  "Gaussian"};

  if ischar(x)
    isvalid = @any(strcmpi(x,models));
  elseif iscell(x)
    for m_i=1:numel(x)
      if !any(strcmpi(x{m_i},models))
	isvalid=false;
       return;
      endif
    endfor
    isvalid=true;
  endif
                       
function isvalid=validate_fixparams(x)
  if ismatrix(x)
    if !isreal(x) || length(x)!=3
      isvalid=false;
      return;
    endif
  elseif iscell(x)
    for m_i=1:numel(x)
      if !isreal(x{m_i}) || size(x{m_i})!=[1 3]
	isvalid=false;
	return;
      endif
    endfor
  endif
  isvalid=true;
  
  
%!test
%! z=[ones(30,1) zeros(30,1)]'(:); %A sequence of 1-0-1-...
%! x=[1:length(z)]';
%! y=zeros(size(z));
%! pos=[x y];
%! [s lags np]=semivariogram(pos,z,1,Inf);
%! assert(s,0.5*z(1:end-1))


%!demo
%! %--------------------------------------------------
%! % Comparison with the variogram from
%! % http://www.ats.ucla.edu/stat/r/faq/variogram.htm
%! %--------------------------------------------------
%
%! f_name="ozone.csv";
%! while exist(f_name,'file') f_name=['_' f_name]; endwhile %Ensure that we don't overwrite an existing file named 'ozone.csv'
%! urlwrite("http://www.ats.ucla.edu/stat/r/faq/ozone.csv",f_name); %Download test data
%! system(["sed -i 's/\"/%\"/' " f_name]); %Comment the first line, which causes problems when loading the data
%! ozone=load(f_name);
%! system(['rm ' f_name]);
%! pos=[ozone(:,4) ozone(:,3)];
%! z=ozone(:,2);
%! [s lags np]=semivariogram(pos,z, 0.15,1.5,'FigNum',figure());
%! title("Ozone semi-variogram");
%


%!demo
%! %--------------------------------------------------
%! % Comparison with the variogram from
%! % http://www.ats.ucla.edu/stat/sas/faq/SAS_variogram_fit.htm
%! %--------------------------------------------------
%
%! data=[
%!   0.7  59.6  34.1   2.1  82.7  42.2   4.7  75.1  39.5
%!   4.8  52.8  34.3   5.9  67.1  37.0   6.0  35.7  35.9
%!   6.4  33.7  36.4   7.0  46.7  34.6   8.2  40.1  35.4
%!  13.3   0.6  44.7  13.3  68.2  37.8  13.4  31.3  37.8
%!  17.8   6.9  43.9  20.1  66.3  37.7  22.7  87.6  42.8
%!  23.0  93.9  43.6  24.3  73.0  39.3  24.8  15.1  42.3
%!  24.8  26.3  39.7  26.4  58.0  36.9  26.9  65.0  37.8
%!  27.7  83.3  41.8  27.9  90.8  43.3  29.1  47.9  36.7
%!  29.5  89.4  43.0  30.1   6.1  43.6  30.8  12.1  42.8
%!  32.7  40.2  37.5  34.8   8.1  43.3  35.3  32.0  38.8
%!  37.0  70.3  39.2  38.2  77.9  40.7  38.9  23.3  40.5
%!  39.4  82.5  41.4  43.0   4.7  43.3  43.7   7.6  43.1
%!  46.4  84.1  41.5  46.7  10.6  42.6  49.9  22.1  40.7
%!  51.0  88.8  42.0  52.8  68.9  39.3  52.9  32.7  39.2
%!  55.5  92.9  42.2  56.0   1.6  42.7  60.6  75.2  40.1
%!  62.1  26.6  40.1  63.0  12.7  41.8  69.0  75.6  40.1
%!  70.5  83.7  40.9  70.9  11.0  41.7  71.5  29.5  39.8
%!  78.1  45.5  38.7  78.2   9.1  41.7  78.4  20.0  40.8
%!  80.5  55.9  38.7  81.1  51.0  38.6  83.8   7.9  41.6
%!  84.5  11.0  41.5  85.2  67.3  39.4  85.5  73.0  39.8
%!  86.7  70.4  39.6  87.2  55.7  38.8  88.1   0.0  41.6
%!  88.4  12.1  41.3  88.4  99.6  41.2  88.8  82.9  40.5
%!  88.9   6.2  41.5  90.6   7.0  41.5  90.7  49.6  38.9
%!  91.5  55.4  39.0  92.9  46.8  39.1  93.4  70.9  39.7
%!  94.8  71.5  39.7  96.2  84.3  40.3  98.2  58.2  39.5
%! ];
%! pos=[data(:,1:2);data(:,4:5);data(:,7:8)];
%! z=[data(:,3);data(:,6);data(:,9)];
%! [s lags np]=semivariogram(pos,z, 6,105,'matheron','FigNum',figure());
%! title('Matheron');
%! [s lags np]=semivariogram(pos,z, 6,105,'CH','FigNum',figure(),'FitModel',{'Exponential','Exponential'},'FixParams',{[NA NA NA],[NA NA 10]});
%! title('Cressie-Hawkins + Exponential fit + Exponential fit with fixed r=10');
%! [s lags np]=semivariogram(pos,z, 6,105,'QN','FigNum',figure(),'FitModel','Gaussian');
%! title('Genton (QN) + Gaussian fit');
%! [s lags np]=semivariogram(pos,z, 6,105,'Dowd','FigNum',figure(),'FitModel','Spherical');
%! title('Dowd + Spherical fit');
%

%{
data=[
 0.7  59.6  34.1   2.1  82.7  42.2   4.7  75.1  39.5 
 4.8  52.8  34.3   5.9  67.1  37.0   6.0  35.7  35.9
 6.4  33.7  36.4   7.0  46.7  34.6   8.2  40.1  35.4
 13.3   0.6  44.7  13.3  68.2  37.8  13.4  31.3  37.8
 17.8   6.9  43.9  20.1  66.3  37.7  22.7  87.6  42.8
 23.0  93.9  43.6  24.3  73.0  39.3  24.8  15.1  42.3
 24.8  26.3  39.7  26.4  58.0  36.9  26.9  65.0  37.8
 27.7  83.3  41.8  27.9  90.8  43.3  29.1  47.9  36.7
 29.5  89.4  43.0  30.1   6.1  43.6  30.8  12.1  42.8
 32.7  40.2  37.5  34.8   8.1  43.3  35.3  32.0  38.8
 37.0  70.3  39.2  38.2  77.9  40.7  38.9  23.3  40.5
 39.4  82.5  41.4  43.0   4.7  43.3  43.7   7.6  43.1
 46.4  84.1  41.5  46.7  10.6  42.6  49.9  22.1  40.7
 51.0  88.8  42.0  52.8  68.9  39.3  52.9  32.7  39.2
 55.5  92.9  42.2  56.0   1.6  42.7  60.6  75.2  40.1
 62.1  26.6  40.1  63.0  12.7  41.8  69.0  75.6  40.1
 70.5  83.7  40.9  70.9  11.0  41.7  71.5  29.5  39.8
 78.1  45.5  38.7  78.2   9.1  41.7  78.4  20.0  40.8
 80.5  55.9  38.7  81.1  51.0  38.6  83.8   7.9  41.6
 84.5  11.0  41.5  85.2  67.3  39.4  85.5  73.0  39.8
 86.7  70.4  39.6  87.2  55.7  38.8  88.1   0.0  41.6
 88.4  12.1  41.3  88.4  99.6  41.2  88.8  82.9  40.5
 88.9   6.2  41.5  90.6   7.0  41.5  90.7  49.6  38.9
 91.5  55.4  39.0  92.9  46.8  39.1  93.4  70.9  39.7
 94.8  71.5  39.7  96.2  84.3  40.3  98.2  58.2  39.5
];
pos=[data(:,1:2);data(:,4:5);data(:,7:8)];
z=[data(:,3);data(:,6);data(:,9)];
%}