#Contexto

Este proyecto tiene como objetivo desarrollar un modelo de interpolación espacial mediante Kriging Ordinario para caracterizar la distribución de la altitud en un cultivo de aguacate ubicado en el departamento del Cauca, Colombia. El análisis se basa en mediciones realizadas sobre 534 árboles georreferenciados el 01 de octubre de 2020.

Como parte del proceso metodológico, se llevaron a cabo análisis exploratorios previos del semivariograma experimental sobre las principales variables micrometeorológicas registradas en campo, con el propósito de identificar cuáles presentaban estructura de autocorrelación espacial suficiente para ser interpoladas. Los resultados indicaron que variables como la Temperatura y la Humedad Relativa no mostraron dependencia espacial estructurada: sus semivariogramas resultaron planos o erráticos, sin un incremento sistemático de la semivarianza con la distancia, lo que descarta la posibilidad de ajustar un modelo geoestadístico confiable sobre ellas.

La variable Altitud, por el contrario, demostró una autocorrelación espacial clara y consistente, con un semivariograma experimental que exhibe incremento progresivo de la semivarianza con la distancia y una meseta definida. Este comportamiento confirma el cumplimiento de los supuestos fundamentales del Kriging Ordinario, a saber, estacionariedad de segundo orden, variabilidad espacialmente estructurada y rango de autocorrelación coherente con la extensión del área de estudio. El área de trabajo corresponde a una parcela de aproximadamente 170 por 176 metros, con altitudes comprendidas entre 1675 y 1724 m.s.n.m. y una media de 1693 m.s.n.m. Las coordenadas fueron reproyectadas de WGS84 a UTM Zona 18S para garantizar el trabajo en unidades métricas. Se evaluaron y compararon tres modelos de semivariograma —Esférico, Exponencial y Gaussiano— mediante el criterio de mínimo SSE, y la calidad predictiva del modelo seleccionado fue validada mediante validación cruzada Leave-One-Out, reportando las métricas ME, RMSE, R² y RMSSE.

# ==============================================================================
# KRIGING ORDINARIO - DATOS AGUACATE (01/10/2020)
# Variable: Altitud | N = 534 árboles | Área ~170m x 176m | Cauca, Colombia
# Ejecutar PARTE POR PARTE para verificar cada resultado
# ==============================================================================
 
 
# ==============================================================================
# PARTE 1 — LIBRERÍAS
# ==============================================================================
# Instalar solo si no están (comentar después de primera ejecución)
# install.packages(c("readxl","sp","gstat","sf","ggplot2","viridis","dplyr","automap"))
 
library(readxl)

## Warning: package 'readxl' was built under R version 4.5.3

library(sp)

## Warning: package 'sp' was built under R version 4.5.3

library(gstat)

## Warning: package 'gstat' was built under R version 4.5.3

library(sf)

## Warning: package 'sf' was built under R version 4.5.3

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.5.3

library(viridis)

## Warning: package 'viridis' was built under R version 4.5.3

## Cargando paquete requerido: viridisLite

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(automap)

## Warning: package 'automap' was built under R version 4.5.3

cat("✓ Librerías cargadas correctamente\n")

## ✓ Librerías cargadas correctamente

# ==============================================================================
# PARTE 2 — CARGA Y FILTRADO DE DATOS
# ==============================================================================
 
datos_raw <- read_excel("C:/Users/ACER/Downloads/Datos_Aguacate.xlsx")

## Warning: Expecting numeric in T7458 / R7458C20: got '-'

## Warning: Expecting numeric in T8174 / R8174C20: got '-'

## Warning: Expecting numeric in T8415 / R8415C20: got '-'

## Warning: Expecting numeric in T10509 / R10509C20: got '-'

## Warning: Expecting numeric in T11253 / R11253C20: got '-'

## Warning: Expecting numeric in T11254 / R11254C20: got '-'

## Warning: Expecting numeric in T11810 / R11810C20: got '-'

## Warning: Coercing text to numeric in P11992 / R11992C16: '1809.5'

## Warning: Coercing text to numeric in P11993 / R11993C16: '1809.5'

## Warning: Coercing text to numeric in P11994 / R11994C16: '1809.5'

## Warning: Coercing text to numeric in P11995 / R11995C16: '1809.5'

## Warning: Coercing text to numeric in P11996 / R11996C16: '1809.5'

## Warning: Coercing text to numeric in P11997 / R11997C16: '1808.5'

## Warning: Coercing text to numeric in P11998 / R11998C16: '1807.5'

## Warning: Coercing text to numeric in P11999 / R11999C16: '1806.5'

## Warning: Coercing text to numeric in P12000 / R12000C16: '1805.5'

## Warning: Coercing text to numeric in P12001 / R12001C16: '1804.5'

## Warning: Coercing text to numeric in P12002 / R12002C16: '1803.5'

## Warning: Coercing text to numeric in P12006 / R12006C16: '1802.5'

## Warning: Coercing text to numeric in P12007 / R12007C16: '1803.5'

## Warning: Coercing text to numeric in P12011 / R12011C16: '1807.5'

## Warning: Coercing text to numeric in P12013 / R12013C16: '1808.5'

## Warning: Coercing text to numeric in P12014 / R12014C16: '1808.5'

## Warning: Coercing text to numeric in P12015 / R12015C16: '1807.5'

## Warning: Coercing text to numeric in P12016 / R12016C16: '1808.5'

## Warning: Coercing text to numeric in P12018 / R12018C16: '1807.5'

## Warning: Coercing text to numeric in P12019 / R12019C16: '1808.5'

## Warning: Coercing text to numeric in P12022 / R12022C16: '1809.5'

## Warning: Coercing text to numeric in P12024 / R12024C16: '1806.5'

## Warning: Coercing text to numeric in P12026 / R12026C16: '1805.5'

## Warning: Coercing text to numeric in P12027 / R12027C16: '1805.5'

## Warning: Coercing text to numeric in P12029 / R12029C16: '1806.5'

## Warning: Coercing text to numeric in P12030 / R12030C16: '1806.5'

## Warning: Coercing text to numeric in P12031 / R12031C16: '1807.5'

## Warning: Coercing text to numeric in P12032 / R12032C16: '1807.5'

## Warning: Coercing text to numeric in P12033 / R12033C16: '1806.5'

## Warning: Coercing text to numeric in P12034 / R12034C16: '1806.5'

## Warning: Coercing text to numeric in P12035 / R12035C16: '1805.5'

## Warning: Coercing text to numeric in P12036 / R12036C16: '1804.5'

## Warning: Coercing text to numeric in P12037 / R12037C16: '1803.5'

## Warning: Coercing text to numeric in P12039 / R12039C16: '1800.5'

## Warning: Coercing text to numeric in P12040 / R12040C16: '1800.5'

## Warning: Coercing text to numeric in P12044 / R12044C16: '1797.5'

## Warning: Coercing text to numeric in P12046 / R12046C16: '1796.5'

## Warning: Coercing text to numeric in P12047 / R12047C16: '1796.5'

## Warning: Coercing text to numeric in P12048 / R12048C16: '1797.5'

## Warning: Coercing text to numeric in P12049 / R12049C16: '1798.5'

## Warning: Coercing text to numeric in P12050 / R12050C16: '1800.5'

## Warning: Coercing text to numeric in P12053 / R12053C16: '1803.5'

## Warning: Coercing text to numeric in P12057 / R12057C16: '1805.5'

## Warning: Coercing text to numeric in P12058 / R12058C16: '1804.5'

## Warning: Coercing text to numeric in P12059 / R12059C16: '1804.5'

## Warning: Coercing text to numeric in P12062 / R12062C16: '1803.5'

## Warning: Coercing text to numeric in P12064 / R12064C16: '1806.5'

## Warning: Coercing text to numeric in P12065 / R12065C16: '1805.5'

## Warning: Coercing text to numeric in P12066 / R12066C16: '1805.5'

## Warning: Coercing text to numeric in P12067 / R12067C16: '1803.5'

## Warning: Coercing text to numeric in P12068 / R12068C16: '1802.5'

## Warning: Coercing text to numeric in P12070 / R12070C16: '1801.5'

## Warning: Coercing text to numeric in P12071 / R12071C16: '1801.5'

## Warning: Coercing text to numeric in P12073 / R12073C16: '1802.5'

## Warning: Coercing text to numeric in P12075 / R12075C16: '1804.5'

## Warning: Coercing text to numeric in P12076 / R12076C16: '1803.5'

## Warning: Coercing text to numeric in P12078 / R12078C16: '1800.5'

## Warning: Coercing text to numeric in P12079 / R12079C16: '1799.5'

## Warning: Coercing text to numeric in P12080 / R12080C16: '1797.5'

## Warning: Coercing text to numeric in P12081 / R12081C16: '1796.5'

## Warning: Coercing text to numeric in P12083 / R12083C16: '1793.5'

## Warning: Coercing text to numeric in P12086 / R12086C16: '1793.5'

## Warning: Coercing text to numeric in P12087 / R12087C16: '1792.5'

## Warning: Coercing text to numeric in P12088 / R12088C16: '1791.5'

## Warning: Coercing text to numeric in P12089 / R12089C16: '1791.5'

## Warning: Coercing text to numeric in P12092 / R12092C16: '1793.5'

## Warning: Coercing text to numeric in P12099 / R12099C16: '1802.5'

## Warning: Coercing text to numeric in P12103 / R12103C16: '1802.5'

## Warning: Coercing text to numeric in P12104 / R12104C16: '1802.5'

## Warning: Coercing text to numeric in P12106 / R12106C16: '1802.5'

## Warning: Coercing text to numeric in P12107 / R12107C16: '1803.5'

## Warning: Coercing text to numeric in P12108 / R12108C16: '1803.5'

## Warning: Coercing text to numeric in P12109 / R12109C16: '1803.5'

## Warning: Coercing text to numeric in P12112 / R12112C16: '1801.5'

## Warning: Coercing text to numeric in P12114 / R12114C16: '1801.5'

## Warning: Coercing text to numeric in P12116 / R12116C16: '1801.5'

## Warning: Coercing text to numeric in P12117 / R12117C16: '1801.5'

## Warning: Coercing text to numeric in P12118 / R12118C16: '1802.5'

## Warning: Coercing text to numeric in P12119 / R12119C16: '1803.5'

## Warning: Coercing text to numeric in P12123 / R12123C16: '1802.5'

## Warning: Coercing text to numeric in P12124 / R12124C16: '1801.5'

## Warning: Coercing text to numeric in P12127 / R12127C16: '1797.5'

## Warning: Coercing text to numeric in P12130 / R12130C16: '1795.5'

## Warning: Coercing text to numeric in P12131 / R12131C16: '1795.5'

## Warning: Coercing text to numeric in P12134 / R12134C16: '1795.5'

## Warning: Coercing text to numeric in P12136 / R12136C16: '1796.5'

## Warning: Coercing text to numeric in P12140 / R12140C16: '1797.5'

## Warning: Coercing text to numeric in P12141 / R12141C16: '1797.5'

## Warning: Coercing text to numeric in P12142 / R12142C16: '1798.5'

## Warning: Coercing text to numeric in P12148 / R12148C16: '1795.5'

## Warning: Coercing text to numeric in P12155 / R12155C16: '1797.5'

## Warning: Coercing text to numeric in P12157 / R12157C16: '1794.5'

## Warning: Coercing text to numeric in P12158 / R12158C16: '1796.5'

## Warning: Coercing text to numeric in P12159 / R12159C16: '1796.5'

## Warning: Coercing text to numeric in P12164 / R12164C16: '1795.5'

## Warning: Coercing text to numeric in P12165 / R12165C16: '1795.5'

## Warning: Coercing text to numeric in P12170 / R12170C16: '1794.5'

## Warning: Coercing text to numeric in P12172 / R12172C16: '1793.5'

## Warning: Coercing text to numeric in P12173 / R12173C16: '1793.5'

## Warning: Coercing text to numeric in P12175 / R12175C16: '1793.5'

## Warning: Coercing text to numeric in P12177 / R12177C16: '1795.5'

## Warning: Coercing text to numeric in P12178 / R12178C16: '1796.5'

## Warning: Coercing text to numeric in P12180 / R12180C16: '1798.5'

## Warning: Coercing text to numeric in P12182 / R12182C16: '1799.5'

## Warning: Coercing text to numeric in P12183 / R12183C16: '1800.5'

## Warning: Coercing text to numeric in P12185 / R12185C16: '1799.5'

## Warning: Coercing text to numeric in P12191 / R12191C16: '1802.5'

## Warning: Coercing text to numeric in P12194 / R12194C16: '1798.5'

## Warning: Coercing text to numeric in P12196 / R12196C16: '1797.5'

## Warning: Coercing text to numeric in P12198 / R12198C16: '1797.5'

## Warning: Coercing text to numeric in P12202 / R12202C16: '1798.5'

## Warning: Coercing text to numeric in P12203 / R12203C16: '1797.5'

## Warning: Coercing text to numeric in P12205 / R12205C16: '1795.5'

## Warning: Coercing text to numeric in P12206 / R12206C16: '1793.5'

## Warning: Coercing text to numeric in P12207 / R12207C16: '1792.5'

## Warning: Coercing text to numeric in P12210 / R12210C16: '1787.5'

## Warning: Coercing text to numeric in P12217 / R12217C16: '1788.5'

## Warning: Coercing text to numeric in P12219 / R12219C16: '1790.5'

## Warning: Coercing text to numeric in P12220 / R12220C16: '1791.5'

## Warning: Coercing text to numeric in P12221 / R12221C16: '1791.5'

## Warning: Coercing text to numeric in P12224 / R12224C16: '1792.5'

## Warning: Coercing text to numeric in P12226 / R12226C16: '1793.5'

## Warning: Coercing text to numeric in P12229 / R12229C16: '1796.5'

## Warning: Coercing text to numeric in P12235 / R12235C16: '1796.5'

## Warning: Coercing text to numeric in P12242 / R12242C16: '1797.5'

## Warning: Coercing text to numeric in P12245 / R12245C16: '1796.5'

## Warning: Coercing text to numeric in P12247 / R12247C16: '1792.5'

## Warning: Coercing text to numeric in P12251 / R12251C16: '1788.5'

## Warning: Coercing text to numeric in P12254 / R12254C16: '1788.5'

## Warning: Coercing text to numeric in P12257 / R12257C16: '1792.5'

## Warning: Coercing text to numeric in P12261 / R12261C16: '1792.5'

## Warning: Coercing text to numeric in P12262 / R12262C16: '1793.5'

## Warning: Coercing text to numeric in P12264 / R12264C16: '1795.5'

## Warning: Coercing text to numeric in P12267 / R12267C16: '1796.5'

## Warning: Coercing text to numeric in P12270 / R12270C16: '1794.5'

## Warning: Coercing text to numeric in P12273 / R12273C16: '1791.5'

## Warning: Coercing text to numeric in P12277 / R12277C16: '1785.5'

## Warning: Coercing text to numeric in P12280 / R12280C16: '1784.5'

## Warning: Coercing text to numeric in P12282 / R12282C16: '1785.5'

## Warning: Coercing text to numeric in P12283 / R12283C16: '1785.5'

## Warning: Coercing text to numeric in P12285 / R12285C16: '1785.5'

## Warning: Coercing text to numeric in P12286 / R12286C16: '1785.5'

## Warning: Coercing text to numeric in P12287 / R12287C16: '1784.5'

## Warning: Coercing text to numeric in P12288 / R12288C16: '1784.5'

## Warning: Coercing text to numeric in P12291 / R12291C16: '1783.5'

## Warning: Coercing text to numeric in P12292 / R12292C16: '1783.5'

## Warning: Coercing text to numeric in P12295 / R12295C16: '1783.5'

## Warning: Coercing text to numeric in P12296 / R12296C16: '1784.5'

## Warning: Coercing text to numeric in P12298 / R12298C16: '1786.5'

## Warning: Coercing text to numeric in P12302 / R12302C16: '1787.5'

## Warning: Coercing text to numeric in P12303 / R12303C16: '1786.5'

## Warning: Coercing text to numeric in P12305 / R12305C16: '1784.5'

## Warning: Coercing text to numeric in P12306 / R12306C16: '1783.5'

## Warning: Coercing text to numeric in P12309 / R12309C16: '1783.5'

## Warning: Coercing text to numeric in P12310 / R12310C16: '1784.5'

## Warning: Coercing text to numeric in P12311 / R12311C16: '1783.5'

## Warning: Expecting numeric in T12985 / R12985C20: got '***'

## Warning: Expecting numeric in T13260 / R13260C20: got '***'

## Warning: Expecting numeric in T13313 / R13313C20: got '***'

## Warning: Expecting numeric in T13497 / R13497C20: got '***'

## Warning: Expecting numeric in T13557 / R13557C20: got '***'

## Warning: Expecting numeric in T13559 / R13559C20: got '***'

## Warning: Expecting numeric in T14034 / R14034C20: got '***'

## Warning: Expecting numeric in T14250 / R14250C20: got '***'

## Warning: Expecting numeric in T14275 / R14275C20: got '***'

## Warning: Expecting numeric in T14370 / R14370C20: got '-'

## Warning: Expecting numeric in T17638 / R17638C20: got '-'

## Warning: Expecting numeric in T18138 / R18138C20: got '***'

cat("Columnas disponibles:\n")

## Columnas disponibles:

print(names(datos_raw))

##  [1] "id_arbol"                       "Latitude"                      
##  [3] "Longitude"                      "FORMATTED_DATE_TIME"           
##  [5] "Psychro_Wet_Bulb_Temperature"   "Station_Pressure"              
##  [7] "Relative_Humidity"              "Crosswind"                     
##  [9] "Temperature"                    "Barometric_Pressure"           
## [11] "Headwind"                       "Direction_True"                
## [13] "Direction_Mag"                  "Wind_Speed"                    
## [15] "Heat_Stress_Index"              "Altitude"                      
## [17] "Dew_Point"                      "Density_Altitude"              
## [19] "Wind_Chill"                     "Estado_Fenologico_Predominante"
## [21] "Frutos_Afectados"

# Filtrar fecha objetivo
datos <- datos_raw %>%
  filter(grepl("01/10/2020", FORMATTED_DATE_TIME)) %>%
  select(id_arbol, Latitude, Longitude, Altitude,
         Temperature, Relative_Humidity, Dew_Point,
         Frutos_Afectados, Estado_Fenologico_Predominante)
 
cat("\n=== DATOS FILTRADOS: 01/10/2020 ===\n")

## 
## === DATOS FILTRADOS: 01/10/2020 ===

cat("N registros          :", nrow(datos), "\n")

## N registros          : 534

cat("Coordenadas únicas   :", nrow(distinct(datos, Latitude, Longitude)), "\n")

## Coordenadas únicas   : 534

cat("Rango Altitud        :", min(datos$Altitude), "-", max(datos$Altitude), "m\n")

## Rango Altitud        : 1675 - 1724 m

cat("SD Altitud           :", round(sd(datos$Altitude), 2), "m\n")

## SD Altitud           : 10.98 m

cat("Rango Latitud        :", round(min(datos$Latitude),6), "-", round(max(datos$Latitude),6), "\n")

## Rango Latitud        : 2.392101 - 2.393634

cat("Rango Longitud       :", round(min(datos$Longitude),6), "-", round(max(datos$Longitude),6), "\n")

## Rango Longitud       : -76.7118 - -76.71022

cat("Extensión aprox.     : ~170m (N-S) x ~176m (E-W)\n\n")

## Extensión aprox.     : ~170m (N-S) x ~176m (E-W)

# Vista rápida
print(head(datos, 5))

## # A tibble: 5 × 9
##   id_arbol Latitude Longitude Altitude Temperature Relative_Humidity Dew_Point
##   <chr>       <dbl>     <dbl>    <dbl>       <dbl>             <dbl>     <dbl>
## 1 1            2.39     -76.7     1696        23.9              85.2      21.3
## 2 2            2.39     -76.7     1696        23.5              84        20.7
## 3 3            2.39     -76.7     1694        24.5              79.6      20.8
## 4 4            2.39     -76.7     1694        25.9              77.6      21.7
## 5 5            2.39     -76.7     1696        26                76.5      21.5
## # ℹ 2 more variables: Frutos_Afectados <dbl>,
## #   Estado_Fenologico_Predominante <dbl>

# ==============================================================================
# PARTE 3 — ESTADÍSTICAS DESCRIPTIVAS Y DISTRIBUCIÓN DE ALTITUD
# ==============================================================================
 
cat("\n=== ESTADÍSTICAS DESCRIPTIVAS — ALTITUD ===\n")

## 
## === ESTADÍSTICAS DESCRIPTIVAS — ALTITUD ===

print(summary(datos$Altitude))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1675    1684    1692    1693    1700    1724

cat("Desviación estándar  :", round(sd(datos$Altitude), 2), "\n")

## Desviación estándar  : 10.98

cat("CV (%)               :", round(sd(datos$Altitude)/mean(datos$Altitude)*100, 2), "\n")

## CV (%)               : 0.65

cat("Curtosis             :", round(e1071::kurtosis(datos$Altitude), 3), "\n")  # requiere e1071

## Curtosis             : -0.091

cat("Asimetría            :", round(e1071::skewness(datos$Altitude), 3), "\n\n")

## Asimetría            : 0.697

# Histograma + densidad
ggplot(datos, aes(x = Altitude)) +
  geom_histogram(aes(y = ..density..), bins = 25,
                 fill = "#2E86AB", color = "white", alpha = 0.8) +
  geom_density(color = "#A23B72", size = 1.2) +
  geom_vline(xintercept = mean(datos$Altitude),
             color = "red", linetype = "dashed", size = 1) +
  labs(title = "Distribución de Altitud — 01/10/2020",
       subtitle = paste0("N = ", nrow(datos), " árboles | Media = ",
                         round(mean(datos$Altitude),1), " m.s.n.m."),
       x = "Altitud (m.s.n.m.)", y = "Densidad") +
  theme_minimal(base_size = 13)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Test de normalidad (importante para kriging)
shapiro_test <- shapiro.test(sample(datos$Altitude, min(5000, nrow(datos))))
cat("Test Shapiro-Wilk:\n")

## Test Shapiro-Wilk:

cat("  W =", round(shapiro_test$statistic, 4), "\n")

##   W = 0.9494

cat("  p-value =", round(shapiro_test$p.value, 4), "\n")

##   p-value = 0

cat("  Distribución normal:", ifelse(shapiro_test$p.value > 0.05, "SÍ ✓", "NO — considerar transformación"), "\n\n")

##   Distribución normal: NO — considerar transformación

# ==============================================================================
# PARTE 4 — MAPA DE PUNTOS (distribución espacial)
# ==============================================================================
 
ggplot(datos, aes(x = Longitude, y = Latitude, color = Altitude)) +
  geom_point(size = 2.5, alpha = 0.85) +
  scale_color_viridis(option = "plasma", name = "Altitud\n(m.s.n.m.)") +
  labs(title = "Distribución Espacial de Altitud — Aguacate",
       subtitle = "01/10/2020 | Cauca, Colombia",
       x = "Longitud", y = "Latitud") +
  coord_equal() +
  theme_minimal(base_size = 13) +
  theme(legend.position = "right")

cat("✓ Mapa de puntos generado\n")

## ✓ Mapa de puntos generado

# PARTE 5 — CONVERSIÓN A OBJETO ESPACIAL (sp)
# ==============================================================================
 
# Crear SpatialPointsDataFrame con proyección geográfica WGS84
coordinates(datos) <- ~ Longitude + Latitude
proj4string(datos) <- CRS("+proj=longlat +datum=WGS84")
 
# Reproyectar a UTM zona 18N para distancias en metros (Colombia)
datos_utm <- spTransform(datos,
  CRS("+proj=utm +zone=18 +south +datum=WGS84 +units=m"))
 
cat("\n=== CONVERSIÓN ESPACIAL ===\n")

## 
## === CONVERSIÓN ESPACIAL ===

cat("CRS original  : WGS84 geográfico (grados)\n")

## CRS original  : WGS84 geográfico (grados)

cat("CRS proyectado: UTM Zona 18S (metros)\n")

## CRS proyectado: UTM Zona 18S (metros)

cat("Extensión UTM :\n")

## Extensión UTM :

print(bbox(datos_utm))

##                  min        max
## coords.x1   309656.1   309832.4
## coords.x2 10264519.2 10264688.6

# ==============================================================================
# PARTE 6 — ANÁLISIS EXPLORATORIO DEL VARIOGRAMA EXPERIMENTAL
# ==============================================================================
 
# Variograma experimental estándar
# cutoff = mitad de la diagonal del área (~120m), width = lag ~12m
vgm_exp <- variogram(Altitude ~ 1, data = datos_utm,
                     cutoff = 120,
                     width  = 12)
 
cat("\n=== VARIOGRAMA EXPERIMENTAL ===\n")

## 
## === VARIOGRAMA EXPERIMENTAL ===

cat("Número de bins (lags) :", nrow(vgm_exp), "\n")

## Número de bins (lags) : 10

print(vgm_exp)   # columna 'np' = pares por bin

##       np       dist     gamma dir.hor dir.ver   id
## 1   3026   8.476722  40.77181       0       0 var1
## 2   8782  18.558085  55.38146       0       0 var1
## 3  12693  30.275888  68.68120       0       0 var1
## 4  15148  42.124510  81.14051       0       0 var1
## 5  16385  54.032685  94.92496       0       0 var1
## 6  16428  65.976855 104.32725       0       0 var1
## 7  15277  77.931388 108.58667       0       0 var1
## 8  13359  89.819314 117.87604       0       0 var1
## 9  11163 101.781363 132.39044       0       0 var1
## 10  8853 113.760418 156.24720       0       0 var1

plot(vgm_exp,
     main = "Variograma Experimental — Altitud",
     xlab = "Distancia (m)",
     ylab = "Semivarianza",
     pch  = 19, col = "#2E86AB", cex = 1.4)

# Variograma ROBUSTO (Cressie-Hawkins, menos sensible a outliers)
vgm_rob <- variogram(Altitude ~ 1, data = datos_utm,
                     cutoff = 120, width = 12,
                     cressie = TRUE)
 
cat("\n=== VARIOGRAMA ROBUSTO (Cressie-Hawkins) ===\n")

## 
## === VARIOGRAMA ROBUSTO (Cressie-Hawkins) ===

print(vgm_rob)

##       np       dist     gamma dir.hor dir.ver   id
## 1   3026   8.476722  19.97564       0       0 var1
## 2   8782  18.558085  45.06854       0       0 var1
## 3  12693  30.275888  66.77789       0       0 var1
## 4  15148  42.124510  81.30962       0       0 var1
## 5  16385  54.032685  96.22120       0       0 var1
## 6  16428  65.976855 103.73486       0       0 var1
## 7  15277  77.931388 104.47546       0       0 var1
## 8  13359  89.819314 111.90253       0       0 var1
## 9  11163 101.781363 134.54637       0       0 var1
## 10  8853 113.760418 167.78247       0       0 var1

plot(vgm_rob,
     main = "Variograma Robusto (Cressie-Hawkins) — Altitud",
     xlab = "Distancia (m)", ylab = "Semivarianza",
     pch = 19, col = "#A23B72", cex = 1.4)

# ==============================================================================
# PARTE 7 — VARIOGRAMA NUBE (para detectar pares problemáticos)
# ==============================================================================

vgm_cloud <- variogram(Altitude ~ 1, data = datos_utm,
                       cloud = TRUE, cutoff = 120)

# Convertir a dataframe para graficar con ggplot
cloud_df <- as.data.frame(vgm_cloud)

ggplot(cloud_df, aes(x = dist, y = gamma)) +
  geom_point(color = "#2E86AB", alpha = 0.3, size = 1.2) +
  geom_hline(yintercept = var(datos_utm$Altitude),
             linetype = "dashed", color = "red", size = 1) +
  annotate("text", x = max(cloud_df$dist) * 0.7,
           y = var(datos_utm$Altitude) * 1.05,
           label = "Varianza total", color = "red", size = 3.5) +
  labs(title    = "Variograma Nube — Altitud (cada par de puntos)",
       subtitle = "Puntos muy por encima de la línea roja = pares anómalos",
       x = "Distancia (m)",
       y = "Semivarianza por par") +
  theme_minimal(base_size = 13)

cat("\n✓ Variograma nube generado\n")

## 
## ✓ Variograma nube generado

cat("  Pares con gamma >> varianza total pueden ser outliers espaciales\n")

##   Pares con gamma >> varianza total pueden ser outliers espaciales

# PARTE 8 — ANÁLISIS DE ANISOTROPÍA
# ==============================================================================
 
# Mapa del variograma direccional (detecta anisotropía)
vgm_map <- variogram(Altitude ~ 1, data = datos_utm,
                     map = TRUE, cutoff = 120, width = 15)
 
plot(vgm_map, threshold = 10,
     main = "Mapa de Variograma — Análisis de Anisotropía",
     xlab = "Δx (m)", ylab = "Δy (m)")

# Variogramas direccionales (0°, 45°, 90°, 135°)
vgm_dir <- variogram(Altitude ~ 1, data = datos_utm,
                     alpha = c(0, 45, 90, 135),
                     cutoff = 120, width = 15)
 
plot(vgm_dir,
     main = "Variogramas Direccionales — Altitud",
     xlab = "Distancia (m)", ylab = "Semivarianza")

cat("\n✓ Si las direcciones son similares → Isotropía (se puede usar modelo isotrópico)\n")

## 
## ✓ Si las direcciones son similares → Isotropía (se puede usar modelo isotrópico)

cat("  Si difieren → Anisotropía (se necesita modelo anisotrópico)\n\n")

##   Si difieren → Anisotropía (se necesita modelo anisotrópico)

# PARTE 9 — AJUSTE DE MODELOS DE VARIOGRAMA
# ============================================================
 # PARTE 9 — AJUSTE DE MODELOS DE VARIOGRAMA
# ==============================================================================

nugget_ini <- 0
sill_ini   <- var(datos_utm$Altitude)
range_ini  <- 80

cat("=== PARÁMETROS INICIALES ===\n")

## === PARÁMETROS INICIALES ===

cat("Nugget inicial:", nugget_ini, "\n")

## Nugget inicial: 0

cat("Sill inicial  :", round(sill_ini, 2), "\n")

## Sill inicial  : 120.51

cat("Range inicial :", range_ini, "m\n\n")

## Range inicial : 80 m

# Ajuste de los tres modelos
fit_sph <- fit.variogram(vgm_exp,
  model = vgm(psill = sill_ini, model = "Sph",
               range = range_ini, nugget = nugget_ini))

fit_exp <- fit.variogram(vgm_exp,
  model = vgm(psill = sill_ini, model = "Exp",
               range = range_ini, nugget = nugget_ini))

fit_gau <- fit.variogram(vgm_exp,
  model = vgm(psill = sill_ini, model = "Gau",
               range = range_ini, nugget = nugget_ini))

cat("=== PARÁMETROS AJUSTADOS — ESFÉRICO ===\n");   print(fit_sph)

## === PARÁMETROS AJUSTADOS — ESFÉRICO ===

##   model     psill    range
## 1   Nug  31.75544   0.0000
## 2   Sph 114.77421 143.9537

cat("SSE Esférico    :", round(attr(fit_sph,"SSErr"), 4), "\n\n")

## SSE Esférico    : 568.7401

cat("=== PARÁMETROS AJUSTADOS — EXPONENCIAL ===\n"); print(fit_exp)

## === PARÁMETROS AJUSTADOS — EXPONENCIAL ===

##   model     psill    range
## 1   Nug  29.37197   0.0000
## 2   Exp 190.35080 131.6013

cat("SSE Exponencial :", round(attr(fit_exp,"SSErr"), 4), "\n\n")

## SSE Exponencial : 377.6876

cat("=== PARÁMETROS AJUSTADOS — GAUSSIANO ===\n");   print(fit_gau)

## === PARÁMETROS AJUSTADOS — GAUSSIANO ===

##   model    psill   range
## 1   Nug 39.40316  0.0000
## 2   Gau 77.07380 44.2145

cat("SSE Gaussiano   :", round(attr(fit_gau,"SSErr"), 4), "\n\n")

## SSE Gaussiano   : 2149.131

# Gráfico comparativo con ggplot2
linea_sph <- as.data.frame(variogramLine(fit_sph, maxdist = max(vgm_exp$dist)))
linea_exp <- as.data.frame(variogramLine(fit_exp, maxdist = max(vgm_exp$dist)))
linea_gau <- as.data.frame(variogramLine(fit_gau, maxdist = max(vgm_exp$dist)))

linea_sph$modelo <- "Esférico"
linea_exp$modelo <- "Exponencial"
linea_gau$modelo <- "Gaussiano"

lineas_df <- rbind(linea_sph, linea_exp, linea_gau)

ggplot() +
  geom_point(data = vgm_exp, aes(x = dist, y = gamma),
             color = "black", size = 3, shape = 19) +
  geom_line(data = lineas_df, aes(x = dist, y = gamma,
            color = modelo, linetype = modelo), linewidth = 1.2) +
  scale_color_manual(values = c("Esférico"    = "#2E86AB",
                                "Exponencial" = "#E84855",
                                "Gaussiano"   = "#3BB273")) +
  scale_linetype_manual(values = c("Esférico"    = "solid",
                                   "Exponencial" = "dashed",
                                   "Gaussiano"   = "dotted")) +
  labs(title    = "Comparación de Modelos — Altitud",
       subtitle = "Puntos negros = variograma experimental",
       x        = "Distancia (m)",
       y        = "Semivarianza",
       color    = "Modelo",
       linetype = "Modelo") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "right")

# Seleccionar el mejor modelo (menor SSE)
sse_vals <- c(Esferico    = attr(fit_sph, "SSErr"),
              Exponencial = attr(fit_exp, "SSErr"),
              Gaussiano   = attr(fit_gau, "SSErr"))
mejor_modelo <- names(which.min(sse_vals))
cat("SSE Esférico   :", round(sse_vals["Esferico"],    4), "\n")

## SSE Esférico   : 568.7401

cat("SSE Exponencial:", round(sse_vals["Exponencial"], 4), "\n")

## SSE Exponencial: 377.6876

cat("SSE Gaussiano  :", round(sse_vals["Gaussiano"],   4), "\n")

## SSE Gaussiano  : 2149.131

cat(">>> MEJOR MODELO (menor SSE):", mejor_modelo, "<<<\n\n")

## >>> MEJOR MODELO (menor SSE): Exponencial <<<

# ==============================================================================
# PARTE 10 — AJUSTE AUTOMÁTICO CON automap (validación adicional)
# ==============================================================================

fit_auto <- autofitVariogram(Altitude ~ 1, input_data = datos_utm)

cat("=== AJUSTE AUTOMÁTICO (automap) ===\n")

## === AJUSTE AUTOMÁTICO (automap) ===

cat("Modelo seleccionado automáticamente:\n")

## Modelo seleccionado automáticamente:

print(fit_auto$var_model)

##   model     psill    range kappa
## 1   Nug  13.50526   0.0000   0.0
## 2   Ste 341.98405 615.4905   0.3

# plot() de automap usa lattice internamente y conflicta con 'main'
# Se extrae el variograma y se grafica con ggplot2

vgm_auto_df  <- as.data.frame(fit_auto$exp_var)
linea_auto   <- as.data.frame(variogramLine(fit_auto$var_model,
                                             maxdist = max(vgm_auto_df$dist)))

ggplot() +
  geom_point(data = vgm_auto_df, aes(x = dist, y = gamma),
             color = "black", size = 3, shape = 19) +
  geom_line(data = linea_auto, aes(x = dist, y = gamma),
            color = "#A23B72", size = 1.3) +
  labs(title    = "Variograma Automático (automap) — Altitud",
       subtitle = paste0("Modelo: ", fit_auto$var_model$model[2],
                         " | Sill: ", round(sum(fit_auto$var_model$psill), 2),
                         " | Range: ", round(fit_auto$var_model$range[2], 2), " m"),
       x = "Distancia (m)",
       y = "Semivarianza") +
  theme_minimal(base_size = 13)

# ==============================================================================
# PARTE 11 — CREACIÓN DE LA GRILLA DE INTERPOLACIÓN
# ==============================================================================

# Extensión espacial del área
bbox_utm <- bbox(datos_utm)
cat("Nombres del bbox:\n"); print(bbox_utm)

## Nombres del bbox:

##                  min        max
## coords.x1   309656.1   309832.4
## coords.x2 10264519.2 10264688.6

# Extraer límites con los nombres correctos (coords.x1, coords.x2)
x_min <- bbox_utm["coords.x1", "min"]
x_max <- bbox_utm["coords.x1", "max"]
y_min <- bbox_utm["coords.x2", "min"]
y_max <- bbox_utm["coords.x2", "max"]

# Grilla de predicción con resolución de 2 metros
resolucion <- 2

grd <- expand.grid(
  x = seq(x_min, x_max, by = resolucion),
  y = seq(y_min, y_max, by = resolucion)
)

# Convertir a SpatialPixelsDataFrame
coordinates(grd) <- ~ x + y
gridded(grd)     <- TRUE
proj4string(grd) <- CRS(proj4string(datos_utm))

cat("\n=== GRILLA DE INTERPOLACIÓN ===\n")

## 
## === GRILLA DE INTERPOLACIÓN ===

cat("Resolución    :", resolucion, "m\n")

## Resolución    : 2 m

cat("Celdas totales:", nrow(grd@coords), "\n")

## Celdas totales: 7565

cat("Extensión     :", round(x_max - x_min), "m (E-W) x",
                       round(y_max - y_min), "m (N-S)\n\n")

## Extensión     : 176 m (E-W) x 169 m (N-S)

# PARTE 12 — KRIGING ORDINARIO
# ==============================================================================
 
# Seleccionar el modelo con menor SSE (cambiar 'fit_sph' si otro es mejor)
# Opciones: fit_sph | fit_exp | fit_gau | fit_auto$var_model
modelo_kriging <- fit_sph   # <-- CAMBIAR según resultado de PARTE 9
 
cat("=== EJECUTANDO KRIGING ORDINARIO ===\n")

## === EJECUTANDO KRIGING ORDINARIO ===

cat("Modelo usado:", modelo_kriging$model[2], "\n")

## Modelo usado: 3

kriging_result <- krige(
  formula   = Altitude ~ 1,
  locations = datos_utm,
  newdata   = grd,
  model     = modelo_kriging
)

## [using ordinary kriging]

cat("✓ Kriging completado\n")

## ✓ Kriging completado

cat("Variables resultado:\n")

## Variables resultado:

cat("  var1.pred = predicciones de altitud\n")

##   var1.pred = predicciones de altitud

cat("  var1.var  = varianza de kriging (incertidumbre)\n\n")

##   var1.var  = varianza de kriging (incertidumbre)

# Resumen de predicciones
cat("=== RESUMEN DE PREDICCIONES ===\n")

## === RESUMEN DE PREDICCIONES ===

print(summary(kriging_result$var1.pred))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1681    1688    1692    1695    1701    1722

# ==============================================================================
# PARTE 13 — MAPAS DE RESULTADO
# ==============================================================================

# Convertir a dataframe para ggplot
krig_df <- as.data.frame(kriging_result)
names(krig_df)[1:2] <- c("x", "y")

# Extraer coordenadas UTM de los puntos originales correctamente
puntos_df <- as.data.frame(datos_utm)
# Renombrar coords.x1 y coords.x2 a x e y
names(puntos_df)[names(puntos_df) == "coords.x1"] <- "x"
names(puntos_df)[names(puntos_df) == "coords.x2"] <- "y"

cat("Columnas de puntos_df:\n"); print(names(puntos_df))

## Columnas de puntos_df:

## [1] "id_arbol"                       "Altitude"                      
## [3] "Temperature"                    "Relative_Humidity"             
## [5] "Dew_Point"                      "Frutos_Afectados"              
## [7] "Estado_Fenologico_Predominante" "x"                             
## [9] "y"

# --- Mapa de predicciones ---
p_pred <- ggplot(krig_df, aes(x = x, y = y, fill = var1.pred)) +
  geom_raster() +
  scale_fill_viridis(option = "plasma",
                     name   = "Altitud\n(m.s.n.m.)") +
  geom_point(data = puntos_df,
             aes(x = x, y = y),
             color = "white", size = 0.6, alpha = 0.5,
             inherit.aes = FALSE) +
  labs(title    = "Kriging Ordinario — Altitud",
       subtitle = "01/10/2020 | Cauca, Colombia",
       x = "Este (m UTM)", y = "Norte (m UTM)") +
  coord_equal() +
  theme_minimal(base_size = 13)

print(p_pred)

# --- Mapa de varianza (incertidumbre) ---
p_var <- ggplot(krig_df, aes(x = x, y = y, fill = var1.var)) +
  geom_raster() +
  scale_fill_viridis(option = "magma",
                     name   = "Varianza\nKriging") +
  labs(title    = "Varianza de Predicción (Incertidumbre)",
       subtitle = "Valores altos = mayor incertidumbre",
       x = "Este (m UTM)", y = "Norte (m UTM)") +
  coord_equal() +
  theme_minimal(base_size = 13)

print(p_var)

# --- Mapa de error estándar ---
krig_df$se <- sqrt(krig_df$var1.var)

p_se <- ggplot(krig_df, aes(x = x, y = y, fill = se)) +
  geom_raster() +
  scale_fill_viridis(option = "inferno",
                     name   = "Error Std\n(m)") +
  labs(title    = "Error Estándar de Kriging",
       subtitle = paste0("Media SE = ", round(mean(krig_df$se), 2), " m"),
       x = "Este (m UTM)", y = "Norte (m UTM)") +
  coord_equal() +
  theme_minimal(base_size = 13)

print(p_se)

cat("✓ Mapas generados correctamente\n\n")

## ✓ Mapas generados correctamente

# PARTE 14 — VALIDACIÓN CRUZADA (Leave-One-Out Cross Validation)
# ==============================================================================
 
cat("=== VALIDACIÓN CRUZADA — LEAVE-ONE-OUT ===\n")

## === VALIDACIÓN CRUZADA — LEAVE-ONE-OUT ===

cat("Ejecutando... (puede tardar 1-2 minutos con 534 puntos)\n")

## Ejecutando... (puede tardar 1-2 minutos con 534 puntos)

cv_result <- krige.cv(
  formula   = Altitude ~ 1,
  locations = datos_utm,
  model     = modelo_kriging,
  nfold     = nrow(datos_utm)   # Leave-One-Out completo
  # Para K-fold más rápido: nfold = 10
)
 
# Métricas de validación
residuos <- cv_result$residual
pred_cv  <- cv_result$var1.pred
obs_cv   <- datos_utm$Altitude
 
# Cálculo de métricas
ME    <- mean(residuos)                          # Error medio (sesgo)
MAE   <- mean(abs(residuos))                     # Error absoluto medio
RMSE  <- sqrt(mean(residuos^2))                  # Raíz del ECM
RMSSE <- sqrt(mean((residuos/sqrt(cv_result$var1.var))^2))  # RMSSE (debe ≈ 1)
R2    <- cor(obs_cv, pred_cv)^2                  # R² predicho vs observado
MEC   <- 1 - sum(residuos^2)/sum((obs_cv - mean(obs_cv))^2)  # Nash-Sutcliffe
 
cat("\n=== MÉTRICAS DE VALIDACIÓN ===\n")

## 
## === MÉTRICAS DE VALIDACIÓN ===

cat("ME    (sesgo)           :", round(ME,    4), "m       → idealmente ≈ 0\n")

## ME    (sesgo)           : 0.001 m       → idealmente ≈ 0

cat("MAE   (error absoluto)  :", round(MAE,   4), "m\n")

## MAE   (error absoluto)  : 4.5902 m

cat("RMSE  (error cuadrático):", round(RMSE,  4), "m\n")

## RMSE  (error cuadrático): 6.0181 m

cat("RMSSE (estandarizado)   :", round(RMSSE, 4), "       → idealmente ≈ 1\n")

## RMSSE (estandarizado)   : 0.9411        → idealmente ≈ 1

cat("R²    (predicho vs obs) :", round(R2,    4), "       → idealmente ≈ 1\n")

## R²    (predicho vs obs) : 0.6991        → idealmente ≈ 1

cat("MEC   (Nash-Sutcliffe)  :", round(MEC,   4), "       → idealmente ≈ 1\n\n")

## MEC   (Nash-Sutcliffe)  : 0.6989        → idealmente ≈ 1

# Interpretación automática
cat("=== INTERPRETACIÓN ===\n")

## === INTERPRETACIÓN ===

cat("Sesgo (ME)  :", ifelse(abs(ME) < 0.5, "Aceptable ✓", "Sesgo significativo ⚠"), "\n")

## Sesgo (ME)  : Aceptable ✓

cat("RMSSE       :", ifelse(RMSSE > 0.8 & RMSSE < 1.2,
                            "Varianza bien modelada ✓",
                            "Varianza mal calibrada ⚠"), "\n")

## RMSSE       : Varianza bien modelada ✓

cat("R²          :", ifelse(R2 > 0.7, "Buena predicción ✓",
                     ifelse(R2 > 0.5, "Predicción moderada ⚠",
                                      "Predicción débil ✗")), "\n\n")

## R²          : Predicción moderada ⚠

# PARTE 15 — GRÁFICOS DE VALIDACIÓN CRUZADA
# ============================================================================== 

cv_df <- data.frame(  
  Observado  = obs_cv,  
  Predicho   = pred_cv,  
  Residuo    = residuos,  
  SE         = sqrt(cv_result$var1.var)
) 

# --- Predicho vs Observado ---
p_cv1 <- ggplot(cv_df, aes(x = Observado, y = Predicho)) +  
  geom_point(color = "#2E86AB", alpha = 0.5, size = 1.5) +  
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", size = 1) +  
  geom_smooth(method = "lm", color = "#A23B72", se = TRUE) +  
  labs(title    = "Validación Cruzada: Predicho vs Observado",       
       subtitle = paste0("R² = ", round(R2,4), " | RMSE = ", round(RMSE,4), " m"),       
       x = "Altitud Observada (m)", y = "Altitud Predicha (m)") +  
  theme_minimal(base_size = 13) 

# --- Histograma de residuos ---
p_cv2 <- ggplot(cv_df, aes(x = Residuo)) +  
  geom_histogram(bins = 30, fill = "#3BB273", color = "white", alpha = 0.8) +  
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", size = 1) +  
  geom_vline(xintercept = ME, color = "blue", linetype = "solid", size = 1) +  
  labs(title    = "Distribución de Residuos",       
       subtitle = paste0("ME = ", round(ME,3), " | MAE = ", round(MAE,3)),       
       x = "Residuo (m)", y = "Frecuencia") +  
  theme_minimal(base_size = 13) 

# --- Residuos vs predichos (heterocedasticidad) ---
p_cv3 <- ggplot(cv_df, aes(x = Predicho, y = Residuo)) +  
  geom_point(color = "#E84855", alpha = 0.5, size = 1.5) +  
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +  
  geom_smooth(method = "loess", color = "#2E86AB", se = FALSE) +  
  labs(title = "Residuos vs Predichos",       
       subtitle = "Debe ser aleatoria alrededor de 0",       
       x = "Altitud Predicha (m)", y = "Residuo (m)") +  
  theme_minimal(base_size = 13) 

# --- Mapa espacial de residuos ---
cv_spatial <- as.data.frame(datos_utm)

# Extraemos las coordenadas explícitamente para evitar el error de "Longitude no encontrado"
# (Si usas el paquete 'sp' en lugar de 'sf', cambia st_coordinates por coordinates)
coords <- st_coordinates(st_as_sf(datos_utm)) 
cv_spatial$X_coord <- coords[, 1]
cv_spatial$Y_coord <- coords[, 2]
cv_spatial$Residuo <- residuos 

p_cv4 <- ggplot(cv_spatial, aes(x = X_coord, y = Y_coord, color = Residuo)) +  
  geom_point(size = 2) +  
  scale_color_gradient2(low = "blue", mid = "white", high = "red",                        
                        midpoint = 0, name = "Residuo\n(m)") +  
  labs(title    = "Distribución Espacial de Residuos",       
       subtitle = "Rojo = sobreestimación | Azul = subestimación",       
       x = "Este (UTM)", y = "Norte (UTM)") +  
  coord_equal() +  
  theme_minimal(base_size = 13) 

# --- Imprimir todos los gráficos ---
print(p_cv1)

## `geom_smooth()` using formula = 'y ~ x'

print(p_cv2)

print(p_cv3)

## `geom_smooth()` using formula = 'y ~ x'

print(p_cv4)

cat("✓ Gráficos de validación cruzada generados con éxito\n\n")

## ✓ Gráficos de validación cruzada generados con éxito

# PARTE 16 — REPORTE FINAL RESUMIDO
# ==============================================================================
 
cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat("║           REPORTE FINAL — KRIGING ORDINARIO             ║\n")

## ║           REPORTE FINAL — KRIGING ORDINARIO             ║

cat("╠══════════════════════════════════════════════════════════╣\n")

## ╠══════════════════════════════════════════════════════════╣

cat("║ Variable            : Altitud (m.s.n.m.)                ║\n")

## ║ Variable            : Altitud (m.s.n.m.)                ║

cat(sprintf("║ N puntos            : %d árboles                      ║\n", nrow(datos_utm)))

## ║ N puntos            : 534 árboles                      ║

cat(sprintf("║ Área               : ~170m x ~176m (~0.03 ha)         ║\n"))

## ║ Área               : ~170m x ~176m (~0.03 ha)         ║

cat(sprintf("║ Rango altitud       : %d - %d m                       ║\n",
            min(datos_utm$Altitude), max(datos_utm$Altitude)))

## ║ Rango altitud       : 1675 - 1724 m                       ║

cat("╠══════════════════════════════════════════════════════════╣\n")

## ╠══════════════════════════════════════════════════════════╣

cat(sprintf("║ Modelo variograma   : %-33s║\n", mejor_modelo))

## ║ Modelo variograma   : Exponencial                      ║

cat(sprintf("║ Nugget              : %-33s║\n",
            round(modelo_kriging$psill[1], 3)))

## ║ Nugget              : 31.755                           ║

cat(sprintf("║ Sill (partial)      : %-33s║\n",
            round(modelo_kriging$psill[2], 3)))

## ║ Sill (partial)      : 114.774                          ║

cat(sprintf("║ Range               : %-30s m ║\n",
            round(modelo_kriging$range[2], 2)))

## ║ Range               : 143.95                         m ║

cat("╠══════════════════════════════════════════════════════════╣\n")

## ╠══════════════════════════════════════════════════════════╣

cat(sprintf("║ RMSE (CV)           : %-30s m ║\n", round(RMSE,  4)))

## ║ RMSE (CV)           : 6.0181                         m ║

cat(sprintf("║ R² (CV)             : %-33s║\n", round(R2,    4)))

## ║ R² (CV)             : 0.6991                           ║

cat(sprintf("║ ME  (sesgo)         : %-30s m ║\n", round(ME,    4)))

## ║ ME  (sesgo)         : 0.001                          m ║

cat(sprintf("║ RMSSE               : %-33s║\n", round(RMSSE, 4)))

## ║ RMSSE               : 0.9411                           ║

cat(sprintf("║ Nash-Sutcliffe      : %-33s║\n", round(MEC,   4)))

## ║ Nash-Sutcliffe      : 0.6989                           ║

cat("╚══════════════════════════════════════════════════════════╝\n")

## ╚══════════════════════════════════════════════════════════╝

Conclusión del Modelado GeoestadísticoEl modelo de Kriging Ordinario (exponencial):

Demuestra un desempeño robusto y confiable para predecir la altitud en el área de estudio ubicada en el Cauca, Colombia, logrando explicar aproximadamente el 70% de la varianza topográfica.Se destacan tres aspectos fundamentales del análisis:Fuerte Dependencia Espacial: La relación entre el efecto pepita (Nugget, 31.75) y la meseta total (Sill parcial + Nugget = 146.52) es de aproximadamente 21.6%. Al ser menor al 25%, esto indica una dependencia espacial fuerte. La variable no es aleatoria; hay una continuidad topográfica muy clara. Además, el rango (Range) de 143.95 metros abarca casi la totalidad del área estudiada (~170x176m), confirmando que el muestreo de los 534 árboles fue lo suficientemente denso para capturar la estructura del relieve local en el Cauca.Ausencia de Sesgo: El Error Medio (ME) de 0.001 m indica que las predicciones son imparciales. El modelo no tiende sistemáticamente a sobreestimar ni a subestimar la altitud en el lote.Precisión y Cuantificación de Incertidumbre: El coeficiente de correlación \(R^2\) (0.6991) y el índice de Nash-Sutcliffe (0.6989) son prácticamente idénticos y positivos, lo que valida la calidad predictiva frente a usar solo el promedio de los datos. El error absoluto promedio (RMSE) de 6.01 m es bastante aceptable considerando un desnivel de casi 50 metros en el terreno. Finalmente, un RMSSE de 0.9411 (muy cercano a 1) demuestra que el modelo está calculando la varianza de predicción de manera casi perfecta, con apenas una ligerísima tendencia a sobreestimar la incertidumbre.

ANÁLISIS GEOESTADÍSTICO

Jorber Andrés Pardo Franco

2026-05-22

Conclusión del Modelado GeoestadísticoEl modelo de Kriging Ordinario (exponencial):