Skip to content
Prev 23237 / 29559 Next

How to plot multiple semi-variogram from a single dataset efficiently in R?

Orpheus,
here's one option to get what you want.
The problem with fitting that many variograms is that if you want to 
specify each variogram individually (i.e. nugget, range, model etc.), 
there is no way to automate this with loops or thelike.
You can, however, use library(automap) to automatically fit a variogram 
to your data. Then a possible solution could look like this:


### in order to use latticeCombineGrid() you need to use ##################
### library(devtools) to install library(Rsenal) ##########################
# install.packages(devtools)
# library(devtools)
# install_github("environmentalinformatics-marburg/Rsenal")
###########################################################################

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

seoul311 <- read.csv("Downloads/seoul1to7.csv")
seoul311 <- na.omit(seoul311)

### first we split seoul311 by time into a list rather than subsetting 
manually
seoul311_splt <- split(seoul311, seoul311$time)

### now we loop (using lapply()) over each seoul311_splt entry and 
calculate
### variogram using autofitVariogram and return the variogram plot
vars <- lapply(seq(seoul311_splt), function(i) {

   dat <- seoul311_splt[[i]]
   coordinates(dat) <- ~LON+LAT
   proj4string(dat) <- "+proj=longlat +datum=WGS84"
   dat <- spTransform(dat, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

   variogram <- autofitVariogram(PM10 ~ 1, dat)

   plt <- 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, 45000), ylim = c(0, 1000))

   return(plt)
})

### now we actually have 23 * 7 variogram plots which we will combine
### into 23 hourly plots using latticeCombineGrid()
hrs <- substr(names(seoul311_splt), 9, 10)

plts_hrs <- lapply(seq(unique(hrs)), function(j) {

   indx <- hrs %in% unique(hrs)[j]
   hr_plt <- vars[indx]

   return(latticeCombineGrid(hr_plt, layout = c(3, 3)))

})

### plot for hour 1
plts_hrs[[1]]

### plot for hour 23
plts_hrs[[23]]

### ...


HTH

Best,
Tim
On 17.08.2015 09:43, Uzzal wrote: