## 8. Kriging estimation

 8.1 Create 2D Mask (download) Create a 2D grid mask from sample data and a maximum distance. Variables to customize: varcolXY<- c(1,2) # columns of the coordinates in the GTD data frame `distanciaXY <- c(200,200) # maximum distance os a samples in X and Y directions to select the cells maskshape <- "rectangle"  # Shape for selection ellipse or rectangle xmincenter <- 25 # X coordinate of the left bottom cell centerymincenter <- 25 # Y coordinate of the left bottom cell centerncellsx <- 55 # number of cells in X directionncellsy <- 65 # number of cells in Y directioncellsizex <- 100 # cell size in X directioncellsizey <- 100 # cell size in Y directioncellvalue <- 1 #value to assign to the selected cells nodata <- -9  #value to assign to the selected cells outside maskcolour <- c(rgb(1,0,0,0.25),rgb(.8,.8,.8,0.5)) #colours within mask and outsidemaintitle <- "Área de estudo solos 50m x 50m" # title of the graphicfileMaskGrid <- "C:/temp/Mask.OUT" # pathname and filename of a file to store the mask grid numbers. Values are writted by Y and then by X` 8.2 Estimation using the inverse of power distance (download) Estimation of a grid mesh by using the estimator inverse of power distance. Variables to customize: `varcolXY <- c(1,2) # columns of the coordinates in the GTD data framevarlistindex <- 5 # column of the variable in the GTD data frame to estimatexmincenter <- 25 # X coordinate of the left bottom cell centerymincenter <- 25 # Y coordinate of the left bottom cell centerncellsx <- 55 # number of cells in X directionncellsy <- 65 # number of cells in Y directioncellsizex <- 100 # cell size in X directioncellsizey <- 100 # cell size in Y directionnodata <- -9 #value to assign to non estimated locationspower <- 2 # power value to assign to the distancemaxsamples <- 8 # maximum number of samples to usemaintitle <- paste("Estimação pelo inverso da potência da distância: ",power) # title of the imagefileMaskGrid <- "C:/temp/Mask.OUT" # file with mask valuesfileOutGrid <- "C:/temp/Cd_IPD.out" # output file with the estimated results. The structure is the same of the mask file` 8.3 Cross validation test for inverse power distance estimator  (download) Cross validation test for estimation using the inverse of power distance, powers can range between a minimim and a maximum value, error and squared error are computed for each power and represented graphically, and the best power can be evaluated. Variables to customize: `varcolXY <- c(1,2) # columns of the coordinates in the GTD data framevarlistindex <- 5 # index of the GTD frame quantitative variable to test for estimationminpower <- 0 # minimum power applied to distance to evaluatemaxpower <- 4 # minimum power applied to distance to evaluateintpower <- 0.1 # increment of the power value to evaluatemaxsamples <- 20 # maximum number of samples to use in each estimationescolhabi <- 12 # index of the estimation to plot a scattergram between true values and estimated values (graphic on the right)titleGRAPH <- "Avaliação do estimador potência da distância" # title of the graphic` 8.4 Ordinary and Simple Kriging demo of results (download) Demo of ordinary and simple kriging, inititate a simple configuration of samples and perform kriging of an unknown location to evaluate kriging weights, and estimated value. Variables to customize: `ktype <- "KS" # KN (ordinary kriging) or KS (simples kriging)titleXY <- c("X","Y") # title of the axis X and Yxloc <- c(-3,-2,-1,0,0,0,1.0,1.0) # X coordinates of the samples yloc <- c(0,0,0,1,2,-2,0.3,-0.3) # Y coordinates of the samplesvalor <- c(12.5,10.2,12.4,9.5,12,11,12,13) # list of hypothetical true valuesxnode <- 0 # X coordinate of the location to estimateynode <- 0 # Y coordinate of the location to estimatemedialocal <- mean(valor) # local mean for simple kriging, can be a value or the mean of the true data##  Variogram models. Models are inputed by proportion of variance (sum of Nugget effect + Sills = 1)#nugget_effect <- 1 # nugget effect valuenstructures <- 0 # number of structures (0, 1 or 2)model_type <- c("Sph","Sph") # model type for each structure, can be Sph, Exp, Gauss, or Powersill <- c(0.1,1.0) # sill of each modelrange <- c(2,6.0) # range of each modelanisotropy <- c(1,1) # anisotropy factor, 1 or highermaindir <- c(0.0,0.0) # in anisotropic, write the azimuth angle of the main direction clockwise, N direction is zero reference` 8.5 kriging cross validation test (download) Cross validation test by using Simple Kriging or Ordinary Kriging. Values are estimated in the samples locations, and errors (simple and squared) are computed. Two graphics show estimation errors, first a scatterplot of true values vs estimated values and then underestimated and overestimated locations are displayed. Variables to customize: `ktype <- "KN" # KN (ordinary kriging) or KS (simples kriging)varcolXY <- c(1,2) # columns of the X and Y coordinated in GTD data framevarlistindex <- 7 # column of the variable to estimate in GTD data framemaxpercent <- 50 # percentage of mean for classification of high deviation (values in black)nclasses <- 10 # number of classes of the histogramnodata <- -9 # value for unestimated locationsmaxsamples <- 10 # maximum number of closest samples to use in the estimationmedialocal <- mean(gtd[,varlistindex]) # global mean for simples kriging estimation, can be the average of the variable or other valuefileOutTVC="C:/temp/Cd_TVC.OUT" # output file results##  Models are inputed by proportion of variance (sum of Nugget effect + Sills = 1)#nugget_effect <- 0.0 # nugget effect valuenstructures <- 1 # number of structures (0, 1 or 2)model_type <- c("Sph","Sph") # model type for each structure, can be Sph, Exp, Gauss, or Powersill <- c(1.0,1.0) # sill of each modelrange <- c(1600,1) # range of each modelanisotropy <- c(2,1) # anisotropy factor, 1 or highermaindir <- c(60.0,0.0) # in anisotropic, write the azimuth angle of the main direction clockwise, N direction is zero reference` 8.6 Grid estimation by using simple or ordinary kriging (download) Estimation of a mesh of locations by using Simple Kriging or Ordinary Kriging. Three graphics show estimation results, the first is the estimated map of values, the second is the kriging variance map and the third displays two histograms, one of the sample data and another of the estimated values. Variables to customize: `ktype <- "KN"  # KN (ordinary kriging) or KS (simples kriging)varcolXY <- c(1,2) # columns of the X and Y coordinated in GTD data framevarlistindex <- 5 # column of the variable to estimate in GTD data framenclasses <- 20 # number of classes of the histogramlinedensity <- TRUE # add line of density in the histogramsxmincenter <- 25 # X coordinate of the left bottom cell centerymincenter <- 25 # Y coordinate of the left bottom cell centerncellsx <- 55 # number of cells in X directionncellsy <- 65 # number of cells in Y directioncellsizex <- 100 # cell size in X directioncellsizey <- 100 # cell size in Y directionthreshold <- -9  # Cut off value to binarize the variable and perform kriging of an indicator variablefileMaskGrid <- "C:/temp/Mask.out" # file with mask valuesfileOutGrid <- "C:/temp/Cd_KN.OUT" # output file with the estimated results. The structure is the same of the mask file^nodata <- -9 # value for unestimated locationsmaxsamples <- 8 # maximum number of closest samples to use in the estimationmedialocal <- mean(gtd[,varlistindex]) # global mean for simples kriging estimation, can be the average of the variable or other valueviewLabels <- FALSE # View true values (TRUE) or not (FALSE)##  Models are inputed by proportion of variance (sum of Nugget effect + Sills = 1)#nugget_effect <- 0.0 # nugget effect valuenstructures <- 1 # number of structures (0, 1 or 2)model_type <- c("Sph","Sph") # model type for each structure, can be Sph, Exp, Gauss, or Powersill <- c(1.0,1.0) # sill of each modelrange <- c(1600,1) # range of each modelanisotropy <- c(2,1) # anisotropy factor, 1 or highermaindir <- c(60.0,0.0) # in anisotropic, write the azimuth angle of the main direction clockwise, N direction is zero reference`  8.7 View 2D grids and sample data (download) Visualization of previously estimated images by inverse of power distance or kriging. Variables to customize: `varcolXY <- c(1,2)  # columns of the X and Y coordinated in GTD data frametitleXY <- c("X","Y") # title of X and Y axisvarlistindex <- 7 # column of the variable to estimate in GTD data framexmincenter <- 25 # X coordinate of the left bottom cell centerymincenter <- 25 # Y coordinate of the left bottom cell centerncellsx <- 55 # number of cells in X directionncellsy <- 65 # number of cells in Y directioncellsizex <- 100 # cell size in X directioncellsizey <- 100 # cell size in Y directionnodata <- -9 # value for unestimated locationsfileNameGrid <- "C:/temp/Cd_KN.OUT" # name of the file to import and viewmaintitle <- "Máscara de delimitação de área de estudo (50x50 m)" # Title of the 2D graphical display` 