Kriging: Contour Estimation from Sparse Data
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: KRIGING */
/* TITLE: Kriging: Contour Estimation from Sparse Data */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: MATRIX */
/* PROCS: IML PLOT SORT PRINT */
/* DATA: */
/* */
/* SUPPORT: LWB UPDATE: SEP2013 */
/* REF: */
/* MISC: TRANSLATED FROM MATRIX TO IML USING MATIML */
/* */
/****************************************************************/
* INPUT RAW DATA ALL AS ONE OBSERVATION ----TEN POINTS;
DATA ONE;
INPUT Y1-Y10 Z1O1-Z1O10 Z2O1-Z2O10;
CARDS;
2 1.5 2.1 2 1.9 1.7 2.2 2.6 2.9 1.7 .3 .5 .3 .5 .7 .9 .6 .4 .2 .9 .2 .3
.6 .6 1.0 1.1 1.2 1.4 1.5 1.8
;
RUN;
* GENERATE A COARSE GRID OF POINTS TO COMPUTE ESTIMATES AT;
DATA TWO;
KEEP Z1 Z2 C1-C10 ;
ARRAY C C1-C10 ;
ARRAY Z1O Z1O1-Z1O10;
ARRAY Z2O Z2O1-Z2O10;
SET ONE;
DO Z1=0 TO 1 BY .1;
DO Z2=0 TO 2 BY .1;
DO OVER C;
*COMPUTE DISTANCES FROM RAW DATA TO GRID VALUES;
C=SQRT((Z1O-Z1)**2+(Z2O-Z2)**2);
*USING THE BELOW ASSUMED COVARIANCE FUNCTION--NOTE C(0)=1;
*ALSO THE MEAN IS ASSUMED TO BE ZERO ;
IF C > 1 THEN C=0;
ELSE C=1-1.5*C+.5*C**3;
END;
OUTPUT;
END;
END;
RUN;
*COMPUTE THE KRIG ESTIMATES WITHIN PROC IML;
PROC IML;
USE TWO ;
READ ALL INTO Z ;
USE ONE ;
READ ALL INTO ZO ;
*COMPUTE THE VARIANCE/COVARIANCE MATRIX OF THE RAW DATA POINTS;
Z1= ZO[1,11:20];
Z2= ZO[1,21:30];
*ALL POSSIBLE DIFFERENCES OF COORDINATES;
A= J(10,1,1);
Z1=Z1@A-(A@Z1)`;
Z2=Z2@A-(A@Z2)`;
*COMPUTE THE COVARIANCES;
CO= SQRT(Z1##2+Z2##2);
CO=CO><1;
CO=1-(1.5*CO)+.5*(CO##3);
*LOAD THE MATRICES Y AND C AND INVERT THE RAW DATA VAR/COV MATRIX;
Y= ZO[1,1:10];
C=( Z[,1:10])`;
Z= Z[,11:12];
COINV= INV(CO);
*COMPUTE THE ESTIMATES OF THE UNKNOWN VALUES M=Y*COINV*C;
B=Y*COINV;
M=B*C;
*GET DIAGONALS OF VARIANCE MATRIX;
S=C#(COINV*C);
S= S[+,];
S2=1-S;
S3=CHOOSE(S2<1e-12 & S2>=-1e-12,0,S2);
S= SQRT(S3);
*RECALL THAT C(0)=1;
*CONCATENATE THE GRID POINTS AND THE ESTIMATES FOR AN OUTPUT DATASET;
Z=Z||(M`)||(S`);
CNAME={'Z1' 'Z2' 'MEAN' 'STD'};
_TMP_ROW = 'ROW1 ' : compress('ROW'+char(nrow(Z)));
CREATE FOUR ( RENAME=(_TMP_ROW=ROW )) FROM Z
[ROWNAME=_TMP_ROW COLNAME=CNAME ];
APPEND FROM Z [ROWNAME=_TMP_ROW];
QUIT;
proc template;
define statgraph ContourPlotParm;
dynamic _X _Y _Z _TITLE;
begingraph;
entrytitle _TITLE;
layout overlay;
contourplotparm x=_X y=_Y z=_Z /
contourtype=fill nhint=12 colormodel=twocolorramp name="Contour";
continuouslegend "Contour" / title=_Z;
endlayout;
endgraph;
end;
run;
TITLE3 'COARSE KRIG ESTIMATES';
proc sgrender data=Four template=ContourPlotParm;
dynamic _TITLE="Mean Field"
_X="z1" _Y="z2" _Z="Mean";
run;
proc sgrender data=Four template=ContourPlotParm;
dynamic _TITLE="Standard Deviation Field"
_X="z1" _Y="z2" _Z="Std";
run;
* DO QUICK AND DIRTY LINEAR INTERPOLATION;
DATA Five;
SET FOUR;
BY Z1 Z2;
KEEP Z1 Z2 MEAN STD;
RETAIN ZL ML SL;
IF FIRST.Z1 THEN DO;
ZL=Z2;
SL=STD;
ML=MEAN;
OUTPUT;
END;
ELSE DO;
INC=(Z2-ZL)/4;
SLOPEM=(MEAN-ML)/4;
SLOPES=(STD-SL)/4;
Z2=ZL;
STD=SL;
MEAN=ML;
DO I=1 TO 4;
Z2=Z2+INC;
MEAN=MEAN+SLOPEM;
STD=STD+SLOPES;
OUTPUT;
END;
ZL=Z2;
SL=STD;
ML=MEAN;
END;
RUN;
PROC SORT data=Five;
BY Z2 Z1;
RUN;
* DO MORE QUICK AND DIRTY INTERPOLATION;
DATA Five;
SET Five;
KEEP Z1 Z2 MEAN STD;
RETAIN ZL ML SL;
BY Z2 Z1;
IF FIRST.Z2 THEN DO;
ZL=Z1;
SL=STD;
ML=MEAN;
OUTPUT;
END;
ELSE DO;
INC=(Z1-ZL)/4;
SLOPEM=(MEAN-ML)/4;
SLOPES=(STD-SL)/4;
Z1=ZL;
STD=SL;
MEAN=ML;
DO I=1 TO 4;
Z1=Z1+INC;
STD=STD+SLOPES;
MEAN=MEAN+SLOPEM;
OUTPUT;
END;
ZL=Z1;
SL=STD;
ML=MEAN;
END;
RUN;
* PLOT THE INTERPOLATED VALUES ;
TITLE3 'KRIG ESTIMATES WITH ASSUMED COVARIANCE FUNCTION';
TITLE4 'C(R)=1-1.5*R+.5*R**3';
proc sgrender data=Five template=ContourPlotParm;
dynamic _TITLE="Mean Field"
_X="z1" _Y="z2" _Z="Mean";
run;
proc sgrender data=Five template=ContourPlotParm;
dynamic _TITLE="Standard Deviation Field"
_X="z1" _Y="z2" _Z="Std";
run;
title;title2;title3;title4;