[R-sig-Geo] Gaussian transformation

[acanion]

Carlo,

I am also interested in this because I'm working with groundwater nitrate
data that is zero-inflated.  

I wrote some code to perform the Gaussian anamorphosis in RGeostats (see the
example below with Meuse data), then do the kriging and simulation in gstat,
and finally back-transform the results using RGeostats.  I am not a serious
geostatistics person, so it would be nice if someone could confirm that this
is ok to do.  

My main issue is that the gaussian anamorphosis fit doesn't really seem to
work with this zero inflated data, altough I've seen references that say it
should 
(e.g, http://www.sciencedirect.com/science/article/pii/S0269749104002738).  
Any thoughts from the community here?

Anyway, here is my example code with the meuse data set.  The final result
is a data frame with x/y coordinates and each simulation result.  If you run
enough simluations, you can calcuate probabilities of exceeding a
concentration at each point.  

#01. Fit Gaussian Anamorphosis 
mdb = meuse[,c("x","y","zinc")] #must be 3 columns: x,y,z
mdb.rgdb = db.create(mdb,ndim=2,autoname=F)
mdb.herm = anam.fit(mdb.rgdb,name="zinc",type="hermitian")
mdb.hermtrans =
anam.z2y(mdb.rgdb,names=db.getname(mdb.rgdb,"z",1),anam=mdb.herm)
zinc.trans = mdb.hermtrans at items$Gaussian.zinc

#02. Variogram fit to gaussian transformed data using gstat 
meuse$zn.trans= zinc.trans
coordinates(meuse) = ~x+y
zn.svgm <- variogram(zn.trans~1,data=meuse)
vgm1 = vgm(1.17,"Sph",1022, 0.0925)
plot(zn.svgm,vgm1)
meuse.gstat <- gstat(id="zn.trans", formula=zn.trans~1, model=vgm1,
data=meuse,nmin=2, nmax=30)

#03. Kriging predictions with 4 simulations for simplicity (gstat)
data(meuse.grid)
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid)=T
nsim = 4
zntrans.ok = predict.gstat(meuse.gstat,meuse.grid,nsim=nsim) #multiple
simulations

#04. Back transform using RGeostats
zn.bts = cbind(coordinates(zntrans.ok),zntrans.ok at data)
zn.bts.db = db.create(zn.bts,autoname = F)
transout = zn.bts
transout[,3:length(transout)]=NA

  #This loops through each  item in the RGeostats DB and backtransform it
for(i in 4:(nsim+4)){
  
  sim = db.getname(zn.bts.db,rank = i)
  simtrans = paste("Raw.",sim,sep="")
  tempdb = anam.y2z(zn.bts.db,names=sim,anam = mdb.herm)
  transout[,i-1] = eval(parse(text =
paste("tempdb at items$",simtrans,sep="")))
  
} 











--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Gaussian-transformation-tp7588126p7588459.html
Sent from the R-sig-geo mailing list archive at Nabble.com.