library(maptools)
library(rgdal)
library(RSAGA)
library(gstat)
SA.aea <- "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs"

# load the gridded and point data:
load(url("http://spatial-analyst.net/book/system/files/WISE_SOC.Rdata"))
SA.bbox <- gridsSA@bbox
writeOGR(GSPD.aea, "GSPD_aea.shp", ".", "ESRI Shapefile")
# subset the points:
rsaga.geoprocessor(lib="shapes_tools", module=14, param=list(SHAPES="GSPD_aea.shp", 
CUT="m_GSPD_aea.shp", METHOD=0, TARGET=0, CUT_AX=SA.bbox[1,1], 
CUT_BX=SA.bbox[1,2], CUT_AY=SA.bbox[2,1], CUT_BY=SA.bbox[2,2]))

# ------------------------------------------------------------
# Variogram modelling
# ------------------------------------------------------------

# original variable:
m_SOC <- readShapePoints("m_GSPD_aea.shp", CRS(SA.aea))
m_SOC.ov <- overlay(gridsSA, m_SOC)
m_SOC.ov$SOC <- m_SOC$SOC
m_SOC.ov <- remove.duplicates(m_SOC.ov)  # there are duplicate points!
sel <- !is.na(m_SOC.ov$SOC)&!is.na(m_SOC.ov$SOC_reg)
SOC.svar <- variogram(log1p(SOC)~1, m_SOC.ov[sel,])
SOC.ovgm <- fit.variogram(SOC.svar, vgm(nugget=0, model="Exp", range=80000, sill=var(log1p(m_SOC.ov$SOC), na.rm=T)))
SOC.ovgm.plt <- plot(SOC.svar, SOC.ovgm, pc="+", cex=1.5, pl=TRUE, col="black", main="Original variable")

# residuals:
res_SOC.svar <- variogram(log1p(SOC)~SOC_reg, m_SOC.ov[sel,])
SOC.rvgm <- fit.variogram(res_SOC.svar, vgm(nugget=var(m_SOC.ov$SOC, na.rm=T)/2, 
model="Exp", range=80000, sill=var(m_SOC.ov$SOC, na.rm=T)/2))
SOC.rvgm.plt <- plot(res_SOC.svar, SOC.rvgm, pc="+", cex=1.5, pl=FALSE, col="black", main="Residuals")
SOC.rvgm

# synchronize the two plots:
SOC.rvgm.plt$x.limits <- SOC.ovgm.plt$x.limits
SOC.rvgm.plt$y.limits <- SOC.ovgm.plt$y.limits
print(SOC.ovgm.plt, split=c(1,1,2,1), more=TRUE)
print(SOC.rvgm.plt, split=c(2,1,2,1), more=FALSE)

# regression-kriging:
m_SOC.rk <- krige(log1p(SOC)~SOC_reg, m_SOC.ov[sel,], gridsSA, SOC.rvgm, nmin=30, nmax=40, block=c(10e3, 10e3))
# back-transform:
m_SOC.rk$SOC_rk <- expm1(m_SOC.rk$var1.pred)

# Correct the map where the density of points is low using the USDA map:
download.file("http://spatial-analyst.net/worldmaps/SOC.zip", destfile=paste(getwd(), "/SOC.zip", sep=""))
unzip("SOC.zip")
# resample to the same grid:
rsaga.esri.to.sgrd(in.grids="SOC.asc", out.sgrd="SOC.sgrd", in.path=getwd())
rsaga.geoprocessor(lib="pj_proj4", 2, param=list(SOURCE_PROJ="\"+proj=longlat +datum=WGS84\"", 
TARGET_PROJ=paste('"', SA.aea ,'"', sep=""), SOURCE="SOC.sgrd", TARGET="m_SOC_USDA.sgrd", TARGET_TYPE=2, 
INTERPOLATION=1, GET_SYSTEM_SYSTEM_NX=m_SOC.rk@grid@cells.dim[[1]], GET_SYSTEM_SYSTEM_NY=m_SOC.rk@grid@cells.dim[[2]], 
GET_SYSTEM_SYSTEM_X=m_SOC.rk@grid@cellcentre.offset[[1]], GET_SYSTEM_SYSTEM_Y=m_SOC.rk@grid@cellcentre.offset[[2]], 
GET_SYSTEM_SYSTEM_D=10000))
rsaga.sgrd.to.esri(in.sgrds="m_SOC_USDA.sgrd", out.grids="m_SOC_USDA.asc", out.path=getwd(), prec=3)
m_SOC.rk$SOC_USDA <- readGDAL("m_SOC_USDA.asc")$band1

# merge the two maps (BCSP formula):
w <- sqrt(m_SOC.rk$var1.var)/sqrt(var(log1p(m_SOC.ov$SOC), na.rm=T))
m_SOC.rk$w <- 1-w/max(w, na.rm=TRUE)
m_SOC.rk$SOC.f <- round((m_SOC.rk$w * m_SOC.rk$SOC_rk+(1-m_SOC.rk$w) * m_SOC.rk$SOC_USDA), 0)

spplot(m_SOC.rk[c("SOC_rk","SOC_USDA", "SOC.f")], col.regions=bpy.colors(), at=seq(0,60,5))