library(sp)
library(spacetime)
library(gstat)
setwd("/home//tchakounte//Bureau//Bureau.6.5.2014//nao")
dataSet <- read.csv("donneesannees.csv", sep=";")
str(dataSet)
# drop the last two empty rows
dataSet <- dataSet[1:39,]
str(dataSet)


# set up the spatial slot
spObj <- dataSet[,c("latitude", "longitude", "FEC")]
spObj
# make the data.frame spatial
coordinates(spObj) <- ~longitude+latitude
proj4string(spObj) <- CRS("+init=epsg:4324")

# Grille de représentation
#rr <- spObj[,"2015-01-01/2018-01-01"]
#rr <- as(rr,"STSDF")
#x1 <- seq(from=6,to=15,by=1)
#x2 <- seq(from=48,to=55,by=1)

x.range <- as.double(range(spObj@coords[,1]))
y.range <- as.double(range(spObj@coords[,2]))
x.range
y.range

#grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.01018847), y=seq(from=y.range[1], to=y.range[2], by=0.000532) )
#data(grd)
#coordinates(grd) <- c("x", "y")
#gridded(grd) <- TRUE

# plot for visual checking
spplot(spObj, "FEC", scales=list(draw=TRUE))

# set up the time object
timeObj <- as.POSIXct(c("2010-01-01","2011-01-01","2012-01-01","2013-01-01","2014-01-01"))
timeObj
# set up the data slot
dataObj <- data.frame(population=as.numeric(unlist(dataSet[c("pop2010", "pop2011", "pop2012", "pop2013", "pop2014")])))
dataObj
# build the spatio-temporal structure
stfdfObj <- STFDF(spObj, timeObj, dataObj)
stfdfObj

# plot for visual checking
stplot(stfdfObj, col.regions=bpy.colors())

# on Google Eareth if you like
library(plotKML)
plotKML(stfdfObj)

#######

DE_gridded <- SpatialPoints(cbind(rep(x.range,length(y.range)), rep(y.range,each=length(x.range))), 
                            proj4string=CRS(proj4string(stfdfObj@sp)))
gridded(DE_gridded) <- TRUE
DE_pred <- STF(sp=as(DE_gridded,"SpatialPoints"), time=stfdfObj@time)

# calculate and plot the empirical spatio-temporal variogram
empVgm <- variogramST(population~1, stfdfObj, tlags=0:3,cutoff=0.05)
plot(empVgm, wireframe=TRUE)
empVgm


# mind the number of points per bin: > 100
#empVgm$np
summary(empVgm)
# metric
sumMetricVgm <- vgmST("sumMetric", space=vgm( 4.4, "Lin", 196.6,  3), time =vgm( 2.2, "Lin",   1.1,  2),
                     joint=vgm( 34.6, "Exp", 136.6, 12),
                     stAni=51.7)
# Prediction
pred = krigeST(formula=population~1, data=stfdfObj, newdata = DE_pred, modelList=sumMetricVgm)
pred
gridded(pred@sp) <- TRUE
stplot(pred)