1

I would like to save the information generated from an ordinary kriging (processed with the function autoKrige from the automap R package), into a file extension readable in ArcGis for desktop 9.3.1 (more specifically, an .asc file - aka ESRI ASCII raster format).

I tried to search for functions inside the raster (raster) and rgdal (writeOGR) packages but without success.

require(gstat)
require(automap)
require(sp)

#Reading data to provide a reproducible example.
data(meuse)
coordinates(meuse) =~ x+y
data(meuse.grid)
gridded(meuse.grid) =~ x+y

#Perform ordinary kriging and store results on object of type "autoKrige" "list" 
kriging_result = autoKrige(zinc~1, meuse, meuse.grid)

plot(kriging_result) #I clipped one panel from the entire output plot.

enter image description here

#Extract coordinates and variable of interest to object of type "data.frame"
kriging_variable = data.frame(kriging_result[[1]])

head(kriging_variable)

#Example of observations for the variable of interest "var1.pred"
x      y       var1.pred  var1.var   var1.stdev
181180 333740  777.5909   96809.44   311.1421
181140 333700  839.7388   78629.08   280.4088
181180 333700  773.9795   84108.76   290.0151
181220 333700  710.5693   90272.17   300.4533
181100 333660  911.2036   61045.30   247.0735
181140 333660  836.4262   65746.42   256.4107

Does anyone have a clue?

3

You can use "writeGDAL" in the rgdal package. Use the "asc" extension to specify an ESRI ASCII raster format. You will need to pull the associated sp object and pass it to writeGDAL. This is vaguely indicated in the help: "This function returns an autoKrige object containing the results of the interpolation (prediction, variance and standard deviation), the sample variogram and the variogram model, The attribute names are krige_output, exp_var, var_model and sserr respectively".

    require(rgdal)

    # Pull sp SpatialPixelsDataFrame 
        kriging.pred <- kriging_result$krige_output

    # Write Kriging estimate and variance
        writeGDAL(kriging.pred["var1.pred"], "KrigingPred.asc")
        writeGDAL(kriging.pred["var1.var"], "KrigingVar.asc")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.