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;