IntR - Interactive GUI for performing geostatistical analysis in R

In 2011 I presented at the UseR conference, held in Warwick (UK), a piece of software for easing the learning curve of R for geostatistical analysis.

It is a very simple attempt to create an interactive interface in R, by using Python as GUI. From the GUI the user can easily create a script to perform some basic geostatistical analysis, such as ordinary kriging and others.

It should be fairly easy to use, even though the user needs to have some basic knowledge of GIS layers and geostatistics in order to be successful. On the other hand, I know a lot of users that are approaching R and perhaps GIS for the first time, and I think this software can be of great help in particular for these users. After a quick phase in which they need to learn some basics, such as what a ASCII grid is, I think this software should ease the learning curve of R itself, which according to different books is pretty steep.

The software is available from here: IntR - Home Page
Here you can find two versions of IntR, one for 2D geostatistics and the other for 3D stuff, plus a sample dataset and manuals.

It is completely free and you can try it and use it as long as you want, you can also modify its source code if you want to improve its functionalities. However, consider it to be in a beta-testing phase so always check its output and be careful in relying of it. I do not accept any responsibilities if the software does not work as expected or fails to produce a correct output.

Below, I attached an image the original abstract with the work I presented at the UseR conference of 2011 (the pdf is available here at page 91):

Below I put the video of the power point presentation I gave in Warwick, plus a short video tutorial (better to watch them at 720p)regarding the use of IntR for Ordinary kriging. You should also take a look at the IntR manual for more info on how to use it.

 

3D Anisotropy parameters

Hello everyone,

I’m trying to create a method for calculating the 5 parameters for the 3D anisotropy in an automatic way.

Basically I created a loop for analysing the range of the horizontal variogram in every direction and extrapolating the maximum (for the angle p). Then it computes the vertical variogram in the direction of maximum continuity and extrapolates the maximum range (for the angle q).

Subsequently, it computes the vertical variogram for the direction of minimum continuity (angle of max continuity + 90), and extrapolates the maximum (for the angle r).

Ultimately, it calculates the ratio between the range of the maximum continuity and the other 2.

I would like to know what the community think of this idea.

Could you please tell me if it makes sense and, if it does, how to improve it?

So far I wrote the following lines:

library(gstat, pos=4)

library(sp, pos=4)

library(lattice)

c1<-read.asciigrid("dtm.asc")

##a covariate that will be used to specify the cutoff of the variogram

data<-read.table("data.txt",sep="",header=T)

##Data format:

## Lat = latitude

## Lon = longitude

## depth = depth

## value = data value

coordinates(data)=~Lat+Lon+depth

##Anisotropy

##First Loop for angle p

mod<-vgm(var(data$value),"Sph",sqrt(areaSpatialGrid(c1)),0)

write.table(matrix(c("Angle","Range"),ncol=2),"hor.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=sqrt(areaSpatialGrid(c1)),width=sqrt(areaSpatialGrid(c1))/15,alpha=c(i),tol.hor=10),mod,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"hor.txt",append=T,row.names=F,col.names=F)}

hor<-read.table("hor.txt",sep="",header=T)

xyplot(Range~Angle,data=hor,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range<-subset(hor,hor$Range==max(hor$Range))

##Second Loop for angle q

mod_depth<-vgm(var(data$value),"Mat",max((data@coords)[,3]),0)

write.table(matrix(c("Angle","Range"),ncol=2),"ver1.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=max((data@coords)[,3]),width=max((data@coords)[,3])/length((data@coords)[,3]),alpha=c(max(max_range$Angle)),beta=c(i),tol.hor=10,tol.ver=10),mod_depth,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"ver1.txt",append=T,row.names=F,col.names=F)}

ver1<-read.table("ver1.txt",sep="",header=T)

xyplot(Range~Angle,data=ver1,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range1<-subset(ver1,ver1$Range==max(ver1$Range))

##Third Loop for angle r

mod_depth<-vgm(var(data$value),"Mat",max((data@coords)[,3]),0)

write.table(matrix(c("Angle","Range"),ncol=2),"ver2.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=max((data@coords)[,3]),width=max((data@coords)[,3])/length((data@coords)[,3]),alpha=c(max(max_range$Angle)+90),beta=c(i),tol.hor=10,tol.ver=10),mod_depth,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"ver2.txt",append=T,row.names=F,col.names=F)}

ver2<-read.table("ver2.txt",sep="",header=T)

xyplot(Range~Angle,data=ver2,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range2<-subset(ver2,ver2$Range==max(ver2$Range))

##Parameters

p=max(max_range$Angle)

q=max(max_range1$Angle)

r=max(max_range2$Angle)

s=max(max_range$Range)/max(max_range1$Range)

t=max(max_range$Range)/max(max_range2$Range)