0

我有一个数据框，其中包含 2 年 4 个站的点日降水量。我想插值到 50m 分辨率并将它们写入 2 个光栅图像。我使用以下代码来实现这一点......

library(ggplot2)
library(gstat)
library(maptools)
library(raster)
library(rgdal)
xcord<-c(100,200,300,400)
ycord<-c(100,200,300,400)
value1<-c(1,2,3,1)
value2<-c(2,5,7,3)
datas<-data.frame(xcord,ycord,value1,value2)
coordinates(datas) = ~xcord+ycord
mesh <- expand.grid(x=seq(0,500,50),y=seq(0,500,50))
coordinates(mesh)=~x+y
gridded(mesh) <- TRUE
oneidw = idw(value1~1,datas,mesh)
spplot(oneidw["var1.pred"], main = " inverse distance weighted    interpolations")

有效。但我想应用一个循环来为另一个变量value2（等等......）而不是手动执行它。我用这个

    sym<-paste("value", 1:2,sep="")
    variable=as.vector(print(sym,quote=FALSE))
    for (i in 3:ncol(datas)){
     one<-idw((print(variable[i],quote=FALSE))~1,datas,mesh)
    }

但我得到了太多空间维度的错误............有人可以帮我解决这个问题......

4

1 回答 1

0

我对spplot不太熟悉，但这对我使用ggplot很有用。

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
library(maps)
library(dplyr)
library(rgdal)

xcord<-c(100,200,300,400)
ycord<-c(100,200,300,400)
value1<-c(1,2,3,1)
value2<-c(2,5,7,3)
datas<-data.frame(xcord,ycord,value1,value2)
new_datas <- select(datas, xcord, ycord)
parse_by <- colnames(datas)[3:4] #change according to designated value columns 

for ( i in parse_by ) { 
  variable <- datas[i]
  new_datas2 <- cbind(new_datas, variable) #combine single variable col w/ coordinates
  colnames(new_datas2)[3] = "variable" #rename so that you can call to in idw formula
  coordinates(new_datas2) = ~xcord+ycord
  mesh <- expand.grid(x=seq(0,500,50),y=seq(0,500,50))
  coordinates(mesh)=~x+y
  gridded(mesh) = TRUE
  plot(mesh) #plot background so ggplot can use later
  points(new_datas2) #points for ggplot to use later
  one<-idw(formula = variable~1, locations = new_datas2, newdata = mesh) #idw formula
  one.output <- as.data.frame(one)
  names(one.output)[1:3] <- c("xcord", "ycord", "var1.pred") #rename for ggplot geom_tile
  ggplot() + geom_tile(data = one.output, alpha = 1, aes(x = xcord, y = ycord, fill = var1.pred)) + 
    geom_point(data = new_datas, aes(x = xcord, y = ycord)) + 
    labs(fill = "inverse distance weighted interpolations") 
  ggsave(paste(i,".png",sep="")) #save as .png according to variable name
}
于 2016-02-16T20:56:31.913 回答