Wednesday 21 September 2011

Variogram fit with RPanel

During the UseR 2011 conference I saw lots of examples of the use of RPanel to create a GUI in R. Yesterday, because I was a bit bored of the work I was doing I started thinking about this and I decided to try this package.

My objective was to create a new panel with all the main setting for fitting a variogram model to an omnidirectional variogram with gstat.

This is the script:

library(rpanel)
library(gstat)

data<-read.table("data.txt",sep="",header=T)
coordinates(data)=~Lat+Lon

##Variogram Fitting
variogram.plot <- function(panel) {
with(panel, {
variogram<-variogram(Oxigen~1,data,cutoff=Cutoff)
vgm.var<-vgm(psill=Sill,model=Model,range=Range,nugget=Nugget)
#fit<-fit.variogram(variogram,vgm.var)
plot(variogram$dist,variogram$gamma,xlab="Distance",ylab="Semivariance")
lines(variogramLine(vgm.var,maxdist=Range))
})
panel
}

var.panel <- rp.control("Variogram",Sill=20,Range=250,Nugget=0,Model="Mat",Cutoff=250)
rp.listbox(var.panel,Model,c("Mat","Sph","Gau"))
rp.slider(var.panel, Cutoff, 0,500,showvalue=T)
rp.slider(var.panel, Sill, 0,500,showvalue=T)
rp.slider(var.panel, Range, 0,1000,showvalue=T)
rp.slider(var.panel, Nugget, 0,15,showvalue=T)
rp.button(var.panel, title = "Fit", action = variogram.plot)


At this address you can find a zip file with a sample dataset that you can use to try this script, however if you know a bit of gstat you can start customizing it straigth away:


This is the screenshot from my R Console:




I also tried to embed the variogram plot into the panel in order to execute the script in batch mode, this is the result:


if(print(require(tcltk))==FALSE){install.packages("tcltk",repos="http://cran.r-project.org")}
if(print(require(tcltk))==TRUE){require(tcltk)}

if(print(require(rpanel))==FALSE){install.packages("rpanel",repos="http://cran.r-project.org")}
if(print(require(rpanel))==TRUE){require(rpanel)}

if(print(require(gstat))==FALSE){install.packages("gstat",repos="http://cran.r-project.org")}
if(print(require(gstat))==TRUE){require(gstat)}



data<-read.table(tclvalue(tkgetOpenFile()),sep="",header=T)
coordinates(data)=~Lat+Lon

grid <- spsample(data, cellsize = 10, type = "regular")
gridded(grid) <- TRUE


##Variogram Fitting
variogram.plot <- function(panel) {
with(panel, {
variogram<-variogram(Oxigen~1,data,cutoff=Cutoff)
vgm.var<-vgm(psill=Sill,model=Model,range=Range,nugget=Nugget)
#fit<-fit.variogram(variogram,vgm.var)
plot(variogram$dist,variogram$gamma,xlab="Distance",ylab="Gamma")
lines(variogramLine(vgm.var,maxdist=Range*2))
})
panel
}

replot.smooth <- function(object) {
rp.tkrreplot(var.panel, plot)
object
}


var.panel <- rp.control("Variogram",Sill=20,Range=250,Nugget=0,Model="Mat",Cutoff=250,size=c(800,600))
rp.listbox(var.panel,Model,c("Mat","Sph","Gau"))
rp.slider(var.panel, Cutoff, 0,sqrt(areaSpatialGrid(grid)),showvalue=T)
rp.slider(var.panel, Sill, 0,var(data$Oxigen)*2,showvalue=T)
rp.slider(var.panel, Range, 0,sqrt(areaSpatialGrid(grid)),showvalue=T)
rp.slider(var.panel, Nugget, 0,var(data$Oxigen),showvalue=T)
rp.button(var.panel, title = "Fit",action=replot.smooth)
rp.tkrplot(var.panel, plot,variogram.plot)
rp.block(var.panel)


It probably need some more work in order to be perfect but at the moment I'm not interested in a perfect automatic script.
If someone has time to spend on it and is able to made it perfectly automatic for every data file, please share the script with the community.

By the way, I also implemented a line from the tcltk package for the selection of the data file interactively.
This is the resulting panel:

2 comments:

  1. Interactive Teaching Tools for Spatial Sampling, published in http://www.jstatsoft.org/v36/i13/paper is also fun!

    ReplyDelete
  2. Hi

    Thank you for the examples. i have a query. How do i make rpanel print the value of a calculation in the rpanel itself – say, just below a button? Assume that the panel has two rp.textentry boxes for variables x and y, and an rp.button and i want the sum of the 2 variables x and y to be printed below the submit button as say Answer: value

    ReplyDelete

Note: only a member of this blog may post a comment.