Universidad Pontificia Javeriana Cali
Informe 3 RMarkdown
Análisis Espacial del Cultivo de Aguacate
Presentado por:
Dainer Orlando Cruz Chate
Unidad 1
Análisis de Información Geográfica y Espacial
Profesor:
David Arango Londoño
Noviembre de 2026
1 Introducción
El análisis espacial permite evaluar si observaciones próximas presentan mayor semejanza que observaciones alejadas. La autocorrelación espacial global puede cuantificarse con el estadístico de Moran [1], mientras que el semivariograma representa la dependencia de segundo orden en función de la separación espacial [2]. Cuando esta dependencia está sustentada, el kriging produce predicciones lineales ponderadas y una medida de incertidumbre basada en un modelo de covarianza [2], [3].
En el cultivo de aguacate, las condiciones microclimáticas pueden variar dentro y entre lotes; por ello, la interpolación debe apoyarse en diagnóstico de calidad, evaluación de tendencia, prueba de autocorrelación y validación cruzada, no únicamente en una superficie visual. Este informe aplica esa secuencia al primer periodo cronológico de la base suministrada. Se prioriza temperatura y, solo si su autocorrelación no es significativa, se evalúan humedad relativa, velocidad del viento y altura, en ese orden.
2 Objetivos
2.1 Objetivo general
Analizar geoestadísticamente la información del primer periodo disponible y generar una predicción espacial mediante kriging para la primera variable candidata con autocorrelación espacial significativa.
2.2 Objetivos específicos
- Explorar variables numéricas y no numéricas, documentar su estructura y calidad.
- Extraer y depurar el primer periodo cronológico.
- Evaluar tendencia y autocorrelación espacial mediante Moran, correlograma y semivariograma.
- Comparar modelos exponencial, gaussiano y esférico mediante validación cruzada.
- Ajustar el modelo seleccionado, interpolar y representar predicción e incertidumbre.
3 Datos y metodología
datos_raw <- readxl::read_excel(archivo, sheet=1)
# La cadena original contiene variaciones de a. m./p. m.; fecha y hora se extraen por separado.
fecha_chr <- stringr::str_extract(as.character(datos_raw$FORMATTED_DATE_TIME), "\\d{2}/\\d{2}/\\d{4}")
hora_chr <- stringr::str_extract(as.character(datos_raw$FORMATTED_DATE_TIME), "\\d{1,2}:\\d{2}:\\d{2}")
datos <- datos_raw %>%
mutate(Fecha = lubridate::dmy(fecha_chr), Hora = hora_chr,
Estado_Fenologico_Predominante = as.character(Estado_Fenologico_Predominante))
primer_periodo <- min(datos$Fecha, na.rm=TRUE)
datos_p1 <- datos %>% filter(Fecha == primer_periodo)Se define periodo como cada fecha de medición
identificable en FORMATTED_DATE_TIME; por tanto, el primer
periodo es 30/07/2019 y contiene 427 registros. Esta decisión es
reproducible y evita imponer una periodicidad no documentada.
3.1 Exploración de todas las variables
resumen_variable <- function(x, nombre) {
tibble(Variable=nombre, Clase=paste(class(x), collapse="/"), N=length(x),
Faltantes=sum(is.na(x)), Unicos=dplyr::n_distinct(x, na.rm=TRUE),
Min=if(is.numeric(x)) signif(min(x,na.rm=TRUE),5) else NA,
Media=if(is.numeric(x)) signif(mean(x,na.rm=TRUE),5) else NA,
Mediana=if(is.numeric(x)) signif(median(x,na.rm=TRUE),5) else NA,
Max=if(is.numeric(x)) signif(max(x,na.rm=TRUE),5) else NA,
Categoria_mas_frecuente=if(!is.numeric(x)) names(sort(table(x),decreasing=TRUE))[1] else NA_character_)
}
resumen_total <- purrr::imap_dfr(datos, ~resumen_variable(.x,.y))
resumen_total %>% knitr::kable(caption="Tabla 1. Resumen exploratorio de todas las variables", digits=3) %>%
kableExtra::kable_styling(full_width=FALSE, position="center")| Variable | Clase | N | Faltantes | Unicos | Min | Media | Mediana | Max | Categoria_mas_frecuente |
|---|---|---|---|---|---|---|---|---|---|
| id_arbol | character | 20271 | 0 | 535 | NA | NA | NA | NA | 1 |
| Latitude | numeric | 20271 | 0 | 1431 | 2.316 | 2.367 | 2.373 | 2.394 | NA |
| Longitude | numeric | 20271 | 0 | 1431 | -76.712 | -76.657 | -76.614 | -76.607 | NA |
| FORMATTED_DATE_TIME | character | 20271 | 0 | 18626 | NA | NA | NA | NA | 01/10/2020 10:11:12 a, m, |
| Psychro_Wet_Bulb_Temperature | numeric | 20271 | 0 | 212 | 11.100 | 18.094 | 18.100 | 27.900 | NA |
| Station_Pressure | numeric | 20271 | 0 | 342 | 798.700 | 811.860 | 807.400 | 830.100 | NA |
| Relative_Humidity | numeric | 20271 | 0 | 864 | 21.600 | 67.782 | 69.300 | 100.000 | NA |
| Crosswind | numeric | 20271 | 0 | 70 | 0.000 | 0.212 | 0.000 | 6.600 | NA |
| Temperature | numeric | 20271 | 0 | 293 | 12.900 | 22.906 | 22.700 | 37.800 | NA |
| Barometric_Pressure | numeric | 20271 | 0 | 234 | 798.700 | 811.820 | 807.400 | 830.000 | NA |
| Headwind | numeric | 20271 | 0 | 102 | -5.300 | 0.061 | 0.000 | 4.600 | NA |
| Direction_True | numeric | 20271 | 0 | 489 | 0.000 | 208.080 | 237.000 | 359.000 | NA |
| Direction_Mag | numeric | 20271 | 0 | 490 | 0.000 | 208.120 | 238.000 | 359.000 | NA |
| Wind_Speed | numeric | 20271 | 0 | 77 | 0.000 | 0.323 | 0.300 | 7.200 | NA |
| Heat_Stress_Index | numeric | 20271 | 0 | 329 | 13.300 | 23.139 | 22.800 | 48.900 | NA |
| Altitude | numeric | 20271 | 0 | 317 | 1648.000 | 1828.400 | 1872.000 | 1960.000 | NA |
| Dew_Point | numeric | 20271 | 0 | 264 | 5.400 | 15.945 | 16.500 | 27.500 | NA |
| Density_Altitude | numeric | 20271 | 0 | 1023 | 2.105 | 2.601 | 2.596 | 3.152 | NA |
| Wind_Chill | numeric | 20271 | 0 | 300 | 12.800 | 22.863 | 22.700 | 37.800 | NA |
| Estado_Fenologico_Predominante | character | 20271 | 22 | 18 | NA | NA | NA | NA | 717 |
| Frutos_Afectados | numeric | 20271 | 0 | 29 | 0.000 | 0.141 | 0.000 | 60.000 | NA |
| Fecha | Date | 20271 | 6078 | 35 | NA | NA | NA | NA | 2019-08-28 |
| Hora | character | 20271 | 0 | 14093 | NA | NA | NA | NA | 10:11:12 |
La tabla integra estadística descriptiva de variables numéricas y, para variables no numéricas, clase, completitud, cardinalidad y categoría modal. Los campos de fecha y hora derivados se conservan para trazabilidad, sin sustituir la columna original.
3.2 Diccionario de datos
descripciones <- c(
id_arbol="Identificador del árbol", Latitude="Latitud geográfica", Longitude="Longitud geográfica",
FORMATTED_DATE_TIME="Fecha y hora original de medición", Psychro_Wet_Bulb_Temperature="Temperatura de bulbo húmedo psicrométrico",
Station_Pressure="Presión de estación", Relative_Humidity="Humedad relativa", Crosswind="Componente transversal del viento",
Temperature="Temperatura", Barometric_Pressure="Presión barométrica", Headwind="Componente frontal del viento",
Direction_True="Dirección respecto al norte verdadero", Direction_Mag="Dirección respecto al norte magnético",
Wind_Speed="Velocidad del viento", Heat_Stress_Index="Índice de estrés térmico", Altitude="Altitud",
Dew_Point="Punto de rocío", Density_Altitude="Altitud de densidad", Wind_Chill="Sensación térmica por viento",
Estado_Fenologico_Predominante="Código del estado fenológico predominante", Frutos_Afectados="Conteo de frutos afectados",
Fecha="Fecha derivada", Hora="Hora derivada")
unidades <- c(Latitude="grados",Longitude="grados",Direction_True="grados",Direction_Mag="grados",
Temperature="°C",Psychro_Wet_Bulb_Temperature="°C",Relative_Humidity="%",Wind_Speed="sin metadato",
Altitude="sin metadato",Station_Pressure="sin metadato",Barometric_Pressure="sin metadato")
diccionario <- tibble(Variable=names(datos), Tipo=vapply(datos, function(x) class(x)[1], character(1)),
Descripcion=unname(descripciones[names(datos)]), Unidad=unname(unidades[names(datos)])) %>%
mutate(Unidad=ifelse(is.na(Unidad),"No especificada",Unidad))
diccionario %>% knitr::kable(caption="Tabla 2. Diccionario de datos (descripciones inferidas del nombre del campo; unidades solo cuando son explícitas)") %>%
kableExtra::kable_styling(full_width=FALSE, position="center")| Variable | Tipo | Descripcion | Unidad |
|---|---|---|---|
| id_arbol | character | Identificador del árbol | No especificada |
| Latitude | numeric | Latitud geográfica | grados |
| Longitude | numeric | Longitud geográfica | grados |
| FORMATTED_DATE_TIME | character | Fecha y hora original de medición | No especificada |
| Psychro_Wet_Bulb_Temperature | numeric | Temperatura de bulbo húmedo psicrométrico | °C |
| Station_Pressure | numeric | Presión de estación | sin metadato |
| Relative_Humidity | numeric | Humedad relativa | % |
| Crosswind | numeric | Componente transversal del viento | No especificada |
| Temperature | numeric | Temperatura | °C |
| Barometric_Pressure | numeric | Presión barométrica | sin metadato |
| Headwind | numeric | Componente frontal del viento | No especificada |
| Direction_True | numeric | Dirección respecto al norte verdadero | grados |
| Direction_Mag | numeric | Dirección respecto al norte magnético | grados |
| Wind_Speed | numeric | Velocidad del viento | sin metadato |
| Heat_Stress_Index | numeric | Índice de estrés térmico | No especificada |
| Altitude | numeric | Altitud | sin metadato |
| Dew_Point | numeric | Punto de rocío | No especificada |
| Density_Altitude | numeric | Altitud de densidad | No especificada |
| Wind_Chill | numeric | Sensación térmica por viento | No especificada |
| Estado_Fenologico_Predominante | character | Código del estado fenológico predominante | No especificada |
| Frutos_Afectados | numeric | Conteo de frutos afectados | No especificada |
| Fecha | Date | Fecha derivada | No especificada |
| Hora | character | Hora derivada | No especificada |
3.3 Calidad del primer periodo
coord_validas <- with(datos_p1, is.finite(Latitude) & is.finite(Longitude) &
between(Latitude,-90,90) & between(Longitude,-180,180))
calidad <- tibble(
Indicador=c("Registros", "Filas exactamente duplicadas", "Identificadores duplicados",
"Coordenadas faltantes o fuera de rango", "Celdas faltantes", "Fechas no interpretables"),
Valor=c(nrow(datos_p1), sum(duplicated(datos_p1)), sum(duplicated(datos_p1$id_arbol)),
sum(!coord_validas), sum(is.na(datos_p1)), sum(is.na(datos_p1$Fecha))))
calidad %>% knitr::kable(caption="Tabla 3. Indicadores de calidad del primer periodo") %>%
kableExtra::kable_styling(full_width=FALSE, position="center")| Indicador | Valor |
|---|---|
| Registros | 427 |
| Filas exactamente duplicadas | 0 |
| Identificadores duplicados | 0 |
| Coordenadas faltantes o fuera de rango | 0 |
| Celdas faltantes | 0 |
| Fechas no interpretables | 0 |
resumen_p1 <- purrr::imap_dfr(datos_p1, ~resumen_variable(.x,.y))
resumen_p1 %>% knitr::kable(caption="Tabla 4. Resumen del primer periodo", digits=3) %>%
kableExtra::kable_styling(full_width=FALSE, position="center")| Variable | Clase | N | Faltantes | Unicos | Min | Media | Mediana | Max | Categoria_mas_frecuente |
|---|---|---|---|---|---|---|---|---|---|
| id_arbol | character | 427 | 0 | 427 | NA | NA | NA | NA | 1 |
| Latitude | numeric | 427 | 0 | 354 | 2.371 | 2.372 | 2.372 | 2.373 | NA |
| Longitude | numeric | 427 | 0 | 348 | -76.608 | -76.607 | -76.607 | -76.607 | NA |
| FORMATTED_DATE_TIME | character | 427 | 0 | 427 | NA | NA | NA | NA | 30/07/2019 1:00:56 p, m, |
| Psychro_Wet_Bulb_Temperature | numeric | 427 | 0 | 63 | 14.300 | 17.026 | 17.000 | 24.200 | NA |
| Station_Pressure | numeric | 427 | 0 | 36 | 802.100 | 803.580 | 803.500 | 805.600 | NA |
| Relative_Humidity | numeric | 427 | 0 | 146 | 36.400 | 43.841 | 43.200 | 65.200 | NA |
| Crosswind | numeric | 427 | 0 | 8 | 0.000 | 0.117 | 0.100 | 0.700 | NA |
| Temperature | numeric | 427 | 0 | 79 | 22.300 | 26.085 | 26.300 | 30.900 | NA |
| Barometric_Pressure | numeric | 427 | 0 | 22 | 802.000 | 803.540 | 803.500 | 805.500 | NA |
| Headwind | numeric | 427 | 0 | 23 | -2.600 | -0.440 | -0.400 | 2.100 | NA |
| Direction_True | numeric | 427 | 0 | 105 | 13.000 | 180.070 | 178.000 | 353.000 | NA |
| Direction_Mag | numeric | 427 | 0 | 106 | 12.000 | 179.570 | 178.000 | 353.000 | NA |
| Wind_Speed | numeric | 427 | 0 | 18 | 0.000 | 0.488 | 0.500 | 2.600 | NA |
| Heat_Stress_Index | numeric | 427 | 0 | 89 | 21.600 | 25.392 | 25.400 | 34.400 | NA |
| Altitude | numeric | 427 | 0 | 36 | 1891.000 | 1910.700 | 1911.000 | 1926.000 | NA |
| Dew_Point | numeric | 427 | 0 | 68 | 9.900 | 12.775 | 12.500 | 22.600 | NA |
| Density_Altitude | numeric | 427 | 0 | 198 | 2.655 | 2.790 | 2.797 | 2.960 | NA |
| Wind_Chill | numeric | 427 | 0 | 80 | 22.200 | 26.044 | 26.300 | 30.900 | NA |
| Estado_Fenologico_Predominante | character | 427 | 0 | 2 | NA | NA | NA | NA | 715 |
| Frutos_Afectados | numeric | 427 | 0 | 19 | 0.000 | 0.960 | 0.000 | 50.000 | NA |
| Fecha | Date | 427 | 0 | 1 | NA | NA | NA | NA | 2019-07-30 |
| Hora | character | 427 | 0 | 427 | NA | NA | NA | NA | 1:00:56 |
Para el análisis espacial se eliminan únicamente registros sin coordenadas válidas, valores no finitos de las variables candidatas y duplicados exactos de posición. No se imputan valores, para evitar introducir estructura espacial artificial.
4 Análisis geoestadístico
candidatas <- c("Temperature","Relative_Humidity","Wind_Speed","Altitude")
base_esp <- datos_p1 %>% filter(coord_validas) %>% distinct(Longitude,Latitude,.keep_all=TRUE)
puntos_ll <- st_as_sf(base_esp, coords=c("Longitude","Latitude"), crs=4326, remove=FALSE)
puntos <- st_transform(puntos_ll, 32618) # WGS84 / UTM zona 18N; distancias en metros
xy <- st_coordinates(puntos); puntos$x <- xy[,1]; puntos$y <- xy[,2]
# Vecindad k=8 estable para un patrón de muestreo irregular.
nb <- spdep::knn2nb(spdep::knearneigh(xy, k=min(8,nrow(puntos)-1)))
lw <- spdep::nb2listw(nb, style="W", zero.policy=TRUE)
pruebas <- purrr::map_dfr(candidatas, function(v){
z <- puntos[[v]]; ok <- is.finite(z)
if(!all(ok)) return(tibble(Variable=v, Moran_I=NA, p_permutacion=NA, N=sum(ok)))
mc <- spdep::moran.mc(z, lw, nsim=999, zero.policy=TRUE, alternative="greater")
tibble(Variable=v, Moran_I=unname(mc$statistic), p_permutacion=mc$p.value, N=length(z))
})
seleccionables <- pruebas %>% filter(!is.na(p_permutacion), p_permutacion < 0.05)
if(nrow(seleccionables)==0) stop("Ninguna variable candidata mostró autocorrelación positiva significativa (p<0.05); no es válido ejecutar kriging.")
variable_sel <- candidatas[candidatas %in% seleccionables$Variable][1]
pruebas %>% mutate(Seleccionada=Variable==variable_sel) %>%
knitr::kable(caption="Tabla 5. Moran global con 999 permutaciones (alternativa: autocorrelación positiva)",digits=4) %>%
kableExtra::kable_styling(full_width=FALSE,position="center")| Variable | Moran_I | p_permutacion | N | Seleccionada |
|---|---|---|---|---|
| Temperature | 0.6939 | 0.001 | 427 | TRUE |
| Relative_Humidity | 0.4385 | 0.001 | 427 | FALSE |
| Wind_Speed | 0.1454 | 0.001 | 427 | FALSE |
| Altitude | 0.9431 | 0.001 | 427 | FALSE |
La variable seleccionada por la regla predefinida es Temperature. La inferencia utiliza 999 permutaciones y nivel de significancia de 5 %. Moran se calcula sobre ocho vecinos cercanos; esta elección evita un umbral arbitrario y mantiene igual número nominal de vecinos por punto.
4.1 Tendencia espacial
form_trend <- as.formula(paste(variable_sel,"~ x + y"))
mod_trend <- lm(form_trend, data=st_drop_geometry(puntos))
p_tend <- anova(mod_trend)[["Pr(>F)"]][1:2]
hay_tendencia <- any(p_tend < 0.05, na.rm=TRUE)
tend_tab <- broom_like <- tibble(Termino=c("x (Este)","y (Norte)"), p_valor=p_tend)
tend_tab %>% knitr::kable(caption="Tabla 6. Prueba de tendencia lineal espacial",digits=5) %>%
kableExtra::kable_styling(full_width=FALSE,position="center")| Termino | p_valor |
|---|---|
| x (Este) | 0.0053 |
| y (Norte) | 0.0000 |
Se detectó al menos un componente lineal espacial significativo; se aplica kriging universal y el semivariograma se estima sobre residuos de tendencia.
4.2 Semivariograma, correlograma y modelos teóricos
formula_krig <- if(hay_tendencia) form_trend else as.formula(paste(variable_sel,"~ 1"))
vg_emp <- gstat::variogram(formula_krig, puntos, cutoff=max(dist(xy))/2, width=max(dist(xy))/20)
var_z <- var(puntos[[variable_sel]],na.rm=TRUE); rango0 <- max(vg_emp$dist)/3
modelos <- c(Exp="Exp", Gau="Gau", Sph="Sph")
ajustes <- purrr::imap(modelos, function(cod,nm){
ini <- gstat::vgm(psill=0.8*var_z, model=cod, range=rango0, nugget=0.2*var_z)
suppressWarnings(gstat::fit.variogram(vg_emp, ini, fit.method=6))
})
curvas <- purrr::imap_dfr(ajustes, function(m,nm){
pr <- variogramLine(m, maxdist=max(vg_emp$dist), n=200)
tibble(dist=pr$dist,gamma=pr$gamma,Modelo=nm)
})
ggplot() + geom_point(data=vg_emp,aes(dist,gamma),size=2) +
geom_line(data=curvas,aes(dist,gamma,color=Modelo),linewidth=1) +
labs(title=paste("Semivariograma experimental y modelos para",variable_sel),x="Distancia (m)",y="Semivarianza") + theme_minimal()
dmax <- max(dist(xy))/2; cortes <- seq(0,dmax,length.out=9)
correg <- purrr::map_dfr(seq_len(length(cortes)-1), function(i){
nbi <- spdep::dnearneigh(xy, cortes[i], cortes[i+1])
if(sum(card(nbi))==0) return(tibble(Distancia=mean(cortes[i:(i+1)]),I=NA_real_,p=NA_real_))
lwi <- nb2listw(nbi,style="W",zero.policy=TRUE)
mt <- moran.mc(puntos[[variable_sel]],lwi,nsim=499,zero.policy=TRUE,alternative="two.sided")
tibble(Distancia=mean(cortes[i:(i+1)]),I=unname(mt$statistic),p=mt$p.value)
})
ggplot(correg,aes(Distancia,I)) + geom_hline(yintercept=0,lty=2) + geom_line() + geom_point(aes(color=p<0.05),size=2) +
scale_color_manual(values=c("grey40","#D55E00"),name="p < 0,05") +
labs(title=paste("Correlograma espacial de",variable_sel),x="Distancia media de clase (m)",y="I de Moran") + theme_minimal()
4.3 Validación cruzada y selección
cv <- purrr::imap_dfr(ajustes, function(m,nm){
out <- gstat::krige.cv(formula_krig, puntos, model=m, nfold=nrow(puntos), verbose=FALSE)
tibble(Modelo=nm, ME=mean(out$residual,na.rm=TRUE),
MAE=mean(abs(out$residual),na.rm=TRUE), RMSE=sqrt(mean(out$residual^2,na.rm=TRUE)),
RMSSE=sqrt(mean((out$residual/sqrt(out$var1.var))^2,na.rm=TRUE)))
})
mejor <- cv$Modelo[which.min(cv$RMSE)]; modelo_final <- ajustes[[mejor]]
cv %>% mutate(Seleccionado=Modelo==mejor) %>% knitr::kable(caption="Tabla 7. Validación cruzada dejar-uno-fuera",digits=4) %>%
kableExtra::kable_styling(full_width=FALSE,position="center")| Modelo | ME | MAE | RMSE | RMSSE | Seleccionado |
|---|---|---|---|---|---|
| Exp | 0.0045 | 0.7285 | 0.9482 | 0.9629 | TRUE |
| Gau | 0.0013 | 0.8668 | 1.1073 | 0.9396 | FALSE |
| Sph | 0.0031 | 0.7807 | 1.0082 | 0.9305 | FALSE |
Se selecciona Exp por el menor RMSE. Como diagnósticos complementarios, ME próximo a cero indica bajo sesgo y RMSSE próximo a uno indica que la varianza de kriging está razonablemente calibrada [2], [3].
5 Predicción espacial por kriging
hull <- st_convex_hull(st_union(puntos))
grid <- st_make_grid(
hull,
n = c(120,120),
what = "centers"
)
grid <- st_sf(
geometry = grid
)
grid <- st_filter(grid, hull)
gxy <- st_coordinates(grid)
grid$x <- gxy[,1]
grid$y <- gxy[,2]
pred <- gstat::krige(formula_krig, puntos, newdata=grid, model=modelo_final)## [using universal kriging]p1 <- ggplot() + geom_sf(data=pred,aes(color=var1.pred),size=1) +
geom_sf(data=puntos,shape=21,fill="white",size=1) + scale_color_viridis_c(name=variable_sel) +
labs(title=paste("Predicción por kriging -",variable_sel),subtitle=paste("Modelo",mejor)) + theme_minimal()
p2 <- ggplot() + geom_sf(data=pred,aes(color=sqrt(var1.var)),size=1) +
geom_sf(data=puntos,shape=21,fill="white",size=1) + scale_color_viridis_c(name="Error estándar") +
labs(title="Incertidumbre de predicción",subtitle="Raíz de la varianza de kriging") + theme_minimal()
print(p1); print(p2)
La primera superficie representa la predicción dentro de la envolvente convexa de los puntos; no se extrapola fuera del dominio muestreado. El mapa de error estándar debe interpretarse junto con la predicción: zonas más alejadas o menos respaldadas por el diseño de muestreo suelen exhibir mayor incertidumbre.
6 Conclusiones
- El primer periodo se identificó objetivamente como 30/07/2019, con 427 registros antes de la depuración espacial.
- La variable finalmente modelada fue Temperature, porque fue la primera candidata en el orden solicitado con Moran positivo significativo a 5 % (I = 0.694; p = 0.0010).
- La evidencia de tendencia justificó kriging universal.
- El modelo Exp obtuvo el menor RMSE de validación cruzada (0.948), por lo que sustenta la superficie final. La incertidumbre cartografiada limita la interpretación y evita tratar la predicción como observación directa.
- Las unidades no explícitas en la base no se presumen. Deben confirmarse con el responsable de captura antes de decisiones agronómicas cuantitativas.
7 Referencias
[1] P. A. P. Moran, “Notes on continuous stochastic phenomena,” Biometrika, vol. 37, no. 1/2, pp. 17–23, 1950, doi: 10.2307/2332142.
[2] N. Cressie, Statistics for Spatial Data, rev. ed. New York, NY, USA: Wiley, 1993, doi: 10.1002/9781119115151.
[3] H. Wackernagel, Multivariate Geostatistics: An Introduction with Applications, 3rd ed. Berlin, Germany: Springer, 2003, doi: 10.1007/978-3-662-05294-5.
[4] E. J. Pebesma, “Multivariable geostatistics in S: the gstat package,” Computers & Geosciences, vol. 30, no. 7, pp. 683–691, 2004, doi: 10.1016/j.cageo.2004.03.012.
[5] R. S. Bivand and D. W. S. Wong, “Comparing implementations of global and local indicators of spatial association,” TEST, vol. 27, pp. 716–748, 2018, doi: 10.1007/s11749-018-0599-x.
8 Anexos
8.1 Anexo A. Parámetros del variograma seleccionado
modelo_final %>% as.data.frame() %>% knitr::kable(caption=paste("Parámetros del modelo",mejor),digits=4) %>%
kableExtra::kable_styling(full_width=FALSE,position="center")| model | psill | range | kappa | ang1 | ang2 | ang3 | anis1 | anis2 |
|---|---|---|---|---|---|---|---|---|
| Nug | 0.0000 | 0.000 | 0.0 | 0 | 0 | 0 | 1 | 1 |
| Exp | 3.2973 | 18.054 | 0.5 | 0 | 0 | 0 | 1 | 1 |
8.2 Anexo B. Reproducibilidad
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tibble_3.3.1 purrr_1.2.1 kableExtra_1.4.1 knitr_1.52
## [5] ggplot2_4.0.3 gstat_2.1-6 spdep_1.4-2 spData_2.3.5
## [9] sf_1.1-3 lubridate_1.9.5 stringr_1.6.0 tidyr_1.3.2
## [13] dplyr_1.2.0 readxl_1.4.5
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.56 bslib_0.10.0 lattice_0.22-6
## [5] vctrs_0.7.1 tools_4.5.0 generics_0.1.4 parallel_4.5.0
## [9] proxy_0.4-29 spacetime_1.3-4 xts_0.14.1 pkgconfig_2.0.3
## [13] KernSmooth_2.23-26 RColorBrewer_1.1-3 S7_0.2.1 lifecycle_1.0.5
## [17] compiler_4.5.0 farver_2.1.2 deldir_2.0-4 FNN_1.1.4.1
## [21] textshaping_1.0.4 stars_0.7-3 htmltools_0.5.9 class_7.3-23
## [25] sass_0.4.10 yaml_2.3.12 pillar_1.11.1 jquerylib_0.1.4
## [29] classInt_0.4-11 cachem_1.1.0 wk_0.9.5 abind_1.4-8
## [33] boot_1.3-31 tidyselect_1.2.1 digest_0.6.39 stringi_1.8.7
## [37] labeling_0.4.3 fastmap_1.2.0 grid_4.5.0 cli_3.6.5
## [41] magrittr_2.0.4 e1071_1.7-17 withr_3.0.2 scales_1.4.0
## [45] sp_2.2-0 timechange_0.4.0 rmarkdown_2.32 igraph_2.3.3
## [49] otel_0.2.0 cellranger_1.1.0 zoo_1.8-15 evaluate_1.0.5
## [53] viridisLite_0.4.3 s2_1.1.11 rlang_1.1.7 Rcpp_1.1.1
## [57] glue_1.8.0 DBI_1.3.0 xml2_1.5.2 svglite_2.2.2
## [61] rstudioapi_0.18.0 jsonlite_2.0.0 R6_2.6.1 systemfonts_1.3.1
## [65] intervals_0.15.5 units_1.0-18.3 Anexo C. Control de cambios y limitaciones
Versión 1.0.0. Se incorpora carga validada, extracción reproducible del primer periodo, diagnóstico de calidad, selección jerárquica de variable, tendencia, Moran por permutación, correlograma, comparación conjunta de tres variogramas, validación cruzada, kriging y mapas de predicción/incertidumbre. No se elimina ninguna columna del archivo de origen; la base espacial es un objeto derivado. Dependencias y sesión quedan registradas. La ejecución está diseñada para RStudio 2025.05.0 Build 496, pero depende de que los paquetes listados y Pandoc estén instalados. La base no contiene metadatos completos de unidades ni límites oficiales del lote; por ello, el diccionario marca esas unidades como no especificadas y la predicción se restringe a la envolvente convexa del muestreo.