0

I am using the automap library where I want to use autokrige to interpolate rainfall data. The resulting interpolated map has different dimensions (different number of rows and colums) than my prediction points map. How is this possible? Please see my code below:

library(automap)
library(raster)
library(rgdal)

r <- read.table("c:\\Active\\2013-001_Qfever\\model\\data\\output_kriging\\raindata_old1.txt", header = TRUE, skip = 2, sep="")
coordinates(r) = ~x+y

predmap <- readGDAL(fname="c:\\Active\\2013-001_Qfever\\model\\NLclone_scalar.map")
kriging_result <- autoKrige(rain~1, r, predmap)

plot(kriging_result)

xyz <- as.data.frame(kriging_result$krige_output)                                                 
head(xyz)
p = xyz[,c('x','y','var1.pred')]

output <- rasterFromXYZ(p, res=c(250,250), crs="+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs", digits=5)

writeRaster(output, "c:\\Active\\2013-001_Qfever\\model\\data\\output_kriging\\test.asc", format="ascii")
plot(output)
Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
user2357183
  • 43
  • 1
  • 3
  • Right now your example is not reproducible as we do not have the data. You can either try and reproduce your situation with datasets available in automap (e.g. the meuse dataset), or you can provide us with the files, or samples thereof. See also http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example. – Paul Hiemstra May 07 '13 at 14:52

1 Answers1

0

Replace the following lines in your code, it will work-

p = xyz[,c('x1','x2','var1.pred')]
ToNoY
  • 1,358
  • 2
  • 22
  • 43