[SOLVED] Isotropic Variogram Map

Here you can report the problems encountered when using a function relative to the variography:
- experimental variography (vario.calc, vario.grid, vmap.calc, ...)
- model fitting (model.auto)

[SOLVED] Isotropic Variogram Map

Postby Celeste » Tue Feb 19, 2013 5:13 am

I created a simulated isotropic dataset and the variogram looks good but the map doesn't look good and maybe I am just using the wrong inputs. It's all the same color basically. Here is my code and a picture.
It won't let me upload a .R or .docx or .txt for my code so here it is, sorry its long.


Code: Select all
# This is 2nd program to Generate n by n (equally spaced) Gaussian Random Field


#Initialize parameter values and matrices

gamma=0.20 #Was 0.20 for isotropy

n=50

m=3 #constant mean value
sig=1 #constant variance
G=matrix(0,n^2,n^2) #Gamma

mu=rep(m,n^2) #Put Mean function Here

V=sig*diag(n^2)  #Variances, V matrix
I=diag(n^2)   #Identity matrix

#Create Gamma matrix
for (r in 1:(n^2-1))
{
  G[r,(r+1)]=gamma
  G[(r+1),r]=gamma
}
for (r in 1:(n^2-n))
{
  G[r,(r+n)]=gamma
  G[(r+n),r]=gamma
}

#now remove n and n+1 neighbors, etc.

for ( r in 1:(n-1))
{
  G[r*n,r*n+1]=0
  G[r*n+1,r*n]=0
}

C=I-G
C=solve(C)


#Now how strong is the correlation?

#Consider the variance of the mean, if gamma=0, Varmean=sig
Varmean=sum(C)/n^2

#Now create (I-Gamma)^{-1}*V
C=C%*%V

#Symmetrize C to allow for Cholesky decomposition
C=(C+t(C))/2

#Now generate the data

L=chol(C)
eps=rnorm(n^2,mean=0,sd=1)
z=mu+L%*%eps

#Now place vector z onto nXn grid
Z=matrix(0,n,n)
for (i in 1:n)
{
  for (j in 1:n)
  {
    Z[i,j]=z[(i-1)*n+j,1]
  }
}


###Convert to Database (http://r.789695.n4.nabble.com/Converting-matrix-data-to-a-list-td3050208.html)
mat<-as.data.frame.table(as.matrix(Z))
colnames(mat) <- c("X", "Y","AU") #(http://stackoverflow.com/questions/6081439/changing-column-names-of-a-data-frame-in-r)
mat
#change letters to numbers for X
mat$X<-as.factor(mat$X)
mat$X
levels(mat$X)<-1:n
mat$X
as.numeric(mat$X)
mat
#change letters to numbers for Y
mat$Y<-as.factor(mat$Y)
mat$Y
levels(mat$Y)<-1:n
mat$Y
as.numeric(mat$Y)
mat

p <- ggplot(mat, aes(X,Y,colour=AU)) +geom_point()+scale_colour_gradientn(colours=c('blue', 'cyan', 'yellow', 'red'))
p<- p+ labs(title = "Plan View for Au Real Dataset with Color Scale")+ xlab("Easting") + ylab("Northing")
p



library(RGeoS)
library(e1071)
library(gstat)
library(fBasics)
library(R2HTML)
library(xtable)
library(geoR)
coords.aniso(coords, aniso.pars, reverse = FALSE)

#Database Creation
goldisodata.db<-db.create(mat,flag.grid=FALSE,ndim=2,autoname=F)
goldisodata.db

plot(goldisodata.db,pch=21,bg.in="black")
# Omnivariogram
#gold1.varioiso0<-vario.calc(goldisodata.db,lag=1,nlag=10,dirvect=c(0,30,60,90,120,150),calcul="vg")
gold1.varioiso0<-vario.calc(goldisodata.db,lag=1,nlag=10,calcul="vg")
gold1.varioiso0
plot(gold1.varioiso0,npairdw=TRUE,npairpt=TRUE)
title(main="OmniVariogram for Isotropic Dataset",xlab="Distance (feet)",
ylab="Variance")

# Experimental Variograms
gold1.varioisoE<-vario.calc(goldisodata.db,lag=1,nlag=10,dirvect=c(0,30,60,90,120,150),calcul="vg")
gold1.varioisoE
plot(gold1.varioisoE,npairdw=TRUE,npairpt=TRUE)
title(main="Experimental Variograms for Isotropic Dataset",xlab="Distance (feet)",
ylab="Variance")

goldiso.model<-model.auto(gold1.varioiso0,struct=c("Nugget Effect","Exponential"),
title="Modelling Variogram with Nugget Effect and Exponential Structure for Fake Isotropic data",
xlab="Variance",ylab="Distance")

#Isotropic
cat("Simtub\n")
#gridi <- db.create(flag.grid=T,nx=c(100,100,100))
gridi<-db.create(flag.grid=T,nx=c(100,100))
grid2i <- simtub(,gridi,goldiso.model,nbtuba=1000)

cat("Variogram Maps\n")
vmapiso <- vmap.grid(grid2i,nx=100,ny=100,calcul="vg")
plot(vmapiso,title="Variogram Map for the Fake Isotropic Dataset")
Attachments
isotropy.JPG
Isotropy Vario Map
isotropy.JPG (44.45 KiB) Viewed 2801 times
Celeste
 
Posts: 20
Joined: Sun Oct 14, 2012 8:10 pm

Re: Isotropic Variogram Map

Postby Didier Renard » Wed Feb 20, 2013 10:02 am

Hi Celeste

Back to question on the Variogram map ... as I see.
I ran your (long) script.
The principle is to generate a model based on an automatic fit performed on a gold data base.
This gold data set itself has been simulated (I think).

Then, using this (isotropic) model, you perform a non-conditional simulation on a grid of 50x50 units (say meters)
Finally you calculate the variogram map over this grid (the variogram map dimension is 100x100m) and visualize it. This is the figure
that you attached to your post.

I checked the contents of your variogram. It is composed of:
- a nugget effect (sill=0.454)
- an exponential with a range of 3.242 and a sill of 0.831
Therefore, on a total sill of 1.284, you have 35% of nugget.

I first visualize the variogram (traditional visualization) over a distance of 100m:

plot(goldiso.model,xlim=c(0,100))

You will see that the range of 3m is very small compared to size of your VMAP (i.e. 100m).
Then I also visualized the theoretical variogram map that you obtain with your model:

plot(model.grid(goldiso.model,nx=c(100,100)))

You hardly see anything but a small dot at the center of the plot.
In order to zoom, you can try:

plot(model.grid(goldiso.model,nx=c(10,10)))

Then you understand that the center is a very narrow hole (or a pick in covariance) whose extension is
given by the range of your model.
So, as a conclusion, I will say that the variogram map that you obtained experimentally from your simulated data
REFLECTS PERFECTLY the contents of your mail: it is essentially flat (equal to the total sill) when you
have passed the distance from the center equal to the range (i.e. 3m).

A final view which may help you understand the VMAP output:

plot(vmapiso,name.persp=5,theta=30,phi=20)

Hope this will help.
Didier Renard
 
Posts: 337
Joined: Thu Sep 20, 2012 4:22 pm

Re: Isotropic Variogram Map

Postby Celeste » Wed Feb 20, 2013 9:12 pm

I will give it a try tonight. Thank you :)
Celeste
 
Posts: 20
Joined: Sun Oct 14, 2012 8:10 pm


Return to Variography

Who is online

Users browsing this forum: No registered users and 2 guests

cron