# (c) 2016 Dr. Karsten Schmidt - Geoscientific Senior Researcher -
# Collaborative Research Centre 1070 - ResourceCultures
# eScience Centre - University of Tübingen


library("sp")
library("deldir")
library("ggplot2")
library("rgeos")
library("geoR")
library("maptools")
library("rgdal")

# Functions
# Carson's Voronoi polygons function
voronoipolygons <- function(x) {
  require(deldir)
  require(sp)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  z <- deldir(crds[,1], crds[,2])
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                          y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
                                                                                       function(x) slot(x, 'ID'))))
}

# Settings
setwd("C:/ARBEIT/SCRIPTS/LEC")
PointDataSet = "Testdatensatz.shp"
Layer = "Testdatensatz"

#Final grid resolution
cellsize <- 500
# CountourLine steppings in m
ElevationBands = 1000

# Point dataset
ogrListLayers(PointDataSet)
Shape=readOGR(PointDataSet, layer=Layer)

#Reduce spatial redundancies
Shape <- remove.duplicates(Shape)

# Thiessen polygon tessellation
VoronoiPolygons <- voronoipolygons(Shape)
VoronoiPolygonNodes <- fortify(VoronoiPolygons) 

df.sp<-VoronoiPolygonNodes
coordinates(df.sp)<-~long+lat
VoronoiPolygonNodes.sp <- remove.duplicates(df.sp)

VoronoiPolygonNodes.sp.df <- as.data.frame(VoronoiPolygonNodes.sp)
VoronoiPolygonNodes.sp.df$NearestNodeDistance <- -9999

#Report
plot(VoronoiPolygons, lty ="dashed")
points(VoronoiPolygonNodes.sp, pch=16, col="red")
points(Shape, pch = 16, col="darkblue")

#Calculate distance between voronoi nodes and arch. sites 
glDistance <- gDistance(VoronoiPolygonNodes.sp, Shape, byid = TRUE) # will be 0, all inside


##Create GRID
# Set cellsize
xmin <- min(VoronoiPolygonNodes$long) - cellsize
xmax <- max(VoronoiPolygonNodes$long) + cellsize
ymin <- min(VoronoiPolygonNodes$lat) - cellsize
ymax <- max(VoronoiPolygonNodes$lat) + cellsize

# Create grid (location where an interpolation shouild be calculated for)
grd <- expand.grid(x=seq(from=xmin, to=xmax, by=cellsize), y=seq(from=ymin, to=ymax, by=cellsize))

for (pts in 1:length(VoronoiPolygonNodes.sp.df[,1])) {
  points(Shape[which.min(glDistance[,pts]),], pch=16, cex=2, col="darkblue")
  points(VoronoiPolygonNodes.sp[pts,],pch=16, cex=2, col="red")

  #Sys.sleep(0.15)
  VoronoiPolygonNodes.sp.df[pts,length(VoronoiPolygonNodes.sp.df[1,])] <- min(glDistance[,pts])/1000 # km
}


###### Ordinary Kriging using geoR-Package
m.coords <- as.matrix(cbind(VoronoiPolygonNodes.sp.df[1],VoronoiPolygonNodes.sp.df[2]))
v.values <- as.vector(VoronoiPolygonNodes.sp.df$NearestNodeDistance)
bin1 <- variog(coords=m.coords, data=v.values)

# Show Variogram 
plot(bin1, pch= 19, cex.lab= 1.5)
ols <- variofit(bin1, ini=grd, fix.nug=TRUE, wei="npairs")
summary(ols)
lines(ols, lty=2)

km <- krige.conv(coords=m.coords, data=v.values, loc=grd, krige=krige.control(obj.model=ols))
image(km, axes=F, xlab="", ylab= "", col =  terrain.colors(40))
axis(2, cex = 2)
axis(1, cex = 2)

# Add contour lines to graph
contour(km, add=TRUE, nlev=15, lty = "dashed", labcex = 1)

# Export Grid and Contour Lines (ESRI-Format)
ExportGrid <- expand.grid(x=seq(from=xmin, to=xmax, by=cellsize), y=seq(from=ymin, to=ymax, by=cellsize))
coordinates(ExportGrid) <- ~x+y
gridded(ExportGrid) <- TRUE
ExportGrid$predict <- km$predict

ExportGrid$predict <- as.numeric(ExportGrid$predict)
writeGDAL(ExportGrid["predict"],"KrigingMap.tif",drivername="GTiff", type="Float64", options=NULL)

im<- as.image.SpatialGridDataFrame(ExportGrid)
cl <- contourLines(im,ElevationBands)
shp <- ContourLines2SLDF(cl)

SLDF <- ContourLines2SLDF(cl)
writeOGR(SLDF, ".", "KrigingMap_ContourLines", driver="ESRI Shapefile")