3D Anisotropy parameters

March 3, 2011
By

(This article was first published on R Video tutorial for Spatial Statistics, and kindly contributed to R-bloggers)

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)

To leave a comment for the author, please follow the link and comment on his blog: R Video tutorial for Spatial Statistics.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.