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")