如何在 R 中有效地从单个数据集中绘制多个半变异函数？答案

【问题标题】：How to plot multiple semi-variogram from a single dataset efficiently in R?如何在 R 中有效地从单个数据集中绘制多个半变异函数？
【发布时间】：2015-10-24 02:15:59
【问题描述】：

我有一个名为 seoul1to7 的数据框，其中包含 2012 年 3 月 1 日至 3 月 7 日的每小时 PM10 浓度数据。please,download.在此数据集中，时间采用 yyyymmddhr 格式。例如，2012030101 表示 2012 年 3 月 1 日凌晨 1 点

数据看起来像：

      ID       time PM10      LAT     LON
1 111121 2012030101   42 37.56464 126.976
2 111121 2012030102   36 37.56464 126.976
3 111121 2012030103   46 37.56464 126.976
4 111121 2012030104   40 37.56464 126.976
.
.

我的最终目标是绘制每小时的半变异函数。例如，对于 2012 年 3 月 1 日上午 1 点（2012030101），有 107 个 PM10 数据。我想绘制 2012030101 到 2012030723 的半变异函数（总共 7*24 半变异函数）。我用 R 写了一些代码：

seoul1to7<-read.csv("seoul1to7.csv", row.names=1)
rownames(seoul1to7)<-NULL

seoul311<-subset(seoul1to7, time==2012030101)
seoul312<-subset(seoul1to7, time==2012030102)
.
.
.
seoul3723<-subset(seoul1to7,time==2012030724)

起初，我尝试通过subset() 函数制作我想要的 (7*24) 数据帧然后我想为每个数据帧绘制半变异函数。例如，我通过以下代码为seoul311（对于 2012030101）绘制了半变异函数：

library(sp)
library(gstat)
library(rgdal)
seoul311<-read.csv("seoul311.csv",row.name=1)
seoul311<-na.omit(seoul311)

coordinates(seoul311)=~LON+LAT
proj4string(seoul311) =  "+proj=longlat +datum=WGS84" 
seoul311<-spTransform(seoul311, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

#plot Omnidirectional Variogram
seoul311.var<-variogram(PM10~1,data=seoul311,cutoff=66000, width=6000)
seoul311.var
plot(seoul311.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 311")

#Model fit
model.311<- fit.variogram(seoul311.var,vgm(psill=250,model="Gau",range=40000,nugget=100),
                          fit.method = 2)
plot(seoul311.var,model=model.311, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")


#Directional Variogram
seoul311.var1<-variogram(PM10~1,data=seoul311,width=6000,cutoff=66000,
                         alpha=seq(0,135,45),tol.hor=15)
seoul311.var1
plot(seoul311.var1,model=model.311, cex=1.1,pch=16,col=1,
     main="ANisotropic Variogram for PM10")    


#anisotropy corrected variograms
model.3112.anis<- fit.variogram(seoul311.var1,vgm(250,"Gau",40000,100,anis=c(45,0.80)),
                                fit.method = 2)

#Final isotropic variogram for kriging
plot(seoul311.var,model=model.3112.anis, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Final Isotropic Variogram")

但我知道我的代码效率很低！我正在写 (7*24) 次 subset(seoul1to7, time==2012030101) 这段代码。然后再次（7 * 24）次绘制半变异函数的代码！我认为这是非常不合适的方式。

那么，我怎样才能从我的数据集seoul1to7 中非常有效地绘制这些 (7*24) 半变异函数（通过使用循环或任何其他函数）？如果您需要任何进一步的信息，请告诉我。

【问题讨论】：

您是否查看了variogramST 的时空变异函数？
其实我这里只是做空间分析。因为时间是固定的。所以我没有检查variogramST 函数。相反，我尝试通过循环或通过应用家庭来发现我可以做到或不做到。但没有成功！
我应该编辑我的问题并使其简短吗？我只是对这个问题的工作感到震惊。
或者我应该在其他论坛发布这个问题吗？请告诉我？
a<-lapply(unique(seoul1to7$time), function(x) subset(seoul1to7, time==x))。通过此代码，我可以将所有数据框作为列表获取。然后我可以通过这种方式调用a[1], a[2]... 来查看每个数据帧。但是如果我想为每个数据帧（161 个变异函数）绘制半变异函数，我应该写半变异函数代码（我上面显示的）161 次吗？有没有其他办法？

标签： r loops gstat automap

【解决方案1】：

library(sp)
library(gstat)
library(rgdal)
library(automap)

seoul1to7<-read.csv("seoul1to7.csv", row.names=1)
seoul1to7 <- na.omit(seoul1to7)
seoul1to7_split<-split(seoul1to7,seoul1to7$time) #a list contain 161 dataframe
seq(seoul1to7_split)

### now we loop (using lapply()) over each seoul1to7_split entry and calculate 
### variogram using autofitVariogram and return the variogram plot

vars<-lapply(seq(seoul1to7_split), function(i)
{
  dat<-seoul1to7_split[[i]]  # for list element [[]]
  coordinates(dat)<-~LON+LAT
  proj4string(dat) <- "+proj=longlat +datum=WGS84" 
  dat <- spTransform(dat, CRS("+proj=utm +north +zone=52 +datum=WGS84"))
  variogram<-autofitVariogram(log(PM10)~1,dat)
  plot<- plot(variogram,plotit=FALSE, asp=1)

  ### in case you do not want to fix xlim and ylim to be identical 
  ### for each plot just comment out the following line or change
  ### values as you see fit
  #plt <- update(plt, xlim = c(-1000, 35000), ylim = c(0, 1000))

  return(plot)
})

### now we actually have 23 * 7 variogram plots which we will combine 
### into 23 hourly plots using latticeCombineGrid()
library(raster)
library(devtools)
#install.packages("Rcpp")
install_github("environmentalinformatics-marburg/Rsenal")
library(Rsenal)

names(seoul1to7_split)
hours<-substr(names(seoul1to7_split),9,10) #substructuring 9th to 10th digit of every element
hours
unique(hours)
class(unique(hours))   #character
seq(unique(hours))

var7_per_plot<-lapply(seq(unique(hours)), function(j)
  {
  index<- hours %in% unique(hours)[j]  
  plot.hours <- vars[index]

  return(latticeCombineGrid(plot.hours, layout =c(3,3)))
})
var7_per_plot[[1]]
var7_per_plot[[2]]
.
.
.
var7_per_plot[[23]]

特别感谢 Tim Appelhans 教会了我。

【讨论】：