########################################################################
##
## Gj: 22-May-2000
##
## Splus commands used for Short Course on Splus, part 3.
##
########################################################################

## to invoke the spatial model:
module(spatial)
## and to unload:
module(spatial,unload=T)

## decrease printout width (for slides)
options(width=50)

########################################################################
##### EDA #####

### the scallops data ###
attach(scallops)  ## attach the data:

## plot sites on a map:
library('maps')  ## uue the map library in S-PLUS:
trellis.device()
map('usa',xlim=c(-74,-71),ylim=c(38.2,41.5))
points(long,lat)

## look at tcatch...
histogram(tcatch,nint=40)
histogram(tcatch[tcatch<100],nint=15)
log.catch <- log(tcatch+1)
histogram(log.catch,nint=10)

## create a ps-file:
trellis.device(postscript,file='hist.ps',horizontal=T,width=11*0.7,height=8*0.7)
histogram(log.catch,nint=10)
dev.off()


## interpolate log-catch and plot
int.data <- interp(long,lat,log.catch)
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
image(int.data,add=T)
image.legend(int.data,x=-72.75,y=38.8,s=c(2,0.2))
contour(int.data,add=T,col=4)
points(long,lat,pch=16,col=8,cex=0.8)

persp(int.data)

## create postscript:
trellis.device(postscript,file='2dim_scallops.ps',
               horizontal=F,width=6,height=6,color=T)
par(mar=c(0,0,0,0))
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
image(int.data,add=T)
image.legend(int.data,x=-72.75,y=38.8,s=c(2,0.2))
contour(int.data,add=T,col=2)
points(long,lat,pch=16,col=8,cex=0.4)
dev.off()

## trend:
trend.1 <- loess(log.catch ~ long*lat,span=0.75)
trend.2 <- loess(log.catch ~ long*lat,span=0.50)
trend.3 <- loess(log.catch ~ long*lat,span=0.25)
trend.4 <- loess(log.catch ~ long*lat,span=0.15)
anova(trend.1,trend.2,trend.3,trend.4)

trend <- loess(log.catch ~ long*lat)
trend

## create a grid:
long.g <- seq(min(long),max(long),le=50)
lat.g <- seq(min(lat),max(lat),le=50)
sc.grid <- expand.grid(long=long.g,lat=lat.g)

## use only grid-points within convex hull:
def.hull <- chull(long,lat)  ## T/F if site def the hull
my.hull <- list(x=long[def.hull], y=lat[def.hull])
in.hull <- points.in.poly(sc.grid$long,sc.grid$lat,
                          my.hull)

## use my own area:
in.area <- points.in.poly(sc.grid$long,sc.grid$lat,
                          my.area)

## plot the AOI and hull:
attach('/home/gardar/pub/.Data')  ## to access  my.area
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
points(sc.grid$long,sc.grid$lat,col=7,cex=0.4,pch=3)
points(long,lat,pch=16,col=8,cex=0.4)
polygon(my.area,lty=2,d=0,col=2)  ## the AOI
polygon(my.hull,lty=3,d=0,col=3)  ## the convex hull

## ps-file:
trellis.device(postscript,file='aoi.ps',
               horizontal=F,width=6,height=6,color=T)
par(mar=c(0,0,0,0))
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
points(sc.grid$long,sc.grid$lat,col=7,cex=0.4,pch=3)
points(long,lat,pch=16,col=8,cex=0.4)
polygon(my.area,lty=2,d=0,col=2)  ## the AOI
polygon(my.hull,lty=3,d=0,col=3)  ## the convex hull
dev.off()

## predict on grid:
pred <- predict(trend,sc.grid)
## put equal to NA outside of area-of-interest:
pred[!in.area] <- NA  ## NA outside of AOI

## plot trend:
trellis.device(postscript,file='pred_trend.ps',
               horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
persp(x=long.g,y=lat.g,pred,xlab="long",ylab="lat")
dev.off()


## trellis-plotting:
wireframe(pred ~ long * lat, data=sc.grid)
levelplot(pred ~ long * lat, data=sc.grid)


########################################################################
##### VARIOGRAM #####

## remove trend:
log.catch.res <- log(tcatch+1) - predict(trend)

## intrinsic?
int.res <- interp(long,lat,log.catch.res)
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
image(int.res,add=T)
contour(int.res,add=T,col=4)
points(long,lat,pch=16,col=8,cex=0.8)

## postscript file:
trellis.device(postscript,file='intrinsic.ps',
               horizontal=F,width=6,height=6,color=T)
par(mar=c(0,0,0,0))
map('usa',xlim=c(-73.5,-71),ylim=c(38.2,41.1))
image(int.res,add=T)
contour(int.res,add=T,col=4)
points(long,lat,pch=16,col=8,cex=0.4)
dev.off()

### empirical ###
v1 <- variogram(log.catch.res ~ loc(long,lat),method='r')
v1[1:5,]  ## show the first 5 lines
plot(v1,main='After trend-removing')  ## plot it
v2 <- variogram(log(tcatch+1) ~ loc(long,lat),method='r')
plot(v2)

trellis.device(postscript,file='omni_var.ps',horizontal=T,width=11*0.7,height=8*0.7,color=T)
par(mar=c(5,5,3,1))
plot(v1,main='After trend-removing')
dev.off()

## check for anisotropy:
v3 <- variogram(log.catch.res ~ loc(long,lat),
                azimuth=c(0,45,90,135),
                tol.azimuth=11.25, method='r')
plot(v3,main='After trend-removing')

trellis.device(postscript,file='dir_var.ps',horizontal=T,width=11*0.7,height=8*0.7,color=T)
par(mar=c(5,5,3,1))
plot(v3,main='After trend-removing')
dev.off()

v4 <- variogram(log.catch ~ loc(long,lat),
                azimuth=c(0,45,90,135),
                tol.azimuth=45/4, method='r')
plot(v4)

## anisotropy plot:
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=0,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=1*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=2*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=3*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=4*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=5*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=6*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))

## ps-file:
trellis.device(postscript,file='ani_var.ps',horizontal=T,width=11*0.7,height=8*0.9,color=T)
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=2*22.5,ratio=seq(1,by=0.25,le=12),
                method='r',layout=c(3,4))
dev.off()


anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=seq(0,180,le=9),ratio=2,
                method='r',layout=c(3,3))
anisotropy.plot(log.catch.res ~ loc(long,lat),
                angle=seq(0,180,le=9),ratio=2.5,
                method='r',layout=c(3,3),)

## correct for anisotropy:
v.f <- variogram(log.catch.res~loc(long,lat,a=45,r=2),
                 azimuth=c(0,45,90,135),
                 tol.azimuth=45/4, method='r')
plot(v.f,main='Angle=45, Ratio=2')

trellis.device(postscript,file='ani_cor_var.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
plot(v.f,main='Angle=45, Ratio=2')
dev.off()

## final variogram:
v.use <- variogram(log.catch.res ~ loc(long,lat,a=45,r=2),
                 method='r')
plot(v.use)

##### FITTING #####


model.variogram(v.use,fun=spher.vgram,range=0.8,sill=1.75-0.5,nugget=0.5)

spher.wfun <- function(gamma,distance,
                       np,range,sill,nugget)
{
  gamma.hat <- spher.vgram(distance,range,sill,nugget)
  sqrt(np) * (gamma/gamma.hat - 1)
}

## test:
spher.wfun(g=1.5,d=1,np=50,r=4,s=0.2,n=0.8)

## fit using nls:
fit.var <- nls(~spher.wfun(gamma,distance,
                           np,range,sill,nugget),
               data=v.use, start=list(range=0.6,
                             sill=2,nugget=1))

summary(fit.var)

## plot the fit:
plot(v.use)
dist.use <- seq(0,max(v.use$dis),le=100)
lines(dist.use,spher.vgram(dist.use,
      range=0.6819,sill=2.387,nugget=0.7013))
    
trellis.device(postscript,file='fit_var.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
plot(v.use)
dist.use <- seq(0,max(v.use$dis),le=100)
lines(dist.use,spher.vgram(dist.use,
      range=0.6819,sill=2.387,nugget=0.7013)) 
dev.off()

########################################################################
##### KRIGING #####

## use the krige() function:
k.fit <- krige(log.catch.res ~ loc(long,lat,a=45,r=2),
               covfun=spher.cov,
               range=0.6819,sill=2.387,nugget=0.7013)
k.fit

## predict:
k.pred <- predict(k.fit,sc.grid)
summary(k.pred)

## plot fitted:
k.pred$fit[!in.area] <- NA
levelplot(fit ~ long*lat,data=k.pred,
          contour=T,pretty=T)

trellis.device(postscript,file='fit_krig.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
levelplot(fit ~ long*lat,data=k.pred,contour=T,pretty=T)
dev.off()

## plot se:
k.pred$se.fit[!in.area] <- NA
levelplot(se.fit ~ long*lat,data=k.pred,
          panel=function(x,y,...) {
            panel.levelplot(x,y,...)
            panel.xyplot(long,lat,col=4,cex=1,pch=16)
          })

trellis.device(postscript,file='se_krig.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
levelplot(se.fit ~ long*lat,data=k.pred,
          panel=function(x,y,...) {
            panel.levelplot(x,y,...)
            panel.xyplot(long,lat,col=4,cex=0.5,pch=16)
          })
dev.off()


## get fitted values:
k.pred$pred <- predict(trend,sc.grid) + k.pred$fit
k.pred$pred[!in.area] <- NA  ## NA outside of AOI

## plot:
levelplot(pred ~ long*lat,data=k.pred,
          contour=T,pretty=T)

trellis.device(postscript,file='pred_krig.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
levelplot(pred ~ long*lat,data=k.pred,contour=T,pretty=T)
dev.off()

### cross-validation:

## function to do CV on a krige object
cv.krige <- function(obj,fit.data) {
  ## obj is the output from krige()
  n <- nrow(fit.data)
  fit.all <- predict(obj,newdata=fit.data)$fit
  fit <- se <- numeric(n)  
  for(i in 1:n) {
    cat("doing observation",i,"\n")
    new.krige <- update(obj,data=fit.data[-i,])  ## drop i
    pred <- predict(new.krige,newdata=fit.data[i,])
    fit[i] <- pred$fit[1]
    se[i] <- pred$se.fit[1]
  }
  res <- (fit.all - fit)/se
  data.frame(fit.all = fit.all, fit = fit,
             res = res, se = se)
}


## evaluate:
fit.data <- data.frame(log.catch.res,long,lat)
k.fit.cv <- cv.krige(k.fit,fit.data)

## plot:
histogram(k.fit.cv$res)
qn <- qqnorm(k.fit.cv$res,plot=F)
plot(qn$x,qn$y,type='n',xlab="Norm",ylab="CV")
text(qn$x,qn$y,1:148)
qqline(k.fit.cv$res)

trellis.device(postscript,file='cv_qq.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
qn <- qqnorm(k.fit.cv$res,plot=F)
plot(qn$x,qn$y,type="n",xlab="Norm",ylab="CV")
text(qn$x,qn$y,1:148,cex=0.8)
qqline(k.fit.cv$res)
dev.off()

## identify obs 68:
plot(long,lat,type='n')
points(long[68],lat[68],pch=16,col=2,cex=2)
text(long,lat,tcatch,cex=0.8)

trellis.device(postscript,file='obs_68.ps',horizontal=T,width=11*0.7,height=8.5*0.7,color=T)
plot(long,lat,type='n')
points(long[68],lat[68],pch=16,col=2,cex=1)
text(long,lat,tcatch,cex=0.5)
dev.off()