###################################################
########### CURSO KRIGING #########################
###################################################

rm(list = ls())
library(sf)
datos <- st_read("datos_un_mj.gpkg")
grilla <- st_read("grilla_un_mj.gpkg")

#### Distribucion de la muestra #####

summary(datos$vut_2019_usd)
hist(datos$vut_2019_usd)

library(RColorBrewer)
col = rev(brewer.pal(n=5, "Spectral"))

contorno <- st_read("contorno provincia.gpkg")
depto <- st_read("depto_un_mj.gpkg")

library(tmap)
mapa_muestra <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(datos) +
  tm_dots("vut_2019_usd",title = "Muestra VUT (USD)", palette = col ,
          style="quantile", n=5, alpha = 0.8) 

mapa_muestra

#### 1. KRIGING ORDINARIO #############

library (gstat)

variograma <- variogram(vut_2019_usd~1, datos, cressie=T, cutoff = 125000)
plot(variograma)
model_vario <- fit.variogram(variograma, vgm(c("Sph", "Gau", "Exp")))
plot(variograma, model_vario)
model_vario

ko <- krige(vut_2019_usd~ 1 , datos, st_centroid(grilla), model_vario)
grilla$ko <- ko$var1.pred

library(tmap)
library(RColorBrewer)
col = rev(brewer.pal(n=5, "Spectral"))

tmap_mode("view")
mapa_ko_vut <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(grilla) +
  tm_polygons("ko",title = "Prediccion KO", palette = col ,
              style="quantile", n=5, alpha = 0.8) + 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) 

mapa_ko_vut

summary(grilla$ko)
summary(datos$vut_2019_usd)

mapa_ko_var <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(ko) +
  tm_dots("var1.var", title = "Varianza", palette = col ,
          style="quantile", n=5, alpha = 0.8)+ 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5)

mapa_ko_var

#### 2. KRIGING UNIVERSAL #############
#### 2.1. KRIGING CON TENDENCIA #######

par(mfrow=c(1,2))
plot(datos$x,datos$vut_2019_usd, xlab="X", ylab="VUT")
lines(lowess(datos$x,datos$vut_2019_usd), col ="red", lwd=5)
plot(datos$y,datos$vut_2019_usd, xlab="Y", ylab="VUT")
lines(lowess(datos$y,datos$vut_2019_usd), col ="red", lwd=5)
par(mfrow=c(1,1))

lm <- lm(formula = vut_2019_usd~ x+I(y^2), data = datos)
summary(lm)

variograma <- variogram(vut_2019_usd~x+I(y^2), datos, cressie=T, cutoff = 125000)
plot(variograma)
model_vario <- fit.variogram(variograma, vgm(c("Sph", "Gau", "Exp")))
plot(variograma, model_vario)

kut <- krige(vut_2019_usd~ x+I(y^2), datos, st_centroid(grilla), model_vario)
grilla$kut <- kut$var1.pred
summary(grilla$kut)

tmap_mode("view")
mapa_kut_vut <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(grilla) +
  tm_polygons("kut",title = "Prediccion KU-T", palette = col ,
              style="quantile", n=5, alpha = 0.8) + 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) 

mapa_kut_vut

mapa_kut_var <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(kut) +
  tm_dots("var1.var", title = "Varianza", palette = col ,
          style="quantile", n=5, alpha = 0.8)+ 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(datos) + 
  tm_dots()

mapa_kut_var

#### 2.2. REGRESSION KRIGING ##########
names(datos)

library(PerformanceAnalytics) 
chart.Correlation(st_drop_geometry(datos), histogram=TRUE, pch=19) 

lm <- lm(formula =  vut_2019_usd ~ morganica + d_puerto + d_urbaniz + d_cacopio, 
        data = datos, na.action = na.omit)
grilla$pred_reg <- predict(object = lm, newdata = grilla)

datos$error_lm <- residuals(lm)

tmap_mode("view")
mapa_residuos <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(datos) +
  tm_dots("error_lm",title = "Residuos LM", palette = col ,
              style="quantile", n=5, alpha = 0.8) + 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) 

mapa_residuos

variograma <- variogram(residuals(lm)~1,datos, cressie=T, cutoff = 125000)
plot(variograma)
model_vario <- fit.variogram(variograma, vgm(c("Sph", "Gau", "Exp")))
plot(variograma, model_vario)
model_vario

kur <- krige(residuals(lm)~1,datos,  st_centroid(grilla), model_vario)

grilla$kur <- grilla$pred_reg + kur$var1.pred

summary(grilla$kur)

tmap_mode("view")

mapa_kur_vut <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(grilla) +
  tm_polygons("kur",title = "Prediccion KU-R", palette = col ,
              style="quantile", n=5, alpha = 0.8) + 
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) 

mapa_kur_vut

####### Comparacion entre predicciones ########

summary(datos$vut_2019_usd)
summary(grilla$ko)
summary(grilla$kut)
summary(grilla$kur)


###################################################
########### CURSO KKNN ############################
###################################################

rm(list=ls())

library(sf)
datos <- st_read("datos_un_mj.gpkg")
grilla <- st_read("grilla_un_mj.gpkg")

#### Distribucion de la muestra #####

summary(datos$vut_2019_usd)
hist(datos$vut_2019_usd)

library(RColorBrewer)
col = rev(brewer.pal(n=5, "Spectral"))

contorno <- st_read("contorno provincia.gpkg")
depto <- st_read("depto_un_mj.gpkg")

library(tmap)
mapa_muestra <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(datos) +
  tm_dots("vut_2019_usd",title = "Muestra VUT (USD)", palette = col ,
          style="quantile", n=5, alpha = 0.8) 

mapa_muestra


########## KKNN y KNN ##################
library(caret)
fitControl <- trainControl(method = "cv", number = 5)

library(doParallel)
cores=detectCores()
cl <- makeCluster(cores[1]) #not to overload your computer registerDoParallel(cl)
registerDoParallel(cl)

###### Diferencia en los parametros estimados
# Solo con las coordenadas
set.seed(9)
train_knn <- train(vut_2019_usd ~ x+y,
               data= datos,
               method = "knn",
               trControl = fitControl)
grilla$pred_knn <- predict(train_knn, newdata=grilla)

train_knn

set.seed(9)
train_kknn <- train(vut_2019_usd ~ x+y,
               data= datos,
               method = "kknn",
               trControl = fitControl)
grilla$pred_kknn <- predict(train_kknn, newdata=grilla)

# Diferencias
train_knn
train_kknn

library(RColorBrewer)
col = rev(brewer.pal(n=5, "Spectral"))
contorno <- st_read("contorno provincia.gpkg")
depto <- st_read("depto_un_mj.gpkg")
grilla2 = grilla

tmap_mode("view")
mapa_prediccion <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(grilla) + 
  tm_polygons("pred_knn", title = "Prediccion KNN", palette = col ,
              style="fixed", n=5, alpha = 0.8, breaks=c(2000,6000,8000,10000,12000,16000))+
  tm_shape(grilla2) + 
  tm_polygons("pred_kknn", title = "Prediccion KKNN", palette = col ,
              style="fixed", n=5, alpha = 0.8, breaks=c(2000,6000,8000,10000,12000,16000))
  
mapa_prediccion

### Continuamos trabajando con kknn 

##### Utilizando otras variables 

set.seed(9)
train_vbles <- train(vut_2019_usd ~ morganica + d_puerto + d_urbaniz + d_cacopio,
               data= datos,
               method = "kknn",
               trControl = fitControl)

grilla$prediccion3 <- predict(train_vbles, newdata = grilla)
summary(grilla$prediccion3)
hist(grilla$prediccion3)

# Utilizando variables y las coordenadas del punto
set.seed(9)
train_vblesc <- train(vut_2019_usd ~ morganica + d_puerto + d_urbaniz + d_cacopio + x + y,
               data= datos,
               method = "kknn",
               trControl = fitControl)
grilla$prediccion4 <- predict(train_vblesc, newdata=grilla)
summary(grilla$prediccion3)
summary(grilla$prediccion4)

# Mapas
library(RColorBrewer)
col = rev(brewer.pal(n=5, "Spectral"))
grilla2 = grilla

tmap_mode("view")
mapa_prediccion <-  tm_basemap(c(Satelite = "Esri.WorldImagery", Politico = "OpenStreetMap.DE")) +
  tm_shape(contorno) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(depto) + 
  tm_fill("white", alpha = 0)+
  tm_borders("black", lwd = 1.5) +
  tm_shape(grilla) + 
  tm_polygons("prediccion3", title = "Prediccion Vbles", palette = col ,
              style="fixed", n=5, alpha = 0.8, breaks=c(2000,6000,8000,10000,12000,16000))+
  tm_shape(grilla2) +
  tm_polygons("prediccion4", title = "Prediccion Vbles + Coord", palette = col ,
              style="fixed", n=5, alpha = 0.8, breaks=c(2000,6000,8000,10000,12000,16000))

mapa_prediccion