1

Two datasets are used:

  1. spatial data in 3 columns (x, y, data)

  2. grids data in 2 columns (x, y)

automap package autoKrige does the calculation of kriging and can be plotted with no x and y tick marks and labels:

plot(kriging_result)
automapPlot(kriging_result$krige_output, "var1.pred", sp.layout = list("sp.points", shimadata), main="OK without grids", xlab="x", ylab="y")

And when I use ggplot2 package, it shows error, however it calculates the kriging:

mydata<-read.table("D:/.../mydata.txt",header=T,sep=",")
#Renaming desired columns:
x<-mydata[,1]
y<-mydata[,2]
waterelev<-mydata[,3]
library(gstat)
coordinates(mydata)=~x+y
library(ggplot2)
theme_set(theme_bw())
library(scales)
library(automap)
grids<-read.table("D:/.../grids.txt",header=T,sep=",")
gridded(grids)=~x+y
kriging_result = autoKrige(log(waterelev)~1, mydata)
#This line turns the log(data) back to the original data:
kriging_result$krige_output$var1.pred<-exp(kriging_result$krige_output$var1.pred)
library(reshape2)
ggplot_data = as.data.frame(kriging_result$krige_output)
ggplot(ggplot_data, aes(x = x, y = y, fill = var1.pred)) + 
   geom_raster() + coord_fixed() + 
   scale_fill_gradient(low = 'white', high = muted('blue'))

The error:

Error: Aesthetics must either be length one, or the same length as the dataProblems:x, y

2 Answers2

3

The automapPlot function is a wrapper around spplot, so any fix that works for spplot would also be applicable to automapPlot. You can start with the spplot documentation for a place to start.

However, I normally use the ggplot2 package for any plotting work, including spatial data. The following is an example which reproduces the result of automapPlot:

library(ggplot2)
theme_set(theme_bw())
library(scales)
library(automap)
data(meuse)
coordinates(meuse) =~ x+y
data(meuse.grid)
gridded(meuse.grid) =~ x+y
kriging_result = autoKrige(zinc~1, meuse, meuse.grid)

# Cast the Spatial object to a data.frame
library(reshape2)
ggplot_data = as.data.frame(kriging_result$krige_output)

ggplot(ggplot_data, aes(x = x, y = y, fill = var1.pred)) + 
    geom_raster() + coord_fixed() + 
    scale_fill_gradient(low = 'white', high = muted('blue'))

enter image description here

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
  • Thanks a lot for your answer. But, when I apply these lines for my data, it says: "Error: Aesthetics must either be length one, or the same length as the dataProblems:x, y" What should I do to fix that? – user3918211 Aug 15 '14 at 15:28
  • Dear Mr. @PaulHiemstra, I found out that geom_raster() is probably causing this error: "Error: Aesthetics must either be length one, or the same length as the dataProblems:x, y" What should I do to fix that? – user3918211 Aug 15 '14 at 15:43
  • Before I can provide any feedback, you need to add code+data which reproduces your problem to your question above. – Paul Hiemstra Aug 15 '14 at 19:52
  • Thank you Mr. @PaulHiemstra. I have added the code, but since I didn't know how to attach a file, I just explained the data format I use. – user3918211 Aug 15 '14 at 21:52
  • Dear Mr. @PaulHiemstra I just learned how to share a file. Here is the link: https://onedrive.live.com/redir?resid=66BD1F9292741792!111&authkey=!AAdRpPlhf6Abf_s&ithint=file%2crar ...It should be noted that the address of the files in code are not specified. – user3918211 Aug 17 '14 at 04:31
1

I used the dataset you provided, and, I post the result like an image:

ggplot result

I just run your code until ggplot_data is produced casting from an spPointsDataFrame (gridded), to a data.frame (with reshape2). Then I built the ggplot object step by step but I found no errors:

g=ggplot(data=ggplot_data,aes(x=x1,y=x2,fill=var1.pred))
g=g+geom_raster()
g=g+coord_fixed()
g=g+scale_fill_gradient(low = 'white', high = muted('blue'))
print(g)

Is this what you expected?

Fabio
  • 518
  • 1
  • 4
  • 10
  • Thank you so much @Fabio. It worked, but I have 5 problems: 1. Since I produced 2 kriging_results (with and without using grids), I don't know which one is plotted by this code (I want to map both of them) 2. x and y axis labels can't be changed with xlab and ylab to x,y. 3. I can't interpolate with universal, simple & block kriging, is it possible with the data I have? And how? The 4th problem is in the next comment!! – user3918211 Aug 17 '14 at 16:03
  • 4. I need to plot the standard error of both kriging_results, but I do not know how, and how to plot them with x and y labels and tick marks 5. I don't know how to evaluate the IDW interpolation and how to plot the SE or ... of it, or how to compare IDW and kriging. Please forgive me, I am new to R, and the code I attached is the first code I have written. Could you please add the needed code to mine and send me the link, or write it here? Thank you very much for all your help. I am sorry again. – user3918211 Aug 17 '14 at 16:06
  • Oh, I just realized you meant the code above, not the one I shared its link. So, problems 1 and 2 are solved for me. Sorry!! – user3918211 Aug 17 '14 at 16:35
  • And, one last thing I realized is that the gradient it shows, doesn't have the contrast needed to get different shades. When we plot(kriging_result), it shows the levels better. Is there a way to add more colors in the gradient?? – user3918211 Aug 17 '14 at 17:02
  • sure it is: use the scale_colour_gradientn(colours = rainbow(7), breaks = c(700,800,900,1000), labels = format(breaks)) – Fabio Aug 17 '14 at 20:10
  • Just to keep track of the issues you want to solve could you please post both a new question, one at the time and in order of importance, with a link here? Thank you. (By the way, try to sign the answers you find useful...othewise this question starts becoming a mess) – Fabio Aug 17 '14 at 20:13
  • Thank you for your time, @Fabio. I will add questions somewhere else. And I signed your answer as the accepted one. :) – user3918211 Aug 18 '14 at 06:48
  • Just add them to StackOverflow, but pèlease link them here too so I can check them quickly – Fabio Aug 18 '14 at 07:44
  • I have added it @Fabio. Thanks. link: http://stackoverflow.com/questions/25363092/how-to-do-cross-validation-for-block-kriging – user3918211 Aug 18 '14 at 12:17