options(scipen=3)
options(digits=15)
rm(list=ls())
require(sp)
require(gstat)
require(maptools)
require(automap)
require(moments)

cellsize <- 50

#                 1           2               3             4
out=paste("Variable","Model","C0","C","Spatial","Range","Kappa","R2","p","SSE","CV.me","CV.rmse","CV.msne","CV.r","CV.rr",sep=",")
setwd("C:\\Rcode")
data1 <- read.csv("yongji2000-GK.csv", header = TRUE, sep = ",")
data1$OLSENP <- data1$OLSENP/2.29    # convert P2O5 to P
data1$EXTK   <- data1$EXTK/1.2       # convert K2O  to K
names(data1)
class(data1)
coordinates(data1) <- ~X+Y
proj4string(data1) <- CRS("+proj=tmerc +lat_0=0 +lon_0=111 +k=1.0 +x_0=500000 +y_0=0 +ellps=krass +units=m +no_defs +towgs84=22,-118,30.5,0,0,0,0")

#read in a shape file of the boundary
aSHAPE <- readShapePoly("yjbdtm.shp", IDvar="YJBDTM_ID", proj4string=CRS(proj4string(data1))) 

#generate the grid coordinates (in meters)
xLL <- round(slot(aSHAPE,"bbox")[1]) - 100
xUR <- round(slot(aSHAPE,"bbox")[3]) + 100
yLL <- round(slot(aSHAPE,"bbox")[2]) - 100
yUR <- round(slot(aSHAPE,"bbox")[4]) + 100
xN <- round((xUR-xLL)/cellsize)
yN <- round((yUR-yLL)/cellsize)
data.grid = SpatialGrid(GridTopology(c(xLL,yLL),c(cellsize,cellsize),c(xN,yN)))
proj4string(data.grid) <- CRS(proj4string(aSHAPE))

#OLSENP
data <- data1
   OLSENP_MAX      <- 46 #max(data$OLSENP)*1.1   # 45
   OLSENP_MIN      <- 0 #min(data$OLSENP)*0.9     # 0.5
   OLSENP_DIFF     <- OLSENP_MAX - OLSENP_MIN
   data$OLSENPt    <- log(((data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)/(1-(data$OLSENP-OLSENP_MIN)/OLSENP_DIFF))

   kurtosis(data$OLSENP)
   kurtosis(data$OLSENPt)
   kurtosis(log(data$OLSENP))
   skewness(data$OLSENP)
   skewness(data$OLSENPt)
   skewness(log(data$OLSENP))
   
   #comment: log-transformation is the best
 
   x11(width=8, height=6.5)
   par(mfrow=c(1,2))
   hist(log(data$OLSENP), main="", xlab="log(OLSENP)", n=25)
   boxplot(log(data$OLSENP), xlab="log(OLSENP)")

   #estimate semivariogram model                          1  2   3   4   5   6    7   8   9
#   soil.v <- autofitVariogram(log(OLSENP)~1, data)$exp_var#[c(-1,-3,-4), ]
   soil.v <- variogram(log(OLSENP)~1, data, boundaries=c(100,200,300,400,600,900,1200,1400,1600,1800,2000,4000,6000,8000,10000,15000,20000))
#   v      <- variogram(log(OLSENP)~1, data, boundaries=c(100,200,300,400,600,900,1000,1500,2000))
#   soil.v <- variogram(log(OLSENP)~1, data, boundaries=c(100,200,300,400,600,900,1000,1500,2000))[-6, ]
   
   #soil.v <- autofitVariogram(log(OLSENP)~1, data)$exp_var[c(-1,-3,-4), ]
   vmodel <- vgm(0.5, "Sph", 15000, add.to=vgm(0.4, "Exp", 100, nugget=0.2))
#   vmodel <- vgm(0.5, "Exp", 1000, 0.01)
   v.sph <- fit.variogram(object=soil.v, model=vmodel, fit.sills=T, fit.ranges=T, fit.method=7)
   v.sph
   x11(width=8, height=6.5)
   par(mfrow=c(1,1))
   plot(soil.v, v.sph, col="black", ylim=c(0,0.5), cex=1, lwd=2.0, main="log(OLSENP)", xlab="Separation distance (m)", ylab="Semivariance", plot.numbers=F)

model <- v.sph$model[2]
c0    <- v.sph$psill[1]
c1    <- v.sph$psill[2]
range <- v.sph$range[2]
kappa <- v.sph$kappa[2]
fitted_gamma <- variogramLine(v.sph, dist_vector=soil.v$dist)

tcor <- cor.test(soil.v$gamma, fitted_gamma$gamma)
R2   <- (tcor$estimate[[1]])^2
R2
p    <- tcor$p.value[[1]]
sse  <- attributes(v.sph)$SSErr

#f <- function(idp, formula, data, ...) 
#     sqrt(mean(krige.cv(formula, data, set=list(debug=0,idp=idp), ...)$residual^2))
#optimize(f, interval=c(0.01,4), formula=OLSENPt~1, data=data)
#$minimum
#[1] 1.31
#$objective
#[1] 0.87

#cross-validation
#cv <- krige.cv(log(OLSENP)~1, data, model=NULL, set=list(idp=1.31), nmax=6, nfold=10)
cv <- krige.cv(log(OLSENP)~1, data, model=v.sph, nfold=10)
summary(cv)
# mean error, ideally 0
cv.me <- round(mean(cv$residual),6)
cv.me
# MSPE, ideally small
round(mean(cv$residual^2),6)
# RMSE, ideally small
cv.rmse <- round(sqrt(mean(cv$residual^2)),6)
cv.rmse
# Mean square normalized error, ideally close to 1
cv.msne <- round(mean(cv$zscore^2),6)
cv.msne
# correlation observed and predicted, ideally 1
cv.r <- round(cor(cv$observed, cv$observed-cv$residual),6)
cv.r
# correlation predicted and residual, ideally 0
cv.rr <- round(cor(cv$observed-cv$residual, cv$residual),6)
cv.rr

out[length(out)+1]=paste("log(OLSENP)", model, c0, c1, round(c0/(c0+c1),2), range, kappa, R2, p, sse, cv.me, cv.rmse, cv.msne, cv.r, cv.rr, sep=",")

# idw
#data.gstat <- gstat(id="OLSENPt", formula=OLSENPt~1, data=data, nmax=20, model=NULL, set=list(idp=1.31))

#ordinary kriging
data.gstat <- gstat(id="log_OLSENP", formula=log(OLSENP)~1, data=data, model=v.sph)

#Predicting surface
data.grid <- predict(data.gstat, data.grid)
names(data.grid)
class(data.grid)

#backward transformation
#data.grid$OLSENP.pred <- exp(data.grid$OLSENPt.pred)/(1+exp(data.grid$OLSENPt.pred))*OLSENP_DIFF+OLSENP_MIN
data.grid$OLSENP.pred <- exp(data.grid$log_OLSENP.pred)

#clip the whole grid using Hongtong boundary shape file.
data.grid.clip           <- data.grid[!is.na(overlay(data.grid, aSHAPE)), ]
gridded(data.grid.clip) <- TRUE

#output to ascii ARC grid
data.grid.clip$OLSENP.pred <- round(data.grid.clip$OLSENP.pred, 2)
write.asciigrid(data.grid.clip["OLSENP.pred"], paste("automap\\OLSENP_OK.asc", sep=""), attr=1, na.value=-9999)

#output display
l1 <- list("sp.polygons", col="black", aSHAPE)
l2 <- list("SpatialPolygonsRescale", layout.north.arrow(), offset=c(xUR- 3000, yLL+2000), scale=3000)
l3 <- list("SpatialPolygonsRescale", layout.scale.bar(),   offset=c(xUR-17000, yLL+2000), scale=10000, fill=c("transparent","black"))
l4 <- list("sp.text", c(xUR-17000, yLL+4000), "0"    )
l5 <- list("sp.text", c(xUR- 7000, yLL+4000), "10 km")
x11(width=8, height=6.5)
par(mfrow=c(1,1))
spplot(data.grid.clip["OLSENP.pred"], col.regions=bpy.colors(256), sp.layout = list(l1,l2,l3,l4,l5))


#write the output file
write(out,file="automap\\Models-gstat-variofit-county-automap.csv")