setwd("~/Documents/GitLab/Valizas-River-Bathymetry-Model/")

# Manipulación de datos
library("data.table")
library("dplyr")

# Reestructurar datos
library("reshape2")

# Conexión a Bases de Datos
library("RPostgreSQL")
library("postGIStools")

# Manipulación de objetos espaciales
library("sf")
library("sp")

# Trigonometría espacial
library("geosphere")

# Interfaz a GEOS
library("rgeos")

# Geoestadística
library("gstat")
library("geoR")
library("geoRglm")

library("automap")
library("gdalUtils")
library("SpatialPosition")
library("rgdal")

# Algunos estadísticos
library("moments")

# Transformaciones Box-Cox
library("forecast")

# Autocorrelaciones - Indices de Moran y Geary
library("ncf")

# Dependencia espacial
library("spdep")

# Selección de modelos por el método de Akaike
library("MASS")

# Importancia relativa predictores
library("relaimpo")

# Correlogramas
library("corrgram")

# Procesamiento en paralelo y en clusters
library("parallel")

# Gráficos
library("ggplot2")
library("gridExtra")
library("plotrix")
library("RColorBrewer")

# Esquema de colores Monokai 
darkgray <- "#272822"
gray <- "#4D4D4D"
lightgray <- "#75715E"
green <- "#81B023"
orange <- "#C97C16"
purple <- "#8F66CC"
red <- "#C72259"
blue <- "#53A8BD"

# Establecer una conexión ODBC

# Info de la Base de datos
dbName <- "batiValizasDB"
dbHost <- "localhost"
dbPort <- 5432
dbUser <- ""
dbPassword <- ""
dbString <- paste("PG:dbname='", dbName, "' host='", dbHost, "' port='", dbPort, 
    "' user='", dbUser, "' password='", dbPassword, "'", sep = "")

# Cargar datos desde el servidor PostGIS
conn <- dbConnect(dbDriver("PostgreSQL"), dbname = dbName, host = dbHost, port = dbPort, 
    user = dbUser, password = dbPassword)

# Ver tablas en la base de datos
dbListTables(conn)

riverBati <- data.table::fread(input = "tablas/LineaFondoTodo.csv", sep = ",", 
    header = TRUE, showProgress = TRUE, data.table = FALSE, stringsAsFactors = FALSE)

# Sistemas de referencias de coordenadas
epsg.32721 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# Cargar capa espacial del Arroyo Valizas (polígono)
arroyoValizasPoly <- read_sf("capasVectoriales/ArroyoValizasPoly-UTM21.shp")

# Cargar capa espacial de la Línea de Referencia (línea)
refLineSmooth <- read_sf("capasVectoriales/ArroyoValizas-RefLine-Smooth.shp")

# Agregar ID
refLineSmooth$spatial_id <- 1:nrow(refLineSmooth)

# Cambiar el nombre de las columnas a minúsculas
names(refLineSmooth) <- tolower(names(refLineSmooth))

# Plot
plot(st_transform(arroyoValizasPoly, 4326), bgc = darkgray, col = blue, 
    border = blue, lwd = 1, xlab = "Longitud", ylab = "Latitud", graticule = TRUE, 
    col_graticule = lightgray, axes = TRUE, main = "Arroyo Valizas y línea de referencia", 
    cex.axis = 0.75, lat = seq(-34.2, -34.4, by = -0.01))
plot(st_transform(refLineSmooth, 4326), col = red, lwd = 2, add = TRUE)
legend("bottomright", legend = c("Arroyo Valizas", "Línea de referencia"), 
    pch = c(15, NA), lty = c(NA, 1), col = c(blue, red), lwd = c(NA, 3), 
    pt.cex = c(2, 1), bg = "white")

# Remover estados de línea con errores (Line_status = 0)
riverBatiF <- subset(x = riverBati, subset = riverBati$Line_status == 1)

# Crear objeto sf
riverBatiF <- st_as_sf(riverBatiF, coords = c("Longitude", "Latitude"))
st_crs(riverBatiF) <- wgs.84  # assign CRS to points

# Transformar coordenadas
riverBatiF <- st_transform(riverBatiF, crs = epsg.32721)

# Eliminar puntos afuera del polígono del río
batiRiverIntersect <- st_intersects(riverBatiF, arroyoValizasPoly)
riverBatiF2 <- riverBatiF[which(!is.na(batiRiverIntersect)), ]

# Eliminar puntos duplicados
riverBatiF3 <- sp::remove.duplicates(obj = as(riverBatiF2, "Spatial"), 
    remove.second = TRUE)
riverBatiF3 <- st_as_sf(riverBatiF3)

# Add temporary unique ID to spatial DF
riverBatiF3$spatial_id <- 1:nrow(riverBatiF3)

# Set column names to lower case
names(riverBatiF3) <- tolower(names(riverBatiF3))

# Create Regular grid
mygrid5m <- st_make_grid(x = arroyoValizasPoly, cellsize = 5, what = "centers")

# Grid points adentro de polígono de Arroyo
gridRiverIntersect <- st_intersects(mygrid5m, arroyoValizasPoly)
mygrid5m <- mygrid5m[which(!is.na(gridRiverIntersect) & !(gridRiverIntersect == 
    0))]

# Add temporary unique ID to spatial DF
dfAux <- data.frame(spatial_id = 1:length(mygrid5m))
dfAux$geom <- st_geometry(mygrid5m)
mygrid5m <- st_as_sf(dfAux)

# Plot
plot(st_transform(arroyoValizasPoly, 4326), bgc = darkgray, col = blue, 
    border = blue, lwd = 1, xlab = "Longitud", ylab = "Latitud", graticule = TRUE, 
    col_graticule = lightgray, axes = TRUE, main = "Mapa conceptual", cex.axis = 0.75, 
    lat = seq(-34.2, -34.4, by = -5e-04), xlim = c(-53.852, -53.8494), 
    ylim = c(-34.359, -34.358))
plot(st_transform(riverBatiF3, 4326), col = orange, cex = 0.5, add = TRUE)
plot(st_transform(refLineSmooth, 4326), col = red, lwd = 3, add = TRUE)
plot(st_transform(mygrid5m, 4326), col = darkgray, cex = 0.5, pch = 19, 
    add = TRUE)
legend("bottomright", legend = c("Arroyo Valizas", "Batimetría", "Línea de referencia", 
    "Grilla de 5 metros"), pch = c(15, 19, NA, 19), lty = c(NA, NA, 1, NA), 
    col = c(blue, orange, red, darkgray), lwd = c(NA, 3, 3, NA), pt.cex = c(2, 
        0.5, 1, 0.5), bg = "white")

st_write(obj = riverBatiF3, layer = "riverdepthb", driver = "PostgreSQL", 
    dsn = dbString, layer_options = c("geometry_name=geom, OVERWRITE=YES"))
st_write(obj = mygrid5m, layer = "mygrid5mb", driver = "PostgreSQL", dsn = dbString, 
    layer_options = c("geometry_name=geom, OVERWRITE=YES"))
st_write(obj = refLineSmooth, layer = "reflinesmoothb", driver = "PostgreSQL", 
    dsn = dbString, layer_options = c("geometry_name=geom, OVERWRITE=YES"))

# Ver tablas en la base de datos
dbListTables(conn)

# Points Transformed Bati
riverdepthsn <- st_read(dsn = dbString, layer = "riverdepthsn")

# Points Transformed River Grid 5m
rivergridsn <- st_read(dsn = dbString, layer = "rivergridsn")

ggDatosBati <- ggplot(data = st_set_geometry(riverdepthsn, NULL))

# Histograma Profundidad
ggHistBati <- ggDatosBati + geom_histogram(aes(x = depth), binwidth = 0.25, 
    fill = blue, col = darkgray, alpha = 1) + labs(title = "", x = "Profundidad (m)", 
    y = "Frecuencia") + scale_x_continuous(limits = c(1, 7.25), breaks = seq(1, 
    7.25, 0.5)) + scale_y_continuous(limits = c(0, 6000), breaks = seq(0, 
    6000, 500)) + theme(plot.margin = unit(c(-1, 1, 1, 1), "mm"))

# Boxplot profundidad
ggBoxPlotBati <- ggDatosBati + geom_boxplot(aes(y = depth, x = rep("Prof", 
    length(depth))), fill = blue, col = darkgray, alpha = 1) + labs(title = "", 
    x = "", y = "") + scale_y_continuous(limits = c(1, 7.25), breaks = seq(1, 
    7.25, 0.5), position = "top") + coord_flip() + theme(plot.margin = unit(c(1, 
    1, -1, 1), "mm"))

# Ver plot
grid.arrange(ggBoxPlotBati, ggHistBati, ncol = 1, nrow = 2, heights = c(1, 
    3))

ResumenEstadisticas <- function(df) {
    
    if (class(df) == "data.frame") {
        
        # Número de columnas del data.frame
        ncol = ncol(df)
        
    } else {
        
        # Coerce vector to data.frame
        df = data.frame(variable = df)
        ncol = ncol(df)
        
    }
    
    # Estadísticas
    maximo = numeric()
    minimo = numeric()
    
    media = numeric()
    mediana = numeric()
    moda = vector()
    
    varianza = numeric()
    desvioEstandar = numeric()
    cv = numeric()
    asimetria = numeric()
    kurtosis = numeric()
    
    NAs = numeric()
    clases = character()
    
    # Función para calcular la moda
    Moda <- function(x) {
        ux = unique(x)
        ux[which.max(tabulate(match(x, ux)))]
    }
    
    # Iteración que recorre todas las columnas del data.frame
    for (i in 1:ncol) {
        
        # Condición de que la columna para las estadísticas tenga que ser
        # numérica
        if (is.numeric(df[, i])) {
            
            mediana[i] = median(df[, i], na.rm = TRUE)
            media[i] = mean(df[, i], na.rm = TRUE)
            moda[i] = Moda(x = df[, i])
            varianza[i] = var(df[, i], na.rm = TRUE)
            desvioEstandar[i] = sd(df[, i], na.rm = TRUE)
            cv[i] = 100 * (sd(df[, i], na.rm = TRUE)/mean(df[, i], na.rm = TRUE))
            asimetria[i] = moments::skewness(df[, i], na.rm = TRUE)
            kurtosis[i] = moments::kurtosis(df[, i], na.rm = TRUE)
            minimo[i] = min(df[, i], na.rm = TRUE)
            maximo[i] = max(df[, i], na.rm = TRUE)
            NAs[i] = length(which(is.na(df[, i]) == TRUE))
            clases[i] = class(df[, i])
            
        } else {
            # En caso de que no sea de tipo numérica asignar NAs a las 
            # estadísticas pero calcular el número de NAs y tipo de datos
            mediana[i] = NA
            media[i] = NA
            moda[i] = Moda(x = df[, i])
            varianza[i] = NA
            desvioEstandar[i] = NA
            cv[i] = 100 * NA
            asimetria[i] = NA
            kurtosis[i] = NA
            minimo[i] = NA
            maximo[i] = NA
            NAs[i] = length(which(is.na(df[, i]) == TRUE))
            clases[i] = class(df[, i])
            
        }
    }
    
    # Generar data.frame con todas las estadísticas
    summary = data.frame(mínimo = minimo, máximo = maximo, media = media, 
        mediana = mediana, moda = moda, varianza = varianza, ds = desvioEstandar, 
        cv = cv, asimetría = asimetria, kurtosis = kurtosis, NAs = NAs, 
        clase = clases)
    
    # Asignar el nombre de las variables (columnas) como nombre de las
    # filas para el objeto de salida
    rownames(summary) <- colnames(df)
    
    # Retorno de la función
    return(summary)
    
}

# Crear objeto con las principales estadísticas y mostrarlo en consola
(resumenriverdepthsn <- ResumenEstadisticas(df = riverdepthsn$depth))

head(riverdepthsn)

str(riverdepthsn)

# Profundidad
shapiro.test(riverdepthsn$depth[sample(x = 1:nrow(riverdepthsn), size = 5000, replace = FALSE)]) 

# Transformar profundidad

# Logaritmo
ProfLog <- log(riverdepthsn$depth)
shapiro.test(ProfLog[sample(x = 1:nrow(riverdepthsn), size = 5000, replace = FALSE)]) # Probar normalidad de la transformación

# Raíz cuadrada
ProfSqrt <- sqrt(riverdepthsn$depth)
shapiro.test(ProfSqrt[sample(x = 1:nrow(riverdepthsn), size = 5000, replace = FALSE)]) # Probar normalidad de la transformación

# Box-Cox
lambdaProf <- BoxCox.lambda(x = riverdepthsn$depth) # Hallar Lambda
ProfBC <- BoxCox(x = riverdepthsn$depth, lambda = lambdaProf) # Transf.vector
shapiro.test(ProfBC[sample(x = 1:nrow(riverdepthsn), size = 5000, replace = FALSE)]) # Probar normalidad de la transformación 

# Q-Q plot profundidad
ggQQprof <- ggDatosBati + 
  geom_qq(aes(sample = depth)) +
  geom_abline(intercept = mean(riverdepthsn$depth, na.rm = TRUE), 
              slope = sd(riverdepthsn$depth, na.rm = TRUE), col = "red") +
  labs(title = "Profundidad (m)", x = "Cuantiles teóricos", 
       y = "Cuantiles de la muestra")

ggQQprofLog <- ggDatosBati + 
  geom_qq(aes(sample = log(depth))) +
  geom_abline(intercept = mean(log(riverdepthsn$depth), na.rm = TRUE), 
              slope = sd(log(riverdepthsn$depth), na.rm = TRUE), col = "red") + 
  labs(title = "Log Prof.", x = "Cuantiles teóricos", 
       y = "Cuantiles de la muestra")

ggQQprofSqrt <- ggDatosBati + 
  geom_qq(aes(sample = sqrt(depth))) +
  geom_abline(intercept = mean(sqrt(riverdepthsn$depth), na.rm = TRUE), 
              slope = sd(sqrt(riverdepthsn$depth), na.rm = TRUE), col = "red") +
  labs(title = "Raíz cuadrada Prof.", x = "Cuantiles teóricos", 
       y = "Cuantiles de la muestra")

ggQQprofBC <- ggDatosBati + 
  geom_qq(aes(sample = ProfBC)) +
  geom_abline(intercept = mean(ProfBC, na.rm = TRUE), 
              slope = sd(ProfBC, na.rm = TRUE), col = "red") +
  labs(title = "Box-Cox Prof", x = "Cuantiles teóricos", 
       y = "Cuantiles de la muestra")

# Ver plots 
grid.arrange(ggQQprof, ggQQprofLog, ggQQprofSqrt, ggQQprofBC, 
             ncol = 2, nrow = 2)

# Profundidad vs transformación Longitud (N)
ggTend_ProfVsLong <- ggDatosBati +
  geom_point(aes(y = depth, x = riverdepthsn$n)) + 
  geom_smooth(aes(y = depth, x = riverdepthsn$n)) + 
  labs(title = "", x = "n", y = "Profundidad (m)")

# Profundidad vs transformación Latitud (S)
ggTend_ProfVsLat <- ggDatosBati +
  geom_point(aes(y = depth, x = riverdepthsn$s)) + 
  geom_smooth(aes(y = depth, x = riverdepthsn$s)) + 
  labs(title = "", x = "s", y = "Profundidad (m)")

# Ver plots
grid.arrange(ggTend_ProfVsLong, ggTend_ProfVsLat,
             ncol = 1, nrow = 2)

# Profundidad vs transformación Longitud (N)
cor.test(x = riverdepthsn$n, 
         y = riverdepthsn$depth, 
         method = "pearson")

# Profundidad vs transformación Latitud (S)
cor.test(x = riverdepthsn$s, 
         y = riverdepthsn$depth, 
         method = "pearson")

corr <- data.frame(no = 1:20)
corr$gls <- correlog(x = temp.data$x, y = temp.data$y, z = model.gls$residuals, increment = 2, resamp = 100, quiet = T)$correlation
corr$gls.spat <- correlog(x = temp.data$x, y = temp.data$y, z = model.gls.spat$residuals, increment = 2, resamp = 100, quiet = T)$correlation

plot(corr$gls ~ corr$no, pch = 1, xlab = 'Distance', ylab = 'Moran I')
points(corr$gls.spat ~ corr$no, pch = 16)
legend('bottomleft', c('GLS', 'spatial GLS'), pch = c(1,16), cex = 0.7)

# Con librería geoR:

# Objetos geodata
geodataZ <- as.geodata(obj = st_set_geometry(riverdepthsn, NULL), coords.col = 3:4, data.col = 2)

plot(geodataZ)

# Tomar submuestra aleatoria 
geodataZ.sample <- sample.geodata(geodataZ, size = 40000, replace = FALSE)

# Crear variograma empírico
varZ <- variog(geodataZ.sample, max.dist = 500, breaks = seq(from = 0, to = 500, by = 5))
gc()

# Plot variogramas empíricos
plot(varZ, pch = 19, col = blue, ylab = "semivarianza", xlab = "distancia", main = "Variograma empírico de Z")

# gstat data
gstatZ.sample <- gstat(id = "depth", formula = depth ~ 1, data = as(riverdepthsn, "Spatial"))

# Crear variograma empírico omnidireccional
varOmniZ <- variogram(gstatZ.sample, width = 5, cutoff = 500)

# Plotear variograma
plot(varOmniZ, xlab = "distancia", ylab = "semivarianza",
     main = "Variograma empírico - Z", pch = 19, col = blue)

# Compute the 10 degree variogram for the following directions:
var1 <- variogram(gstatZ.sample, alpha = 0, tol.hor = 11.25, cutoff = 300, width = 5)
var2 <- variogram(gstatZ.sample, alpha = 22.5, tol.hor = 11.25, cutoff = 90, width = 2)
var3 <- variogram(gstatZ.sample, alpha = 45, tol.hor = 11.25, cutoff = 45, width = 3)
var4 <- variogram(gstatZ.sample, alpha = 67.5, tol.hor = 11.25, cutoff = 25, width = 2)
var5 <- variogram(gstatZ.sample, alpha = 90, tol.hor = 11.25, cutoff = 22, width = 2)
var6 <- variogram(gstatZ.sample, alpha = 112.5, tol.hor = 11.25, cutoff = 39, width = 3)
var7 <- variogram(gstatZ.sample, alpha = 135, tol.hor = 11.25, cutoff = 35, width = 3)
var8 <- variogram(gstatZ.sample, alpha = 157.5, tol.hor = 11.25, cutoff = 90, width = 2)

# Plotear variograma
par(mfrow = c(2,4))

plot(x = var1$dist, y = var1$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 0 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var1, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var1, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var2$dist, y = var2$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 22.5 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var2, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var2, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var3$dist, y = var3$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 45 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var3, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var3, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var4$dist, y = var4$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 67.5 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var4, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var4, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var5$dist, y = var5$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 90 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var5, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var5, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var6$dist, y = var6$gamma, xlab = "distancia", ylab = "semivarianza", 
     main = "Variograma emp. 112.5 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var6, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var6, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var7$dist, y = var7$gamma, xlab = "distancia", ylab = "semivarianza",
     main = "Variograma emp. 135 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var7, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var7, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

plot(x = var8$dist, y = var8$gamma, xlab = "distancia", ylab = "semivarianza",
     main = "Variograma emp. 157.5 grados", pch = 19, col = blue)
lines(x = variogramLine(fit.variogram(var8, model = vgm("Cir")), maxdist = 400)$dist,
      y = variogramLine(fit.variogram(var8, model = vgm("Cir")), maxdist = 400)$gamma, col = red)

par(mfrow = c(1,1))

# submuestrear submuestra
geodataZ.sample2 <- sample.geodata(geodataZ.sample, 1000, replace = FALSE)

# Matriz de posibles parámetros iniciales de rango y meseta para C1
iniCovParsZ <- cbind("sigma" = seq(0.5, 1.5, length.out = 20),
                     "phi" = seq(10, 200, length.out = 20)
)

iniCovParsZ <- as.matrix(iniCovParsZ) # convertir a matriz

# Ajuste modelo esférico
likfitVarZ.sph <- likfit(geodataZ.sample2,
                         ini.cov.pars = iniCovParsZ,
                         limits = pars.limits(phi = c(10, 200), sigmasq = c(0.5, 1.5)),
                         fix.nugget = FALSE,
                         fix.psiA = FALSE,
                         fix.psiR = FALSE,
                         cov.model = "spherical",
                         lik.method = "REML")

# Ajuste modelo exponencial
likfitVarZ.exp <- likfit(geodataZ.sample2,
                         ini.cov.pars = iniCovParsZ,
                         limits = pars.limits(phi = c(10, 200), sigmasq = c(0.5, 1.5)),
                         fix.nugget = FALSE,
                         fix.psiA = FALSE,
                         fix.psiR = FALSE,
                         cov.model = "exponential",
                         lik.method = "REML")

# Ajuste modelo circular
likfitVarZ.cir <- likfit(geodataZ.sample2,
                         ini.cov.pars = iniCovParsZ,
                         limits = pars.limits(phi = c(10, 200), sigmasq = c(0.5, 1.5)),
                         fix.nugget = FALSE,
                         fix.psiA = FALSE,
                         fix.psiR = FALSE,
                         cov.model = "circular",
                         lik.method = "REML")

# Ajuste modelo Matern
likfitVarZ.mat <- likfit(geodataZ.sample2,
                         ini.cov.pars = iniCovParsZ,
                         limits = pars.limits(phi = c(10, 200), sigmasq = c(0.5, 1.5)),
                         fix.kappa = FALSE,
                         fix.nugget = FALSE,
                         fix.psiA = FALSE,
                         fix.psiR = FALSE,
                         cov.model = "matern",
                         lik.method = "REML")

summary(likfitVarZ.sph)
summary(likfitVarZ.exp)
summary(likfitVarZ.cir)
summary(likfitVarZ.mat)

# Plot variograma empírico + modelos

# Crear variograma empírico
varZ.sample2 <- variog(geodataZ.sample2, max.dist = 200, breaks = seq(from = 0, to = 200, by = 10))

plot(varZ.sample2, main = "Z", ylab = "semivarianza", 
     xlab = "distancia",  pch = 19, col = blue, ylim = c(0,1.5))
lines.variomodel(likfitVarZ.sph, col = "red", lty = 2, max.dist = 200)
lines.variomodel(likfitVarZ.exp, col = "red", lty = 3, max.dist = 200)
lines.variomodel(likfitVarZ.cir, col = "red", lty = 4, max.dist = 200)
lines.variomodel(likfitVarZ.mat, col = "red", lty = 5, max.dist = 200)
grid()
legend('bottomright', lty = c(2,3,4,5), legend = c('Esf','Exp','Cir','Mat'), col = rep("red", 3))

# Selección de modelo por criterio de información de Akaike
AICmodelosZ <- data.frame("Modelo" = c('Esférico', 'Exponencial', 'Circular', 'Matern'),
                          "AIC" = c(likfitVarZ.sph$AIC, likfitVarZ.exp$AIC, likfitVarZ.cir$AIC, likfitVarZ.mat$AIC))

print(AICmodelosZ)

# gstat: anis = c("ángulo eje mayor - 0 grados norte", "eje menor / eje mayor = valor entre 0 y 1"), ámgulo en grados

# geoR: anis = c("ángulo eje mayor - 0 grados norte", "eje mayor / eje menor = mayor o igual a 1"), ángulo en radianes, parece devolver el rango del eje menor y no del mayor(ojo!)

# Radians to degrees function 
rad2deg <- function(rad) {(rad * 180) / (pi)}
anisoAng <- rad2deg(likfitVarZ.cir$aniso.pars[1])
anisoRel <- 1/likfitVarZ.cir$aniso.pars[2]
ejeMenor <- likfitVarZ.cir$phi
ejeMayor <- 1/anisoRel * ejeMenor

varMejorTeoZ <- vgm(psill = likfitVarZ.cir$sigmasq, 
                    model = "Cir", 
                    range = ejeMayor,
                    nugget = likfitVarZ.cir$nugget,
                    anis = c(anisoAng, anisoRel)
)

# Plotear variograma experimental y ajuste modelo teórico
plot(x = varOmniZ$dist, y = varOmniZ$gamma, main = "Mejor ajuste de modelo - Z",
     col = blue, ylab = "semivarianza", xlab = "distancia",
     ylim = c(0, 1.5), xlim = c(0, 200), pch = 19)
lines(variogramLine(object = varMejorTeoZ, maxdist = 200), col = "red", lwd = 1, lty = 1)
legend("bottomright", c("Modelo Circular"), lty = 1, col = "red")

# Clip Spatial Objects
ClipSpObj <- function(xmin, xmax, ymin, ymax, spatialObject) {
  
  vectorClip = which(st_coordinates(spatialObject)[,1] >= xmin & st_coordinates(spatialObject)[,1] < xmax & st_coordinates(spatialObject)[,2] >= ymin & st_coordinates(spatialObject)[,2] < ymax)  
  
  spatialObjectClipped = spatialObject[vectorClip,]
  
  return(spatialObjectClipped)
  
}

# Clip Spatial Objects by Y distance
ClipSpObjByDistY <- function(spatialObject, intervalY) {
  
  listClippedObjs <- list()
  
  xmin <- st_bbox(spatialObject)[1]
  xmax <- st_bbox(spatialObject)[2]
  ymin <- st_bbox(spatialObject)[3]
  ymax <- st_bbox(spatialObject)[4]
  
  intervals <- seq(from = ymin, to = ymax, by = intervalY)
  
  for (i in 1:(length(intervals) - 1)) {
    
    listClippedObjs[[i]] <- ClipSpObj(xmin = xmin, xmax = xmax, ymin = intervals[i], ymax = intervals[i + 1], spatialObject = spatialObject)
    
  }
  
  return(listClippedObjs)
  
}

# Clip depthdata and GridSn Spatial Objects by Y distance and including all the points within Y distance from depthData.  
ClipSpObjByDist <- function(samplingDepthData, predictionGridSn, distY) {
  
  listClippedObjs <- list()
  
  xmin <- st_bbox(samplingDepthData)[1]
  xmax <- st_bbox(samplingDepthData)[2]
  ymin <- st_bbox(samplingDepthData)[3]
  ymax <- st_bbox(samplingDepthData)[4]
  
  intervals <- seq(from = ymin, to = ymax, by = distY)
  
  for (i in 1:(length(intervals) - 1)) {
    
    message(paste("Clipping object", i))
    
    if (i == 1) {
      
      sublist = list()
      
      sublist$"depthDataClipped" = ClipSpObj(xmin = xmin, xmax = xmax, ymin = ymin, ymax = intervals[i + 2], spatialObject = samplingDepthData)
      sublist$"gridSnClipped" = ClipSpObj(xmin = xmin, xmax = xmax, ymin = ymin, ymax = intervals[i + 1], spatialObject = predictionGridSn)
      
      listClippedObjs[[i]] = sublist
      
    } 
    else if (i == (length(intervals) - 1)) {
      
      sublist = list()
      
      sublist$"depthDataClipped" = ClipSpObj(xmin = xmin, xmax = xmax, ymin = intervals[i - 1], ymax = ymax, spatialObject = samplingDepthData)
      
      sublist$"gridSnClipped" = ClipSpObj(xmin = xmin, xmax = xmax, ymin = intervals[i], ymax = ymax, spatialObject = predictionGridSn)
      
      listClippedObjs[[i]] = sublist
      
    }
    else {
      
      sublist = list()
      
      sublist$"depthDataClipped" <- ClipSpObj(xmin = xmin, xmax = xmax, ymin = intervals[i - 1], ymax = intervals[i + 2], spatialObject = samplingDepthData)
      
      sublist$"gridSnClipped" <- ClipSpObj(xmin = xmin, xmax = xmax, ymin = intervals[i], ymax = intervals[i + 1], spatialObject = predictionGridSn)
      
      listClippedObjs[[i]] = sublist
      
    }
  }
  
  return(listClippedObjs)
  
}

# Lista con objetos espaciales cortados
clippedSpObj <- ClipSpObjByDist(samplingDepthData = riverdepthsn, predictionGridSn = rivergridsn, distY = 350)

# 5 meters grid
mygrid5m <- st_read(dsn = dbString, layer = "mygrid5m")
st_crs(mygrid5m) <- 32721

# Machines used
numberOfMachines <- 2

# Vector of hosts, usernames and number of cores
# E.g.: Machine1 = host: localhost, user: guzman and nCores: 7
users <- c('guzman', 'glopez', 'lortega')
hosts <- c('10.5.1.109', '10.5.1.197', '10.5.1.152')
nCores <- c(7, 4, 8)

# Machines's Adresses
machineAddresses <- lapply(X = 1:numberOfMachines, FUN = function(x) list(host = hosts[x], user = users[x], ncore = nCores[x]))

# Specifications
spec <- lapply(machineAddresses, function(machine) {
  rep(list(list(host = machine$host, user = machine$user)), machine$ncore)
})

spec <- unlist(spec, recursive = FALSE)

# Calculate the number of cores
no_cores <- sum(nCores[1:numberOfMachines])

# START CLUSTER
print("Making cluster...")
(parallelCluster <- parallel::makeCluster(type = 'PSOCK', master = hosts[1], spec = spec, homogeneous = FALSE, port = 58178, outfile = ""))
message("done! \n", appendLF = TRUE)

# Write Parallelized function
RunFunctionOnParallelizedCluster <- function(parallelCluster,
                                             varlist,
                                             parts,
                                             functionParallel) {
  message("\n", appendLF = TRUE)
  
  # Load library
  library('parallel')
  
  # Objects to export to clusters
  message("Exporting required objects to clusters...")
  clusterExport(cl = parallelCluster, varlist = varlist)
  message("done! \n", appendLF = TRUE)
  
  # Parallel Cluster Function
  message("Running parallel function...")
  system.time(
    functionOut <- parLapplyLB(
      cl = parallelCluster,
      X = 1:parts,
      fun = function(x)
        functionParallel(x = x)
    )
  )
  message("done! \n", appendLF = TRUE)
  
  # Retorno
  return(functionOut)
  
}

print(parallelCluster)

listaLocations <- lapply(X = 1:length(clippedSpObj), function(x, clippedSpObj) {
  return(clippedSpObj[[x]]$depthDataClipped)
}, clippedSpObj)

listaNewdata <- lapply(X = 1:length(clippedSpObj), function(x, clippedSpObj) {
  return(clippedSpObj[[x]]$gridSnClipped)
}, clippedSpObj)

# Iniciar cluster

# Lista de variables requeridas para ejecutar la función
FunctionParallelKrige <- function(x) {
  
  # Cargar librerías en cluster
  library('gstat')
  library('sf')
  
  # Temporal data frame
  tmpLocations = as(listaLocations[[x]], "Spatial") 
  tmpNewdata = as(listaNewdata[[x]], "Spatial") 
  
  # Function
  output <- krige(formula = depth ~ 1, locations = tmpLocations, newdata = tmpNewdata, model = varMejorTeoZ, maxdist = 350)
  message("Local Kriging finished!")
  
  # Free memory
  gc()
  
  # Retorno
  return(st_as_sf(output))
  
}

# All the objects required to run the function
varlist1 <- c("listaLocations", "listaNewdata", "varMejorTeoZ", "FunctionParallelKrige")

# Ejecutar función en paralelo
output1 <- RunFunctionOnParallelizedCluster(parallelCluster = parallelCluster, varlist = varlist1, parts = length(listaLocations), functionParallel = FunctionParallelKrige)

# Transformar Kriging a coordenadas originales 
TransformInterpolatedData <- function(krigedData, gridSn, grid) {
  
  # Combine kriged data with transformed grid to obtain ID
  message("Binding transformed grid to Kriged object data to get the spatial_id...")
  gridOverDataKriged = cbind(gridSn, krigedData)
  gridOverDataKriged = gridOverDataKriged[,-(ncol(gridOverDataKriged))]
  message("Finished!")
  
  # Combinar capa de datos interpolados con capa grid original sin coordenadas transformadas
  message("Merging Kriged data with original grid by spatial_id...")
  gridOverDataKriged = merge(x = grid, y = st_set_geometry(gridOverDataKriged, NULL), by.x = "spatial_id", by.y = "spatial_id")
  message("Finished!")
  
  # Retorno
  message("Transformation finished!")
  return(gridOverDataKriged)
  
}

# Lista de variables requeridas para ejecutar la función
FunctionParallelTransformInterpolatedData <- function(x) {
  
  # Cargar librerías en cluster
  library('sf')
  
  # Function
  message(paste("Transforming interpolated data"))
  output <- TransformInterpolatedData(krigedData = output1[[x]], gridSn = listaNewdata[[x]], grid = mygrid5m)
  
  # Retorno
  return(output)
  
}

# All the objects requiered to run the function
varlist2 = c("output1", "listaNewdata", "mygrid5m", "TransformInterpolatedData", "FunctionParallelTransformInterpolatedData")

# Ejecutar función en paralelo
output2 <- RunFunctionOnParallelizedCluster(parallelCluster = parallelCluster, varlist = varlist2, parts = length(output1), functionParallel = FunctionParallelTransformInterpolatedData)

# Convertir grilla regular de puntos a polígonos
TransformRegularPointsToPolygons <- function(regularPoints, gridResolution) {
  
  # Convert points centroids to polygon
  # set the radius (meters) for the plots
  message("Setting radius...")
  radius = gridResolution / 2
  
  # define the plot boundaries based upon the plot radius
  message("Defining plot polygons boundaries...")
  yPlus = st_coordinates(regularPoints)[,2] + radius
  xPlus = st_coordinates(regularPoints)[,1] + radius
  yMinus = st_coordinates(regularPoints)[,2] - radius
  xMinus = st_coordinates(regularPoints)[,1] - radius
  
  # Calculate polygon coordinates for each plot centroid 
  message("Polygon coordinates for vertices...")
  square = cbind(xMinus, yMinus, xMinus, yPlus, xPlus, yPlus, xPlus, yMinus, xMinus, yMinus)
  
  # Create List of POLYGONS
  message("Creating list of POLYGONS...")
  listPolys = lapply(X = 1:nrow(square), FUN = function(x) {
    st_polygon(list(do.call("rbind", split(square[x,], cut(1:10, breaks = 5)))))
  }
  )
  
  message("Creating Simple Feature Collection (sfc)...")
  multiPolygon = st_sf(ID = as.character(1:nrow(square)), st_sfc(listPolys))
  
  # Retorno  
  message("Finished!")
  return(multiPolygon)
  
}

# Lista de variables requeridas para ejecutar la función
FunctionParallelTransformRegularPointsToPolygons <- function(x) {
  
  # Cargar librerías en cluster
  library('sf')
  
  # Function
  message(paste("Transforming regular points to polygons..."))
  output <- TransformRegularPointsToPolygons(regularPoints = output2[[x]], gridResolution = 5)
  st_crs(output) <- st_crs(32721)
  
  # Retorno
  return(output)
  
}

# All the objects requiered to run the function
varlist3 = c("output2", "TransformRegularPointsToPolygons", "FunctionParallelTransformRegularPointsToPolygons")

# Ejecutar función en paralelo
output3 <- RunFunctionOnParallelizedCluster(parallelCluster = parallelCluster, varlist = varlist3, parts = length(output2), functionParallel = FunctionParallelTransformRegularPointsToPolygons)

# Combinar polígonos de la grilla regular y Kriging
MergePolygonsWithKrigedData <- function(krigedData, polygons) {
  
  # Misma proyección
  message("Check that 'krigedData' and 'polygons' have same CRS...")
  st_crs(krigedData) <- epsg.32721
  st_crs(polygons) <- epsg.32721
  
  # Combinar kriged data con polys
  message("Intersect 'krigedData' with 'polygons'...")
  polysOverDataKriged = st_intersects(x = polygons, y = krigedData)
  noIntersect <- which(is.na(polysOverDataKriged) | polysOverDataKriged == 0)
  
  if (length(noIntersect) > 0) {
    
    warning("Not intersected elements found!")
    
  }
  
  polysOverDataKriged <- polygons[-noIntersect]
  
  # Recombinar
  message("Creating new 'sf' object...")
  polysOverDataKriged2 <- st_sf(data.frame(krigedData)[,-ncol(krigedData)], st_geometry(polysOverDataKriged))
  
  message("Finished!")
  return(polysOverDataKriged2)
  
}

# Lista de variables requeridas para ejecutar la función
FunctionParallelMergePolygonsWithKrigedData <- function(x) {
  
  # Cargar librerías en cluster
  library('sf')
  
  # Function
  message(paste("Bind Polygons with predicted data..."))
  output <- MergePolygonsWithKrigedData(krigedData = output2[[x]], polygons = output3[[x]])
  
  # Retorno
  return(output)
  
}

# All the objects requiered to run the function
varlist4 = c("epsg.32721", "output2", "output3", "MergePolygonsWithKrigedData", "FunctionParallelMergePolygonsWithKrigedData")

# Ejecutar función en paralelo
output4 <- RunFunctionOnParallelizedCluster(parallelCluster = parallelCluster, varlist = varlist4, parts = length(output3), functionParallel = FunctionParallelMergePolygonsWithKrigedData)

# Merge all the predictions
filter <- -9999

mergeParts <- rbind(output4[[1]], output4[[2]])

for (j in 3:length(output4)) {
  
  if (j != filter) {
    mergeParts <- rbind(mergeParts, output4[[j]])
    message(paste("Simple feature collection ", j, " merged!", sep = ""))
  }
}

# Convert to SpatialPixelsDataFrame
mergePartsPixels <- sp::SpatialPixelsDataFrame(points = as(mergeParts, "Spatial"), data = as(mergeParts, "Spatial")@data)
mergePartsPixels$var1.pred <- -1 * mergePartsPixels$var1.pred

# Convert to raster
mergePartsRaster <- raster::raster(mergePartsPixels[,'var1.pred'])

# Export raster
raster::writeRaster(x = mergePartsRaster, filename = "/home/guzman/Documentos/DINARA/SIG/Arroyo Valizas/BathymetryModel/capasVectoriales/raster-bati-arroyovalizas-01.tif", format = "GTiff", overwrite = TRUE)

# Export polygons
write_sf(obj = mergeParts, dsn = "/home/guzman/Documentos/DINARA/SIG/Arroyo Valizas/BathymetryModel/capasVectoriales/", layer = "polys-bati-arroyovalizas-01", driver = "ESRI Shapefile", overwrite_layer = TRUE)

# Plot raster
raster::plot(mergePartsRaster, colNA = 'black', useRaster = TRUE, interpolate = TRUE, col = oce::oceColors9A(n = 128), alpha = 1)

# Cargar márgenes del Arroyo Valizas
ArroyoValizasLineasParalelas <- st_read(dsn = 'capasVectoriales/', layer = 'ArroyoValizasLineasParalelas')
ArroyoValizasLineasParalelas.32721 <- st_transform(ArroyoValizasLineasParalelas, crs = 32721)

# Plot
xmin <- min(st_bbox(output4[[6]])[1], 
            st_bbox(output4[[6]])[1])
ymin <- min(st_bbox(output4[[6]])[2], 
            st_bbox(output4[[6]])[2])
xmax <- max(st_bbox(output4[[6]])[3], 
            st_bbox(output4[[6]])[3])
ymax <- max(st_bbox(output4[[6]])[4], 
            st_bbox(output4[[6]])[4])

# Plot
plot(arroyoValizasPoly, lwd = 0.01, col = gray,
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Márgenes del Arroyo Valizas a partir de digitalización en QGIS")
plot(ArroyoValizasLineasParalelas.32721, col = c(orange, green), add = TRUE, lwd = 3)
text(x = rep(0.5 * (xmin + xmax), 2), 
     y = c(0.5 * (ymin + ymax) + 25, 0.5 * (ymin + ymax) - 30), 
     c("Margen superior", "Margen inferior"), pos = c(3,1), 
     col = c(orange, green), cex = 1.2)

# Función para encontrar puntos a una determinada distancia fija sobre una línea
find <- function(line, distance) {
  
  norm_vec <- function(x) sqrt(sum(x^2))
  
  new_point <- function(p0, p1, di) {
    
    # Finds point in distance di from point p0 in direction of point p1
    v = p1 - p0
    u = v / norm_vec(v)
    
    return(p0 + u * di)
  }
  
  result = line[1, , drop = FALSE]
  
  # for all subsequent points p1, p2 in this data.frame norm_vec(p2 - p1) = distance at all times
  equidistantPoints = line[1, , drop = FALSE]
  line = tail(line, n = -1)
  accDist = 0
  
  
  while (nrow(line) > 0) {
    point = line[1,]
    lastPoint = result[1,]
    
    dist = norm_vec(point - lastPoint)
    
    if ( accDist + dist > distance ) {
      np = new_point(lastPoint, point, distance - accDist)
      equidistantPoints = rbind(np, equidistantPoints) # add np to equidistantPoints
      result = rbind(np, result) # add np to result
      accDist = 0 # reset accDist
    } else {
      # move point from line to result
      line = tail(line, n = -1)
      result = rbind(point, result)
      # update accDist
      accDist = accDist + dist
    }
  }
  allPoints = result[NROW(result):1,] # reverse result
  return(list(newPoints = equidistantPoints, allPoints = allPoints))
}

# Puntos equidistantes para MargenSup (elimina puntos de la línea, simplifica)
margenSup <- ArroyoValizasLineasParalelas[ArroyoValizasLineasParalelas$id == 1,]
margenSup <- st_transform(margenSup, crs = 32721)
margenSup <- find(line = st_coordinates(margenSup)[,c(1,2)], distance = 10)
margenSup <- st_sf(id = 'margenSup', st_sfc(st_linestring(margenSup$newPoints, dim = "XY")))
st_crs(margenSup) <- st_crs(32721)
margenSup <- st_transform(margenSup, 4326)

# Puntos equidistantes para MargenInf (elimina puntos de la línea, simplifica)
margenInf <- ArroyoValizasLineasParalelas[ArroyoValizasLineasParalelas$id == 2,]
margenInf <- st_transform(margenInf, crs = 32721)
margenInf <- find(line = st_coordinates(margenInf)[,c(1,2)], distance = 10)
margenInf <- st_sf(id = 'margenInf', st_sfc(st_linestring(margenInf$newPoints, dim = "XY")))
st_crs(margenInf) <- st_crs(32721)
margenInf <- st_transform(margenInf, 4326)

# Plot
plot(arroyoValizasPoly, lwd = 0.1, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Márgenes conformados con puntos cada 10 metros")
plot(st_transform(margenSup, crs = 32721), add = TRUE, col = orange, lwd = 2)
plot(st_transform(margenInf, crs = 32721), add = TRUE, col = green, lwd = 2)
plot(st_cast(st_transform(margenSup, crs = 32721), "POINT"), add = TRUE, col = orange, pch = 19)
plot(st_cast(st_transform(margenInf, crs = 32721), "POINT"), add = TRUE, col = green, pch = 19)

# Convertir márgenes nuevos a polígonos
polygonizeNewMargins <-
  st_sf(id = "Arroyo_Valizas", st_sfc(st_polygon(list(
    rbind(
      st_coordinates(margenSup)[, c(1, 2)],
      st_coordinates(margenInf)[, c(1, 2)],
      st_coordinates(margenSup)[1, c(1, 2)]
    )
  ))))
st_crs(polygonizeNewMargins) <- st_crs(4326)
polygonizeNewMargins <- st_transform(polygonizeNewMargins, crs = 32721)

# Distancia entre puntos del margen superior e inferior

# Generar puntos a partir de líneas de márgenes
pointsMarSup = data.frame("lon" = rev(st_coordinates(margenSup)[,1]), 
                          "lat" = rev(st_coordinates(margenSup)[,2]))

pointsMarSup <- st_as_sf(pointsMarSup, coords = c("lon", "lat"))
st_crs(pointsMarSup) <- st_crs(4326)
st_transform(pointsMarSup, crs = 32721)

pointsMarInf = data.frame("lon" = st_coordinates(margenInf)[,1], 
                          "lat" = st_coordinates(margenInf)[,2])

pointsMarInf <- st_as_sf(pointsMarInf, coords = c("lon", "lat"))
st_crs(pointsMarInf) <- st_crs(4326)
st_transform(pointsMarInf, crs = 32721) 

# Calcular distancias
distBetweenPoints <- lapply(X = 1:nrow(pointsMarSup), FUN = function(x) st_distance(x = pointsMarSup[x,], y = pointsMarInf[x,]))

# Juntar
distBetweenPoints <- do.call('rbind', distBetweenPoints) # unir partes
distBetweenPoints <- data.frame("lon1" = st_coordinates(pointsMarSup)[,1], 
                                "lat1" = st_coordinates(pointsMarSup)[,2], 
                                "lon2" = st_coordinates(pointsMarInf)[,1], 
                                "lat2" = st_coordinates(pointsMarInf)[,2], 
                                "dist" = distBetweenPoints)

print(head(distBetweenPoints))

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Distancias entre márgenes")
plot(st_transform(pointsMarSup, crs = 32721), add = TRUE, col = orange, lwd = 2, pch = 19)
plot(st_transform(pointsMarInf, crs = 32721), add = TRUE, col = green, lwd = 2, pch = 19)

listDistMeasured <- list()

for (i in 1:nrow(distBetweenPoints)) {
  
  p1 = st_as_sf(x = distBetweenPoints[i, c(1,2)], coords = c("lon1", "lat1"))
  st_crs(p1) =  st_crs(4326) 
  p2 = st_as_sf(x = distBetweenPoints[i, c(3,4)], coords = c("lon2", "lat2"))
  st_crs(p2) =  st_crs(4326) 
  
  listDistMeasuredTemp = geosphere::gcIntermediate(p1 = st_coordinates(p1), 
                                                   p2 = st_coordinates(p2), 
                                                   n = 10)
  
  listDistMeasuredTemp <- st_sf(id = 'dist',
                                st_sfc(st_linestring(listDistMeasuredTemp, dim = "XY")))
  st_crs(listDistMeasuredTemp) <- st_crs(4326)
  
  listDistMeasured[[i]] <- st_transform(listDistMeasuredTemp, crs = 32721)
  plot(st_transform(listDistMeasuredTemp, crs = 32721), col = red, lty = 3, lwd = 3, add = TRUE)
  
}

# Puntos en el medio de la recta entre punto de la línea 1 y punto de la línea 2 
midPoints <- geosphere::midPoint(p1 = st_coordinates(pointsMarSup), 
                                 p2 = st_coordinates(pointsMarInf))
midPointsSF <- st_as_sf(as.data.frame(midPoints), coords = c("lon", "lat"))

# Unir todos los puntos medios
midLine <- st_sf(id = 'midline', st_sfc(st_linestring(midPoints, dim = "XY")))
st_crs(midLine) <- st_crs(4326)


# Splining a polygon.
#
#   The rows of 'xy' give coordinates of the boundary vertices, in order.
#   'vertices' is the number of spline vertices to create.
#              (Not all are used: some are clipped from the ends.)
#   'k' is the number of points to wrap around the ends to obtain
#       a smooth periodic spline.
#
#   Returns an array of points. 
# 
spline.line <- function(xy, vertices, k = 3, ...) {
  # Assert: xy is an n by 2 matrix with n >= k.
  
  # Wrap k vertices around each end.
  n <- dim(xy)[1]
  if (k >= 1) {
    data <- rbind(xy[(n - k + 1):n,], xy, xy[1:k, ])
  } else {
    data <- xy
  }
  
  # Spline the x and y coordinates.
  data.spline <- spline(1:(n + 2 * k), data[,1], n = vertices, ...)
  x <- data.spline$x
  x1 <- data.spline$y
  x2 <- spline(1:(n + 2 * k), data[,2], n = vertices, ...)$y
  
  # Retain only the middle part.
  cbind(x1, x2)[k < x & x <= n + k, ]
}

midLineSpline <- spline.line(st_coordinates(midLine), vertices = 1000, k = 3)
midLineSpline <- st_sf(id = 'midline_spline', st_sfc(st_linestring(midLineSpline, dim = "XY")), crs = 4326)

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Distancias entre márgenes")
plot(st_transform(pointsMarSup, crs = 32721), add = TRUE, col = orange, lwd = 2, pch = 19)
plot(st_transform(pointsMarInf, crs = 32721), add = TRUE, col = green, lwd = 2, pch = 19)
plot(st_transform(midLine, crs = 32721), add = TRUE, col = red, lwd = 1)
plot(st_transform(midLineSpline, crs = 32721), add = TRUE, col = purple, lwd = 1)

# Escribir midLine a un shapefile
# write_sf(obj = midLine, dsn = "../capasVectoriales/", layer = "midLine", driver = "ESRI Shapefile", overwrite_layer = TRUE)

# Edit centerline in QGIS and smooth it using the plugin "Generalizer" (Algorithm: Ckaiken, level: 2, weight: 3)

# Cargar línea media 
centerLine <- read_sf(dsn = "capasVectoriales/", layer = "midLine")
centerLineUTM <- st_transform(centerLine, 32721)

# Puntos equidistantes sobre línea media 
equalSpacedPoints <- find(line = st_coordinates(centerLineUTM)[,c(1,2)], distance = 10)

# Convertir a LINESTRING
equalSpacedCenterLine <- st_sf(id = 'equalspacedpoints', st_sfc(st_linestring(equalSpacedPoints$newPoints, dim = "XY")))
st_crs(equalSpacedCenterLine) <- st_crs(32721)

# Convertir a PUNTOS
equalSpacedCenterPoints <- st_sf(id = 'equalspacedpoints', st_sfc(st_multipoint(equalSpacedPoints$newPoints, dim = "XY")))
st_crs(equalSpacedCenterPoints) <- st_crs(32721)
equalSpacedCenterPoints <- st_cast(equalSpacedCenterPoints, "POINT")
equalSpacedCenterPointsWGS84 <- st_transform(equalSpacedCenterPoints, 4326)

# Líneas perpendiculares desde línea media

# Puntos centrales entre dos puntos de la línea media 
midEqualSpacedCenterPointsWGS84 <- geosphere::midPoint(p1 = st_coordinates(equalSpacedCenterPointsWGS84)[-nrow(equalSpacedCenterPointsWGS84),c(1,2)], 
                                                       p2 = st_coordinates(equalSpacedCenterPointsWGS84)[-1,c(1,2)])

midEqualSpacedCenterPoints <- st_sf(id = 'midEqualSpacedCenterPoints', st_sfc(st_multipoint(midEqualSpacedCenterPointsWGS84, dim = "XY")))
midEqualSpacedCenterPoints <- st_cast(midEqualSpacedCenterPoints, "POINT")
st_crs(midEqualSpacedCenterPoints) <- st_crs(4326)
midEqualSpacedCenterPoints <- st_transform(midEqualSpacedCenterPoints, 32721)

# Segmentos de línea
segmentsList = lapply(X = 1:(nrow(equalSpacedCenterPoints) - 1), FUN = function(x) {
  
  pts <- st_coordinates(equalSpacedCenterPoints)[x:(x + 1),c(1,2)]
  linestring <- st_linestring(pts, dim = "XY")
  return(linestring)
  
})

segmentsList <- st_sfc(segmentsList)
st_crs(segmentsList) <- st_crs(32721)

# Función para generar líneas perpendiculares de cada segmento de línea
PerpendicularLinesFromSegments <- function(linestring, pointsInLines, distance) {
  
  n = length(linestring)
  
  perpLinesList = list()
  
  for (i in 1:n) {
    
    x1 = st_coordinates(linestring[i])[,c(1,2)][1]
    y1 = st_coordinates(linestring[i])[,c(1,2)][3]
    
    x2 = st_coordinates(linestring[i])[,c(1,2)][2]
    y2 = st_coordinates(linestring[i])[,c(1,2)][4]
    
    x3 = st_coordinates(pointsInLines)[i,c(1,2)][1]
    y3 = st_coordinates(pointsInLines)[i,c(1,2)][2]
    
    vx = x2 - x1
    vy = y2 - y1
    len = sqrt( vx*vx + vy*vy )
    ux = -vy/len
    uy = vx/len
    
    x4 = x3 + distance * ux
    y4 = y3 + distance * uy
    
    x5 = x3 - distance * ux
    y5 = y3 - distance * uy
    
    m = matrix(c(x4,x5,y4,y5), 2, 2)
    
    perpLinesList[[i]] = st_linestring(m, dim = "XY")
    
  }
  
  perpLinesListSL = st_sfc(perpLinesList)
  st_crs(perpLinesListSL) <- st_crs(32721)
  
  # Retorno
  return(perpLinesListSL)
  
}

# Líneas perpendiculares desde cada segmento de la línea media 
perpLines <- PerpendicularLinesFromSegments(linestring = segmentsList, pointsInLines = midEqualSpacedCenterPoints, distance = 100)

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Líneas perpendiculares a línea media")
plot(equalSpacedCenterLine, add = TRUE, col = purple, lwd = 2, pch = 19)
plot(segmentsList[1490], add = TRUE, pch = 19, col = red, lwd = 4)
plot(midEqualSpacedCenterPoints, add = TRUE, col = purple, pch = 'x')
plot(perpLines[1490], add = TRUE, col = red, pch = 19, lwd = 2)

# Sample interpolated polygon data layer
clusterExport(parallelCluster, varlist = c('perpLines', 'find'))
perpLinesPoints5m <- parLapply(cl = parallelCluster, X = 1:length(perpLines), fun = function(x) find(line = as.data.frame(sf::st_coordinates(perpLines[x]))[,c(1,2)], distance = 5))

perpLinesPoints5mList <- list()

for (i in 1:length(perpLines)) {
  
  perpLinesPoints5mListTemp <- st_sf(id = paste(i), st_sfc(st_linestring(as.matrix(perpLinesPoints5m[[i]]$newPoints), dim = "XY")))
  st_crs(perpLinesPoints5mListTemp) <- st_crs(32721)
  perpLinesPoints5mList[[i]] <- perpLinesPoints5mListTemp
  
}

perpLinesPoints5mList <- do.call("rbind", perpLinesPoints5mList)

# Convertir líneas perpendiculares a puntos y seleccionar los puntos dentro del Arroyo y para cada parte 
clusterExport(parallelCluster, varlist = c('output3','perpLinesPoints5mList'))
perpLinesParts <- parLapply(cl = parallelCluster, X = 1:length(output3), fun = function(x) sf::st_intersection(output3[[x]], st_cast(perpLinesPoints5mList, "POINT")))

# Extraer datos de los polígonos interpolados 
ExtractDataFromInterpolatedPolys <- function(perpLinesPoints, intDataPoly) {
  
  # Intersección entre línea perpendicular y polígonos interpolados
  intersectPolyPerpLine = st_intersection(x = perpLinesPoints, y = intDataPoly)
  
  # Retorno
  return(intersectPolyPerpLine)
  
}

# Iniciar cluster

# Función para Paralelizar:
FunctionParallelExtractDataFromInterpolatedPolys <- function(x) {
  
  # Cargar librerías en cluster
  library('sf')
  
  # Function
  message("Extracting data from interpolated polygons...")
  output <- ExtractDataFromInterpolatedPolys(perpLinesPoints = perpLinesParts[[x]], intDataPoly = output4[[x]])
  message("Finished!")
  
  # Retorno
  return(output)
  
}

# Lista de variables requeridas para ejecutar la función
varlist6 <- c("perpLinesParts", "output4", "ExtractDataFromInterpolatedPolys", "FunctionParallelExtractDataFromInterpolatedPolys")

# Ejecutar función en paralelo
output6 <- RunFunctionOnParallelizedCluster(parallelCluster = parallelCluster, varlist = varlist6, parts = length(perpLinesParts), functionParallel = FunctionParallelExtractDataFromInterpolatedPolys)

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Extraer datos de profundidad para cada serie de puntos")
plot(output6["id"], add = TRUE, pch = 19)

# Juntar todos los puntos muestreados
output6 <- do.call("rbind", output6)

# Export to shapefile to edit in QGIS
#st_write(obj = output6, dsn = "capasVectoriales/thalweg/", layer = "perp-lines-points-extract", driver = "ESRI Shapefile", overwrite_layer = TRUE)

# Ver perfil transversal
profiles <- split(x = output6, f = output6$id, drop = TRUE) 
profilesList <- lapply(seq_along(profiles), function(x) {
  
  p <- profiles[[x]]  
  distPointInfMargin <- st_distance(st_transform(margenInf, 32721), p)
  p <- p[order(distPointInfMargin),]
  
  p$along.dist <- st_distance(p, p)[,1]
  
  return(p)
  
})

profiles <- do.call("rbind", profilesList)

# Plot (# ver acá de hacer un plot con tres perfiles a la derecha tipo polígono)
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Perfiles transversales")
plot(profiles[profiles$id == 1490,][,9], add = TRUE, pch = 19)
text(x = st_coordinates(profiles[profiles$id == 1490,][,9])[,1],
     y = st_coordinates(profiles[profiles$id == 1490,][,9])[,2],
     labels = round(x = -(profiles[profiles$id == 1490,][,9]$var1.pred), digits = 2), 
     pos = 4)
plot(profiles[profiles$id == 1480,][,9], add = TRUE, pch = 19)
text(x = st_coordinates(profiles[profiles$id == 1480,][,9])[,1],
     y = st_coordinates(profiles[profiles$id == 1480,][,9])[,2],
     labels = round(x = -(profiles[profiles$id == 1480,][,9]$var1.pred), digits = 2), 
     pos = 4)

# Plot
plot(x = profiles[profiles$id == 1480,]$along.dist,
     y = -profiles[profiles$id == 1480,]$var1.pred, 
     xlab = "Distancia (metros)", ylab = "Profundidad (metros)",
     pch = 19, type = "b")

# gg <- ggplot(data.frame("intervalos" = moranRendMaizC2$mean.of.class,
#                         "correlación" = moranRendMaizC2$correlation)) +
#   geom_polygon(aes(x = polyAutocorrPos$x, y = polyAutocorrPos$y, 
#                    group = polyAutocorrPos$id), 
#                    fill = "darkblue", alpha = 0.6) 

# Encontrar los puntos más profundos en cada perfil
riverChannelPointsList <- lapply(seq_along(profilesList), function(x) {
  
  formatExtractedDepth <- profilesList[[x]]
  
  # Buscar el punto más profundo en cada perfil
  minDepthPoint <- formatExtractedDepth[which(-formatExtractedDepth$var1.pred == min(-formatExtractedDepth$var1.pred)),]
  
  # Retornar punto más profundo
  return(minDepthPoint)
  
}) 

# Puntos del canal 
riverChannelPointsList = do.call("rbind", riverChannelPointsList)

# Export to shapefile to edit in QGIS
#st_write(obj = riverChannelPointsList, dsn = "capasVectoriales/thalweg/", layer = "thalweg-points", driver = "ESRI Shapefile", overwrite_layer = TRUE)

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Puntos probables del canal")
plot(output4[[6]][,7], add = TRUE)
plot(riverChannelPointsList, add = TRUE, col = "white", pch = 19, cex = 1)
plot(riverChannelPointsList, add = TRUE, col = "black", pch = 19, cex = 0.75)

# Cargar línea del canal editada en QGIS
mergePartsSMP_edited <- st_read(dsn = 'capasVectoriales/thalweg/', layer = 'thalweg-points-edited-qgis')

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Puntos probables del canal")
plot(output4[[6]][,7], add = TRUE)
plot(mergePartsSMP_edited, add = TRUE, col = "white", pch = 19, cex = 1)
plot(mergePartsSMP_edited, add = TRUE, col = "black", pch = 19, cex = 0.75)

MakeLinestringFromUnsortedPoints <- function(sfObjectPoint, startPointCell) {
  
  # Check if sfObject is a SpatialPoints object
  isSP = class(sfObjectPoint)[1] == "sf"
  
  if (!isSP) {
    
    return(warning("The sf object argument is not a SimpleFeatures object"))
    
  }
  
  # Remove duplicate points
  
  # Eliminar puntos duplicados
  n = nrow(sfObjectPoint)
  
  sfObjectPoint <- sp::remove.duplicates(obj = as(sfObjectPoint, "Spatial"), remove.second = TRUE)
  sfObjectPoint <- st_as_sf(sfObjectPoint)
  
  # Warn if there are duplicate points in data
  if (n > nrow(sfObjectPoint)) {
    
    warning("Duplicate points in your SimpleFeatures object. The duplicated points were removed.")
    
  }
  
  seqCell = integer()
  seqCell[1] = startPointCell
  
  # Now calculate pairwise distances between points with replacement
  for (i in 2:nrow(sfObjectPoint)) {
    
    distancesFromPoint = data.frame('id' = 1:nrow(sfObjectPoint),
                                    'dist' = st_distance(x = sfObjectPoint[seqCell[i - 1],],
                                                         y = sfObjectPoint)[1,])
    
    distancesFromPoint[seqCell,"dist"] <- NA
    distances = distancesFromPoint$dist
    minDist = min(distances, na.rm = TRUE)
    minDistCell = which(distances == minDist)
    seqCell[i] = minDistCell
    
  }
  
  # LINESTRING
  linestring = st_sfc(st_linestring(st_coordinates(sfObjectPoint[seqCell,])), crs = 32721)
  
  
  return(linestring)
  
}

# Convertir a linestring
startPointCell <- which(st_coordinates(mergePartsSMP_edited)[,1] == min(st_coordinates(mergePartsSMP_edited)[,1]))

thalwegLine <- MakeLinestringFromUnsortedPoints(sfObjectPoint = mergePartsSMP_edited, startPointCell = startPointCell)

# Plot
plot(thalwegLine, col = gray,
     xlab = "Longitude", ylab = "Latitude",
     graticule = st_crs(4326), axes = TRUE, 
     main = "Canal del Arroyo", lwd = 2)

# Plot
plot(polygonizeNewMargins, col = gray, 
     xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
     graticule = st_crs(4326), axes = TRUE, main = "Canal del Arroyo")
plot(thalwegLine, col = blue, add = TRUE, lwd = 2)

# Escribir canal a un archivo 
# st_write(obj = thalwegLine, dsn = 'capasVectoriales/thalweg/', layer = "thalweg-line", driver = "ESRI Shapefile", overwrite_layer = TRUE)

# Escribir canal del arroyo a Base de Datos PostGIS
# st_write(obj = thalwegLine, layer = "thalwegline", driver = "PostgreSQL", dsn = dbString, layer_options = "geometry_name=geom", overwrite_layer = TRUE)

# Encontrar n

# Cargar datos de batimetría filtrados desde base de datos
riverdepth <- st_read(dsn = dbString, layer = "riverdepth")

# Punto más cercao a canal de datos de batimetría

# Lista de variables requeridas para ejecutar la función
FunctionParallelSnapPointsToLines <- function(x) {
  
  # Cargar librerías en cluster
  library('maptools')
  library('sf')
  
  riverdepth <- as(riverdepth, "Spatial")
  thalwegLine <- as(thalwegLine, "SpatialLines")
  
  # Parts
  partsRiverdepth <- split(x = 1:length(riverdepth), f = cut(x = 1:length(riverdepth), no_cores))
  
  # Function
  message("Juntando puntos a línea...")
  output <- snapPointsToLines(points = riverdepth[partsRiverdepth[[x]],], lines = thalwegLine, withAttrs = TRUE)
  
  # Retorno
  return(output)
  
}

# All the objects requiered to run the function
varlist5 <- c("riverdepth",
              "thalwegLine",
              "FunctionParallelSnapPointsToLines",
              "no_cores")

# Ejecutar función en paralelo
output5 <- RunFunctionOnParallelizedCluster(
  parallelCluster = parallelCluster,
  varlist = varlist5,
  parts = no_cores,
  functionParallel = FunctionParallelSnapPointsToLines
)

# Combinar
output5merged <- rbind(output5[[1]], output5[[2]])

for (j in 3:length(output5)) {
  
  output5merged <- rbind(output5merged, output5[[j]])
  message(paste("Simple feature collection ", j, " merged!", sep = ""))
  
}

# Convertir a sf
output5merged <- st_as_sf(output5merged)

# Distancia mínima desde línea del canal a los puntos de batimetría
n <- numeric()

for (i in 1:nrow(riverdepth)) {
  
  n[i] <- st_distance(riverdepth[i,], output5merged[i,])
  message(paste("Distance ", i, sep = ""))
  
}

# Seleccionar puntos a un lado y al otro de la línea del canal

# Agregar margen
offsetX <- 0
offsetY <- 100

coordsSide1 <- rbind(c(st_bbox(thalwegLine)['xmax'] + offsetX, 
                       st_bbox(thalwegLine)['ymax'] + offsetY),
                     c(st_bbox(thalwegLine)['xmin'] - offsetX, 
                       st_bbox(thalwegLine)['ymax'] + offsetY),
                     as.data.frame(st_coordinates(thalwegLine))[,c(1,2)],
                     c(st_bbox(thalwegLine)['xmax'] + offsetX, 
                       st_bbox(thalwegLine)['ymax'] + offsetY))

coordsSide2 <- rbind(c(st_bbox(thalwegLine)['xmax'] + offsetX, 
                       st_bbox(thalwegLine)['ymin'] - offsetY),
                     c(st_bbox(thalwegLine)['xmin'] - offsetX, 
                       st_bbox(thalwegLine)['ymin'] - offsetY),
                     as.data.frame(st_coordinates(thalwegLine))[,c(1,2)],
                     c(st_bbox(thalwegLine)['xmax'] + offsetX, 
                       st_bbox(thalwegLine)['ymin'] - offsetY))

# Convertir coordenadas a Polígonos 
polySide1 <- st_sf("id" = 'side1', st_sfc(st_polygon(list(as.matrix(coordsSide1))), crs = 32721))
polySide2 <- st_sf("id" = 'side2', st_sfc(st_polygon(list(as.matrix(coordsSide2))), crs = 32721))

# Seleccionar puntos en el Lado 1
pointsInSide1 <- st_intersection(riverdepth, polySide1)

# Seleccionar puntos en el Lado 2
pointsInSide2 <- st_intersection(riverdepth, polySide2)

# Seleccionar puntos sobre la línea del canal
pointsOverLine <- st_intersects(riverdepth, thalwegLine)
sampleSpatialPointsOverLine <- riverdepth[pointsOverLine, ] # NULL

# Valores negativos
negCells <- which(st_intersects(riverdepth, polySide2, sparse = FALSE))
n[negCells] <- n[negCells] * -1

# Add n coordinates attribute to depth points
riverdepth$n <- n