Saturday, 3 November 2018

gdal - Exporting kriging prediction map in .asc format?



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?



Answer



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")

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...