library(gstat)
library(RColorBrewer)

data(jura)

coordinates(jura.pred)<-~Xloc+Yloc

gridded(jura.grid) <- ~Xloc+Yloc

# indicator cosimulation for the 9 deciles of lead

q <- quantile(jura.pred$Pb, seq(.1,.9,.1))
# These are the actual values of the quantiles
q

# Create 9 spatial datasets for eaqch quantile
# beta is the mean of each variable.  Normally, we would not know the
# true mean of a dataset, but in this case we do, the mean of the 
# first indicator is 10%, for the second indicator is 20%, etc.

jura.i <- gstat(id = "Pb1", formula = I(Pb < q[1])~1, 
       data = jura.pred, nmax = 7, beta = .1, set = list(order = 4, zero = 1e-5))

jura.i <- gstat(jura.i, "Pb2", formula = I(Pb < q[2])~1, 
       data = jura.pred, nmax = 7, beta = .2)

jura.i <- gstat(jura.i,"Pb3", formula = I(Pb < q[3])~1, 
       data = jura.pred, nmax = 7, beta = .3)

jura.i <- gstat(jura.i,"Pb4", formula = I(Pb < q[4])~1, 
       data = jura.pred, nmax = 7, beta = .4)

jura.i <- gstat(jura.i,"Pb5", formula = I(Pb < q[5])~1, 
       data = jura.pred, nmax = 7, beta = .5)

jura.i <- gstat(jura.i,"Pb6", formula = I(Pb < q[6])~1, 
       data = jura.pred, nmax = 7, beta = .6)

jura.i <- gstat(jura.i,"Pb7", formula = I(Pb < q[7])~1, 
       data = jura.pred, nmax = 7, beta = .7)

jura.i <- gstat(jura.i,"Pb8", formula = I(Pb < q[8])~1, 
       data = jura.pred, nmax = 7, beta = .8)

jura.i <- gstat(jura.i,"Pb9", formula = I(Pb < q[9])~1, 
       data = jura.pred, nmax = 7, beta = .9)

# Create a semivariogram model
#	"One size fits all"
jura.i <- gstat(jura.i, model=vgm(1, "Sph", .3, 1), fill.all=T)

# Estimate the variogram of each indicator
x <- variogram(jura.i, cutoff=2)

# Fit the semivariogram for each indicator
# Note, this changes the nuggets and sills, but not the ranges
#  For the range, "one size fits all"
jura.fit = fit.lmc(x, jura.i)

plot(x,model=jura.fit)

# Simulate each indicator
z <- predict(jura.fit,newdata=jura.grid,nsim=2,indicators=TRUE)

# Instead of 9 indicators, I want one variable with each class value
# i.e. 1 for first class, 2 for second class, etc.

z$sim1 <- z$Pb1.sim1+z$Pb2.sim1+z$Pb3.sim1+z$Pb4.sim1+z$Pb5.sim1+z$Pb6.sim1+z$Pb7.sim1+z$Pb8.sim1+z$Pb9.sim1
z$sim1 <- 9-z$sim1
spplot(z,'sim1',at=seq(.5,9.5),col.regions=brewer.pal(9,'YlGnBu'))

z$sim2 <- z$Pb1.sim2+z$Pb2.sim2+z$Pb3.sim2+z$Pb4.sim2+z$Pb5.sim2+z$Pb6.sim2+z$Pb7.sim2+z$Pb8.sim2+z$Pb9.sim2
z$sim2 <- 9-z$sim2

# the data values are 
#	 0     1     2     3     4     5, etc, so I will have class breaks at
# -0.5   0.5   1.5   2.5   3.5   4.5

# spplot(z,'sim2',at=seq(-.5,9.5),col.regions=brewer.pal(9,'YlGnBu'))
# Unfortunately, brewer.pal does not have 10 class values
#  I could use the default color scheme, but those are garish.

# Plot the highest 9 classes (By default, the omitted class will be white)
spplot(z,'sim2',at=seq(.5,9.5),col.regions=brewer.pal(9,'YlGnBu'))




#################################
# indicator cosimulation for the lead at 25, 50, 75, 100, 125

data(jura)

coordinates(jura.pred)<-~Xloc+Yloc

gridded(jura.grid) <- ~Xloc+Yloc

qf <- ecdf(jura.pred$Pb)

jura.i <- gstat(id = "Pb1", formula = I(Pb < 25)~1, 
       data = jura.pred, nmax = 7, beta = qf(25), set = list(order = 4, zero = 1e-5))

jura.i <- gstat(jura.i, "Pb2", formula = I(Pb < 50)~1, 
       data = jura.pred, nmax = 7, beta = qf(50))

jura.i <- gstat(jura.i, "Pb3", formula = I(Pb < 75)~1, 
       data = jura.pred, nmax = 7, beta = qf(75))

jura.i <- gstat(jura.i, "Pb4", formula = I(Pb < 100)~1, 
       data = jura.pred, nmax = 7, beta = qf(100))

jura.i <- gstat(jura.i, "Pb5", formula = I(Pb < 125)~1, 
       data = jura.pred, nmax = 7, beta = qf(125))



# Note, the range for every indicator must be the same!
jura.i <- gstat(jura.i, model=vgm(1, "Sph", .3, 1), fill.all=T)

x <- variogram(jura.i, cutoff=2,width=.15)

# Fit the model
jura.fit = fit.lmc(x, jura.i)
#Plot the moel
plot(x,model=jura.fit)

z <- predict(jura.fit,newdata=jura.grid,nsim=4,indicators=TRUE)

# Create a new variable with the indicator class value
# i.e. 1 for first class, 2 for second class, etc.
z$sim1 <- z$Pb1.sim1+z$Pb2.sim1+z$Pb3.sim1+z$Pb4.sim1+z$Pb5.sim1
z$sim1 <- 5-z$sim1

z$sim2 <- z$Pb1.sim2+z$Pb2.sim2+z$Pb3.sim2+z$Pb4.sim2+z$Pb5.sim2
z$sim2 <- 5-z$sim2
z$sim3 <- z$Pb1.sim3+z$Pb2.sim3+z$Pb3.sim3+z$Pb4.sim3+z$Pb5.sim3
z$sim3 <- 5-z$sim3
z$sim4 <- z$Pb1.sim4+z$Pb2.sim4+z$Pb3.sim4+z$Pb4.sim4+z$Pb5.sim4
z$sim4 <- 5-z$sim4

# the data values are 
#	 0     1     2     3     4     so I will have class breaks at
# -0.5   0.5   1.5   2.5   3.5   4.5
spplot(z,c('sim1','sim2','sim3','sim4'),at=seq(-.5,5.5),col.regions=brewer.pal(6,'YlGnBu'))
# Unfortunately, I do not know how to label the plot to have the lead values labelled!