1 Introducción

En este cuaderno se ilustran dos técnicas de interpolación espacial para estimar el carbono orgánico del suelo (SOC) a una profundidad de 15-30 cm en el departamento de Sucre, Colombia. Las técnicas aplicadas son:

IDW (Inverse Distance Weighted), un método determinístico que asigna mayor peso a las muestras cercanas al punto a predecir.

OK (Ordinary Kriging), un método geoestadístico probabilístico que considera la autocorrelación espacial mediante un modelo de variograma.

Ambos métodos permiten generar una superficie continua de SOC a partir de datos de muestreo obtenidos de la base de datos SoilGrids 250 m.

2 Limpieza de memoria e instalación de librerías

Primero, se limpia la memoria del entorno de R con rm(list = ls()) para asegurar que no haya objetos residuales de sesiones anteriores que puedan interferir con la ejecución del código. Luego, se verifica la instalación de las librerías necesarias y se procede a cargarlas.La librería curl, que se usa para interactuar con URLs y realizar solicitudes HTTP

rm(list=ls())
library(sp)
library(terra)
## terra 1.8.54
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(curl)
## Using libcurl 8.14.1 with Schannel

h <- new_handle() es el paso inicial para configurar de manera detallada cómo R interactuará con un recurso web, permitiendo un control preciso sobre el comportamiento de las solicitudes de red.

h <- new_handle()
handle_setopt(h, http_version = 2)

3 lectura del archivo raster

Necesitamos leer un conjunto de datos para simular datos del mundo real. Para ello, leemos la capa de COS descargada de ISRIC utilizando la librería terra. El resultado muestra las características del objeto SpatRaster cargado, incluyendo sus dimensiones, resolución, extensión, sistema de referencia de coordenadas (CRS) y el nombre de la capa.

archivo <- ("soc_igh_15_30_sucre.tif")
(soc <- rast(archivo))
## class       : SpatRaster 
## size        : 832, 553, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8466250, -8328000, 921500, 1129500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30_sucre.tif 
## name        : soc_igh_15_30_sucre

Ahora, convertimos los datos de COS a porcentaje. El factor de escala para el COS en el sitio web de SoilGrids es 10, donde los datos se proporcionan en dg/kg (decigramos por kilogramo). Para convertir a porcentaje (g/100g), se debe dividir por 10.

soc.perc <-  soc/10

El CRS (Sistema de Referencia de Coordenadas) de los datos del mundo real es Interrupted_Goode_Homolosine. Parece que necesitamos una transformación de este CRS al conocido CRS WGS84, que es ampliamente utilizado para datos geográficos. La transformación produce un nuevo objeto SpatRaster con las siguientes características:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## size        : 856, 636, 1  (nrow, ncol, nlyr)
## resolution  : 0.002184037, 0.002184037  (x, y)
## extent      : -75.8015, -74.41245, 8.276935, 10.14647  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30_sucre 
## min value   :            8.534835 
## max value   :          144.756393

A continuación, convertimos la capa SpatRaster en un objeto stars, un formato de datos espacial que facilita el trabajo con arreglos espacio-temporales.

stars.soc = st_as_stars(geog.soc)

4 Muestreo del mapa

Para simular un escenario de muestreo, obtenemos una muestra de aproximadamente 500 sitios de los datos del mundo real utilizando un muestreo aleatorio. Las principales características del objeto samples son:

Clase (class): Es un objeto SpatVector de tipo points, lo que indica que contiene puntos geográficos.

Dimensiones (dimensions): Contiene 500 geometrías (puntos) y 1 atributo asociado.

Extensión (extent): Define el área geográfica cubierta por los puntos de la muestra, en coordenadas de longitud y latitud (xmin, xmax, ymin, ymax).

CRS (coord. ref.): El sistema de referencia de coordenadas es WGS84 (+proj=longlat +datum=WGS84 +no_defs), consistente con la transformación realizada previamente.

Nombre (names): La capa se llama soc_igh_15_30.

Tipo (type) y Valores (values): Contiene valores numéricos para el atributo principal (COS), con ejemplos como 24.8, 45.83, 45.85, etc.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) 
#
m  # Print the map

Ahora, convertimos el objeto SpatVector a un objeto simple feature (sf), un formato estándar y moderno para datos espaciales en R, que permite una manipulación y análisis más flexibles. Los datos muestran las columnas soc_igh_15_30_sucre (el valor de COS) y geometry (las coordenadas del punto).

set.seed(123456)

# Random sampling of 500 points
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -75.80041, -74.41354, 8.278027, 10.1432  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30_sucre
##  type        :               <num>
##  values      :               16.77
##                               15.5
##                              14.66
(muestras <- sf::st_as_sf(samples))
## Simple feature collection with 500 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -75.80041 ymin: 8.278027 xmax: -74.41354 ymax: 10.1432
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##    soc_igh_15_30_sucre                   geometry
## 1             16.77218 POINT (-75.23692 10.01652)
## 2             15.50385  POINT (-75.1714 10.05583)
## 3             14.66136 POINT (-74.96829 9.872374)
## 4                   NA  POINT (-74.7346 9.629946)
## 5             30.02661  POINT (-74.7608 9.649603)
## 6             14.23169 POINT (-74.62103 9.985944)
## 7             30.80341 POINT (-74.59482 9.308893)
## 8             21.73104  POINT (-75.00978 9.38315)
## 9             31.91427 POINT (-74.57079 8.699546)
## 10                  NA POINT (-74.45722 8.489879)

Luego, se eliminan los valores NA (valores no disponibles) del conjunto de datos muestras. Esto es crucial para asegurar que los análisis posteriores no se vean afectados por datos faltantes. El resultado muestra las filas restantes después de la eliminación de los NA.

nmuestras <- na.omit(muestras)
longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$soc_igh_15_30_sucre
id <- seq(1,500,1) 
(sitios <- data.frame(id, longit, latit, soc))
##      id    longit     latit       soc
## 1     1 -75.23692 10.016521  16.77218
## 2     2 -75.17140 10.055834  15.50385
## 3     3 -74.96829  9.872374  14.66136
## 4     4 -74.73460  9.629946        NA
## 5     5 -74.76080  9.649603  30.02661
## 6     6 -74.62103  9.985944  14.23169
## 7     7 -74.59482  9.308893  30.80341
## 8     8 -75.00978  9.383150  21.73104
## 9     9 -74.57079  8.699546  31.91427
## 10   10 -74.45722  8.489879        NA
## 11   11 -75.17359  9.785013  20.12842
## 12   12 -75.16704 10.143195  15.08282
## 13   13 -74.84380  8.284579  35.68023
## 14   14 -75.50556  9.745700  23.96885
## 15   15 -75.40291  8.751963  22.81652
## 16   16 -75.35486  8.749779  23.10038
## 17   17 -75.33957  8.756331  19.22761
## 18   18 -75.53177  9.837430  13.75122
## 19   19 -75.04255  8.393781        NA
## 20   20 -75.15830  8.478959        NA
## 21   21 -74.81759  9.398438  18.98332
## 22   22 -74.62976  8.317340  41.68094
## 23   23 -75.48154  9.204059  22.26417
## 24   24 -75.77420  9.767541        NA
## 25   25 -75.04691  9.199691  19.90105
## 26   26 -74.43538  8.778172        NA
## 27   27 -74.52711  9.243372  34.58039
## 28   28 -74.45722  9.883295  12.88040
## 29   29 -75.29153  8.924502  34.81120
## 30   30 -75.17140  8.413437  24.05623
## 31   31 -74.77172  9.129802        NA
## 32   32 -75.53177  9.719492  31.86691
## 33   33 -75.06439  9.536033  12.51832
## 34   34 -75.47498  9.153826  16.69121
## 35   35 -75.66063  9.697651        NA
## 36   36 -75.18014  9.972840  17.49859
## 37   37 -74.83506  9.000943        NA
## 38   38 -75.34176  9.112329  22.86924
## 39   39 -75.49027  9.946632  17.84215
## 40   40 -74.82196  9.387518  30.23796
## 41   41 -74.90932  9.374414  18.97344
## 42   42 -75.16485 10.051466  14.94454
## 43   43 -75.46843  9.470512        NA
## 44   44 -74.74333 10.123539        NA
## 45   45 -74.74552  9.297973  19.54498
## 46   46 -74.82633  8.972551        NA
## 47   47 -74.83506  9.782829  24.24742
## 48   48 -75.35923  9.389702  20.22128
## 49   49 -75.57982  9.688915  30.29294
## 50   50 -75.74581  9.850534        NA
## 51   51 -75.26095  9.990312  16.44066
## 52   52 -75.31555  8.500799  34.00455
## 53   53 -74.53148  8.721387  46.69879
## 54   54 -74.79793  8.537928  35.63427
## 55   55 -75.36360  8.610001  31.92605
## 56   56 -75.29589  8.848061  16.87762
## 57   57 -75.44004  9.007496  21.93176
## 58   58 -75.09715  8.485511  42.53857
## 59   59 -74.83943  9.182219  12.78795
## 60   60 -75.68465  9.946632        NA
## 61   61 -74.88529  8.730123  37.32200
## 62   62 -75.06220  9.975024  19.40531
## 63   63 -75.09715  9.804669  18.49838
## 64   64 -74.44194  9.754436        NA
## 65   65 -75.40291  9.297973        NA
## 66   66 -74.99668  8.494247  31.72160
## 67   67 -74.56861  8.284579  64.65160
## 68   68 -75.29589  8.618737  15.97733
## 69   69 -75.56453  9.533849  36.80084
## 70   70 -75.52959 10.084226 113.31328
## 71   71 -74.97484  9.780645  14.21024
## 72   72 -74.71057  8.472406  26.48750
## 73   73 -75.61258  9.048992  39.49779
## 74   74 -74.97921  9.055544  21.09655
## 75   75 -75.67592  9.365678  35.09657
## 76   76 -75.28497  9.291420  18.10870
## 77   77 -74.81104  9.402806  19.02714
## 78   78 -75.66936  8.880821  31.43431
## 79   79 -74.85472  8.658050  33.73241
## 80   80 -74.76299  9.741332        NA
## 81   81 -74.42228  9.697651        NA
## 82   82 -74.72149  9.442119  26.03761
## 83   83 -74.91587  8.553216  33.98458
## 84   84 -75.38762  8.437462  25.15363
## 85   85 -75.22819  8.834957  31.09999
## 86   86 -74.83069  8.317340  31.44515
## 87   87 -75.48154  9.350389  27.78628
## 88   88 -75.71305  8.596897  14.50197
## 89   89 -75.34613  9.315445  12.53536
## 90   90 -74.87656 10.134459  12.78898
## 91   91 -75.77638  8.806564        NA
## 92   92 -74.66471  9.658339  23.64197
## 93   93 -75.53177  9.308893  23.51379
## 94   94 -74.81104  8.937606  57.20580
## 95   95 -75.25877  9.527297  14.79755
## 96   96 -74.88529  8.723571  34.17416
## 97   97 -74.71931  9.481432  25.10708
## 98   98 -75.39199  9.507640  17.60667
## 99   99 -74.51838 10.014337  15.42885
## 100 100 -75.71741  8.824036  25.74006
## 101 101 -75.34394  8.845877  22.07300
## 102 102 -75.46625  8.992207  18.01147
## 103 103 -75.56016  8.422174  35.77085
## 104 104 -75.64316  9.374414  36.22196
## 105 105 -74.51401  9.348205  15.36568
## 106 106 -75.22600  9.490168  15.27421
## 107 107 -75.12335  8.352284        NA
## 108 108 -75.06002  8.850245  37.56453
## 109 109 -74.55551  9.533849  16.34492
## 110 110 -74.75207  9.857086  13.26641
## 111 111 -75.53832 10.082042  70.02888
## 112 112 -75.25658  9.291420  15.39138
## 113 113 -75.15612  9.913871  18.75975
## 114 114 -75.04910  9.542585  13.24178
## 115 115 -75.12117  8.319524        NA
## 116 116 -75.26313  8.813116  16.14482
## 117 117 -75.45970  9.092673  34.18155
## 118 118 -74.56424  8.308604  54.10000
## 119 119 -75.03599  8.570688  33.85738
## 120 120 -74.80230  9.151642  23.74385
## 121 121 -74.96829  8.920134  32.91829
## 122 122 -75.06220  8.516087  35.27352
## 123 123 -74.66471  9.754436  12.57500
## 124 124 -75.52303  9.555689  36.26036
## 125 125 -74.60574  8.920134  46.31846
## 126 126 -74.46378  9.795933  16.83923
## 127 127 -75.01634  9.311077  18.44295
## 128 128 -75.47935  9.507640  25.97346
## 129 129 -75.36360  8.896110  23.32726
## 130 130 -75.12772  8.507351  43.59589
## 131 131 -74.70402  8.889558  31.04134
## 132 132 -75.60603  9.411543  37.46459
## 133 133 -74.51838  9.730412  16.87222
## 134 134 -75.33302  8.786908  24.70642
## 135 135 -75.48591  9.741332  32.62310
## 136 136 -75.18014  9.726044  23.12026
## 137 137 -74.52930  8.815300  40.97387
## 138 138 -75.79604  8.433094        NA
## 139 139 -74.61884  8.518271  30.86839
## 140 140 -74.80449  8.965999        NA
## 141 141 -74.83943  9.643051  18.00583
## 142 142 -75.28279  9.256476  15.97415
## 143 143 -75.72397  8.588160  14.04438
## 144 144 -75.73925  9.660523        NA
## 145 145 -75.36578  8.802196  19.06566
## 146 146 -75.63005  8.304236  14.42288
## 147 147 -75.30245  8.280211  20.82391
## 148 148 -75.12117  9.077385  17.57601
## 149 149 -74.60355  9.769725  15.83266
## 150 150 -74.47906 10.020889  15.65074
## 151 151 -75.30681  9.313261  37.00136
## 152 152 -74.85035  8.712651  53.81565
## 153 153 -74.72586  8.784724  59.22336
## 154 154 -75.38107  8.815300  18.36914
## 155 155 -75.27842  9.763173  19.03321
## 156 156 -75.70431  8.730123  18.74288
## 157 157 -74.65160  9.470512  14.15269
## 158 158 -75.64097  9.169114  41.97935
## 159 159 -74.56206  8.365389  45.79074
## 160 160 -75.05783 10.127907  14.60997
## 161 161 -75.44659  8.441830  23.83446
## 162 162 -74.93771  9.942264  16.10072
## 163 163 -75.80041  8.668970        NA
## 164 164 -74.83943  8.627473  34.80809
## 165 165 -75.14083  9.970656  19.44724
## 166 166 -74.66689  8.856797        NA
## 167 167 -74.48343  8.948527        NA
## 168 168 -74.78920  9.636498  39.91014
## 169 169 -74.91806  9.376598  27.66999
## 170 170 -75.30463  9.201875  23.12178
## 171 171 -75.72397  8.869901  60.28939
## 172 172 -75.30681  8.629657  12.54826
## 173 173 -75.49901  8.354468  24.78473
## 174 174 -75.29808  8.391597  34.74736
## 175 175 -74.68655  9.514192  22.28006
## 176 176 -74.61447  9.254292        NA
## 177 177 -74.86564  9.234636  16.25056
## 178 178 -74.75644  9.977208  11.42093
## 179 179 -75.25003 10.029625  15.51335
## 180 180 -74.79138  9.350389  50.34549
## 181 181 -74.45722  8.647129        NA
## 182 182 -75.50775  9.446487  32.71515
## 183 183 -74.62539 10.073306  14.50429
## 184 184 -75.42038  9.354758  30.93233
## 185 185 -74.80012  9.223715  16.88898
## 186 186 -74.87001 10.071122        NA
## 187 187 -75.79385 10.132275        NA
## 188 188 -75.05128  9.055544  15.07213
## 189 189 -74.81104  8.419990  29.10621
## 190 190 -75.40073  8.409069  12.06130
## 191 191 -74.51619  9.975024  12.59406
## 192 192 -74.54677  8.832773  41.28361
## 193 193 -74.59045  9.308893  28.38046
## 194 194 -74.48343  9.180035        NA
## 195 195 -74.61884  8.380677  40.94198
## 196 196 -75.23474  9.304525  14.56315
## 197 197 -74.68873 10.029625  14.45933
## 198 198 -74.66252  9.758804  16.97757
## 199 199 -75.33302  9.776277  17.05967
## 200 200 -74.90277  8.350100  27.11870
## 201 201 -75.22164  8.417806  36.93121
## 202 202 -75.70649  9.536033        NA
## 203 203 -74.41354  8.876453        NA
## 204 204 -75.52085  9.064281  34.86385
## 205 205 -75.72833  8.706098  21.22930
## 206 206 -74.43538  8.762883        NA
## 207 207 -74.90932  9.387518  18.76190
## 208 208 -75.24348  8.963815  16.40841
## 209 209 -74.88966  8.284579  33.24475
## 210 210 -74.54240  9.042440  45.77991
## 211 211 -75.67155 10.018705        NA
## 212 212 -75.01634  9.813405  15.16873
## 213 213 -75.18888  8.732307  25.65425
## 214 214 -75.11243  8.957263  23.43262
## 215 215 -74.41354  9.581897        NA
## 216 216 -75.13646  8.935422  28.56810
## 217 217 -75.31773  8.367573  50.86503
## 218 218 -75.73052  8.500799  14.48276
## 219 219 -74.68873  8.813116  45.38155
## 220 220 -75.63005 10.138827        NA
## 221 221 -75.71741  8.771620  29.81711
## 222 222 -75.77638  8.865533        NA
## 223 223 -75.40510  9.999049  11.65598
## 224 224 -74.51182  8.452750        NA
## 225 225 -75.18888  9.415911  13.32109
## 226 226 -75.59948  9.289236  20.59500
## 227 227 -74.77828  9.975024        NA
## 228 228 -75.65408  9.876743        NA
## 229 229 -75.39199  8.834957  20.29601
## 230 230 -75.05347  9.249924  15.76778
## 231 231 -74.62976  8.939790  47.46514
## 232 232 -74.67781  9.658339        NA
## 233 233 -75.68684 10.121355        NA
## 234 234 -74.64287  9.522928  15.55548
## 235 235 -74.43757  9.396254        NA
## 236 236 -74.61666  8.478959  28.96349
## 237 237 -75.36797  9.708572  15.04645
## 238 238 -75.46188  9.734780  17.14963
## 239 239 -74.73023  9.553505  21.78990
## 240 240 -74.69528  9.483616  11.86203
## 241 241 -75.66063  9.564425        NA
## 242 242 -74.74115  8.885190  66.48537
## 243 243 -75.71086  8.278027  33.22520
## 244 244 -75.17796  8.885190  28.18065
## 245 245 -75.70431  9.826510        NA
## 246 246 -75.27187  9.127618  20.86129
## 247 247 -75.70431  9.020600  52.72978
## 248 248 -74.71494  9.110145        NA
## 249 249 -75.43567  8.658050  31.40872
## 250 250 -74.63195  9.975024  12.63383
## 251 251 -75.17796  9.625578  18.50673
## 252 252 -75.06875  8.601265  44.07619
## 253 253 -74.83288  8.566320  42.46313
## 254 254 -75.76983  9.123250        NA
## 255 255 -75.52959  8.852429  22.50928
## 256 256 -75.56671  9.570977  30.20056
## 257 257 -75.49246  8.730123  17.18093
## 258 258 -74.50527  9.647419  15.61578
## 259 259 -75.10370  8.913582  25.00490
## 260 260 -74.57516 10.005601  14.00175
## 261 261 -74.64068  8.409069  22.16337
## 262 262 -75.21945  8.511719        NA
## 263 263 -74.96174  9.937896  14.45909
## 264 264 -75.65844  8.939790  52.47978
## 265 265 -74.41573  8.874269        NA
## 266 266 -75.33739  9.324181  19.20984
## 267 267 -74.72586 10.058018  13.08686
## 268 268 -74.45941  8.760699        NA
## 269 269 -74.58171  9.750068  15.07444
## 270 270 -74.96610  8.527007  50.28327
## 271 271 -74.57516  9.346021        NA
## 272 272 -75.01197  8.470222  38.59378
## 273 273 -75.71305  8.963815  39.85454
## 274 274 -74.68873  8.710467  41.67802
## 275 275 -75.55798  9.343837  25.42660
## 276 276 -74.45941  9.073017        NA
## 277 277 -74.43757  9.767541        NA
## 278 278 -74.59482  9.540401  20.59985
## 279 279 -75.59074  9.603738        NA
## 280 280 -74.92898  9.629946  19.89052
## 281 281 -75.43349  8.406885  13.69670
## 282 282 -74.60355  9.180035        NA
## 283 283 -75.66936  9.852718        NA
## 284 284 -74.91587 10.121355  11.25775
## 285 285 -75.10370 10.106066  13.91882
## 286 286 -75.21508  8.334812  20.30000
## 287 287 -75.09059 10.092962  14.06280
## 288 288 -75.69557  9.819958        NA
## 289 289 -75.07312  9.411543  19.30968
## 290 290 -75.74144  9.651787        NA
## 291 291 -74.90714  9.306709  29.01457
## 292 292 -75.15393  9.809037  18.40184
## 293 293 -75.04691 10.134459  16.48414
## 294 294 -75.18451  9.656155  20.05356
## 295 295 -74.67563  9.153826        NA
## 296 296 -75.10151  8.717019        NA
## 297 297 -75.12991  9.459591  19.12000
## 298 298 -75.45314 10.090778 108.19901
## 299 299 -74.58171  9.175666  45.68884
## 300 300 -75.14520  9.007496  17.49349
## 301 301 -74.76736  9.540401  18.32735
## 302 302 -75.73707  8.998759  41.85026
## 303 303 -75.21072  8.963815  19.11983
## 304 304 -74.63850  9.014048  47.61535
## 305 305 -75.46406  8.636209  19.56232
## 306 306 -75.22600  8.577240  23.62259
## 307 307 -75.70431  8.898294  41.71069
## 308 308 -75.74144  9.671443        NA
## 309 309 -75.36797  9.629946  19.42245
## 310 310 -75.18232  9.162562  22.83064
## 311 311 -75.52085  8.675522  19.57948
## 312 312 -74.96392  8.917950  20.82520
## 313 313 -75.59292  8.500799  19.55676
## 314 314 -75.44222  9.514192  18.43251
## 315 315 -74.50964  8.424358        NA
## 316 316 -74.95518  8.310788  35.59274
## 317 317 -75.54706 10.025257  65.08063
## 318 318 -74.94208 10.020889  14.05648
## 319 319 -75.16704  8.907030  23.44567
## 320 320 -75.04691  9.835246  17.18410
## 321 321 -75.66718  8.616553  13.96818
## 322 322 -74.97266  8.627473  43.63337
## 323 323 -74.92461  9.507640  15.68844
## 324 324 -75.36142  8.315156        NA
## 325 325 -74.92898  8.378493  30.57215
## 326 326 -75.21727  9.450855  25.40955
## 327 327 -74.84817  9.754436  12.83601
## 328 328 -75.15612  9.892031  16.81657
## 329 329 -75.35705  9.431199  19.04481
## 330 330 -74.53585  9.940080  12.50328
## 331 331 -74.62976  9.385334  24.50165
## 332 332 -75.04255  9.540401  14.25257
## 333 333 -74.62976  9.101409        NA
## 334 334 -74.50309  8.664602        NA
## 335 335 -75.05347  9.975024  12.78016
## 336 336 -75.56890  9.018416  35.09132
## 337 337 -75.57545  9.503272  24.92367
## 338 338 -74.92461  9.284868  19.56976
## 339 339 -74.97047  9.459591  19.49793
## 340 340 -75.16267 10.009969  15.96885
## 341 341 -74.99231  9.575345  13.83771
## 342 342 -74.92242  8.861165        NA
## 343 343 -74.72804  8.782540        NA
## 344 344 -75.65189 10.060202        NA
## 345 345 -74.61011  8.516087  36.09144
## 346 346 -75.48372  9.319813  17.80127
## 347 347 -75.58855  8.282395  22.26337
## 348 348 -74.54240  8.376309  27.07471
## 349 349 -75.10370  9.512008  15.47777
## 350 350 -75.53395  9.747884  24.04121
## 351 351 -75.73270  9.177851        NA
## 352 352 -75.55143  8.439646  15.45402
## 353 353 -75.67155  9.249924  28.38600
## 354 354 -75.23474  9.247740  23.45555
## 355 355 -75.18014  9.173482  20.41983
## 356 356 -75.79385  8.944159        NA
## 357 357 -75.49027  9.415911  19.17684
## 358 358 -75.61040  9.920423        NA
## 359 359 -74.67781  9.190955        NA
## 360 360 -74.93553 10.121355  38.27866
## 361 361 -74.77609  9.350389        NA
## 362 362 -75.40073 10.123539  22.82617
## 363 363 -75.62787  8.802196  23.95023
## 364 364 -74.84161  9.088305        NA
## 365 365 -75.26313  9.219347  26.48905
## 366 366 -74.47688  8.714835        NA
## 367 367 -75.56016  8.778172  15.70485
## 368 368 -75.50338  8.968183  46.60524
## 369 369 -75.55143 10.077674  72.41357
## 370 370 -74.83506  8.703914  34.79946
## 371 371 -75.76983  8.509535  15.41593
## 372 372 -75.42694  8.880821  24.82403
## 373 373 -74.66471  9.608106        NA
## 374 374 -74.69092  9.643051  21.72633
## 375 375 -75.09496  9.197507  15.07751
## 376 376 -75.18888  9.236820  14.06141
## 377 377 -74.44630  9.581897        NA
## 378 378 -75.79167  9.588450        NA
## 379 379 -75.53395  9.394070  29.11247
## 380 380 -75.74362  9.505456        NA
## 381 381 -75.72178  8.291131  16.68115
## 382 382 -75.52522  9.616842  47.29915
## 383 383 -75.21290  9.092673  22.19215
## 384 384 -74.43757  8.485511        NA
## 385 385 -74.90058  9.579713  14.47940
## 386 386 -75.79604  9.889847        NA
## 387 387 -75.47498  8.712651  21.30603
## 388 388 -74.45504  8.789092        NA
## 389 389 -74.63631  9.588450  17.84518
## 390 390 -74.63850  9.453039  24.12402
## 391 391 -75.66936  9.027152  43.64952
## 392 392 -75.07531 10.132275  19.92230
## 393 393 -74.98576  8.765067        NA
## 394 394 -74.67126  8.990023        NA
## 395 395 -74.66252  8.858981  47.86135
## 396 396 -75.75454  8.686442  41.32888
## 397 397 -75.04036  8.719203  28.48967
## 398 398 -75.00323 10.103882  18.04936
## 399 399 -75.02289  8.880821        NA
## 400 400 -74.69965  8.876453  33.98787
## 401 401 -75.60821  9.649603        NA
## 402 402 -74.57298  8.382861  60.59375
## 403 403 -75.16048  8.507351        NA
## 404 404 -74.60355  9.289236  35.51940
## 405 405 -74.99886  9.817773  15.15561
## 406 406 -75.33084  9.477064  14.47138
## 407 407 -74.56643  9.031520        NA
## 408 408 -74.49435 10.009969  13.95596
## 409 409 -75.72397  9.402806  25.48588
## 410 410 -74.50964  9.197507        NA
## 411 411 -75.19980  8.478959        NA
## 412 412 -75.48591  9.509824  29.12171
## 413 413 -74.42665  9.378782        NA
## 414 414 -75.39854  8.537928  20.05832
## 415 415 -74.93771  8.719203  78.10825
## 416 416 -75.54924  9.791565  22.49571
## 417 417 -74.99450  9.955368  16.99869
## 418 418 -75.01634  8.778172  45.69968
## 419 419 -74.94645  9.404990  17.41001
## 420 420 -75.51211 10.038361  56.91507
## 421 421 -74.99886 10.044913  13.57961
## 422 422 -75.06220  9.592818  17.74710
## 423 423 -74.66471  8.706098  29.41102
## 424 424 -74.43975  9.204059        NA
## 425 425 -75.70212  8.304236  18.27722
## 426 426 -74.47470  8.959447        NA
## 427 427 -75.24566  9.487984  16.64798
## 428 428 -75.40728  8.968183  19.49796
## 429 429 -75.72615  9.752252        NA
## 430 430 -74.68000  8.695178  44.58877
## 431 431 -74.51182  9.649603  13.77651
## 432 432 -74.74988  8.675522  37.38385
## 433 433 -74.81759  9.439935        NA
## 434 434 -75.53614  9.680179  24.35584
## 435 435 -75.73489  9.826510        NA
## 436 436 -75.77201  9.551321        NA
## 437 437 -74.97921  8.286763  29.40895
## 438 438 -75.69339  8.990023  68.85669
## 439 439 -75.14083  8.592528  49.45717
## 440 440 -75.79385  8.374125  19.05748
## 441 441 -74.93771 10.007785  13.45835
## 442 442 -74.54677  9.415911  22.34438
## 443 443 -74.78265  8.655866  28.28709
## 444 444 -74.96829  9.280500  22.92261
## 445 445 -75.60384  9.568793        NA
## 446 446 -75.73925  8.385045  12.43064
## 447 447 -74.98358  9.608106  16.10545
## 448 448 -74.41354  8.585976        NA
## 449 449 -75.57763  9.503272  23.40678
## 450 450 -74.87001  8.944159        NA
## 451 451 -75.57982  8.304236        NA
## 452 452 -75.34176  9.865822  17.26930
## 453 453 -74.58390  9.236820        NA
## 454 454 -75.39418  9.396254  16.07709
## 455 455 -75.21945  8.430910  44.89405
## 456 456 -75.29808  9.385334  15.06358
## 457 457 -74.44849  9.902951  17.55591
## 458 458 -75.35049  8.806564  24.36405
## 459 459 -75.38326  9.898583  15.54631
## 460 460 -75.42475  8.813116  32.71609
## 461 461 -74.44849  9.302341        NA
## 462 462 -75.30463  8.579424  21.93456
## 463 463 -74.88966  9.739148  15.40869
## 464 464 -74.73023  9.953184  15.91918
## 465 465 -75.60821 10.099514        NA
## 466 466 -74.46159  9.147274        NA
## 467 467 -74.51838  8.555400        NA
## 468 468 -74.54677  9.446487  23.55996
## 469 469 -75.77638  9.151642        NA
## 470 470 -75.38326  8.570688  25.07263
## 471 471 -75.79385  8.898294        NA
## 472 472 -74.74988  8.863349  45.43716
## 473 473 -74.78483  8.551032  33.15522
## 474 474 -75.68902  9.688915        NA
## 475 475 -75.79822  9.837430        NA
## 476 476 -75.32647  8.834957  18.12231
## 477 477 -75.33739  8.968183  30.26302
## 478 478 -75.22819  8.581608  15.49244
## 479 479 -75.14083  9.090489  16.15890
## 480 480 -74.85472  8.834957        NA
## 481 481 -74.67563  9.145090        NA
## 482 482 -74.83943 10.040545  30.56700
## 483 483 -75.74799  9.605922        NA
## 484 484 -75.67373  9.529481        NA
## 485 485 -74.97921  8.629657  30.72937
## 486 486 -74.96829  8.527007  52.87319
## 487 487 -75.01852  9.677995  15.06964
## 488 488 -75.19106  9.140722  32.31061
## 489 489 -74.58827  9.629946  21.90367
## 490 490 -75.50993  8.509535  19.17318
## 491 491 -75.74581  9.817773        NA
## 492 492 -75.71086  8.616553  12.66419
## 493 493 -75.27842  9.088305  28.89117
## 494 494 -74.52056  9.846166  12.06861
## 495 495 -74.98794 10.101698  12.76425
## 496 496 -75.58200  9.383150  41.83559
## 497 497 -74.98139  9.350389  17.79132
## 498 498 -75.37015  9.304525  41.37276
## 499 499 -75.35923 10.123539  13.11658
## 500 500 -75.35923  9.804669  15.42140
sitios <- na.omit(sitios)
head(sitios)
##   id    longit     latit      soc
## 1  1 -75.23692 10.016521 16.77218
## 2  2 -75.17140 10.055834 15.50385
## 3  3 -74.96829  9.872374 14.66136
## 5  5 -74.76080  9.649603 30.02661
## 6  6 -74.62103  9.985944 14.23169
## 7  7 -74.59482  9.308893 30.80341
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m  # Print the map

5 Interpolación

Una vez que los datos de muestreo están limpios y preparados, se procede con las tareas de interpolación.

5.1 Creación del objeto gstat

Para interpolar, primero necesitamos crear un objeto de la clase gstat utilizando la función del mismo nombre. Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial, incluyendo la definición del modelo y los datos de calibración. La función gstat interpreta el tipo de modelo de interpolación deseado basándose en sus argumentos:

Sin modelo de variograma: IDW

Modelo de variograma, sin covariables: Kriging Ordinario

Se utilizan tres parámetros principales de la función gstat:

formula: La “fórmula” de predicción que especifica las variables dependientes e independientes (también conocidas como covariables).

data: Los datos de calibración (también conocidos como datos de entrenamiento).

model: El modelo de variograma (usado para Kriging).

5.2 IDW interpolación

Para interpolar el COS utilizando el método IDW, creamos el objeto gstat g1, especificando solo la formula y data.

g1 = gstat(formula = soc_igh_15_30_sucre ~ 1, data = nmuestras)

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función predict para realizar la interpolación y estimar los valores de COS. La función predict acepta un objeto raster (stars) como plantilla y un objeto gstat como modelo. El raster se utiliza para especificar las ubicaciones donde se realizarán las predicciones y, en Kriging Universal, para especificar los valores de las covariables.

# a simple copy
rrr = aggregate(geog.soc, 4)

Creamos un objeto raster (stars) con valores de celda iguales a 1. Este objeto se llama rrr.

rrr
## class       : SpatRaster 
## size        : 214, 159, 1  (nrow, ncol, nlyr)
## resolution  : 0.00873615, 0.00873615  (x, y)
## extent      : -75.8015, -74.41245, 8.276935, 10.14647  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30_sucre 
## min value   :            10.20426 
## max value   :           106.70802

Inicialmente, rrr es un objeto SpatRaster con las mismas dimensiones, resolución, extensión y CRS que los datos originales de soc_igh_15_30_sucre, mostrando un rango de valores de COS.

Definimos nuevos valores y nombres para rrr. Esto a menudo implica reasignar los valores de las celdas a un valor constante (como 1, como se menciona) y cambiar el nombre de la capa, lo que lo convierte en una cuadrícula de referencia para la interpolación.

values(rrr) <-1
names(rrr) <- "valor"
rrr
## class       : SpatRaster 
## size        : 214, 159, 1  (nrow, ncol, nlyr)
## resolution  : 0.00873615, 0.00873615  (x, y)
## extent      : -75.8015, -74.41245, 8.276935, 10.14647  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Después de la redefinición, rrr se convierte en un objeto SpatRaster donde todas las celdas tienen un valor de 1, y la capa se renombra a valor. Este rrr modificado actúa como una cuadrícula de predicción o plantilla para la interpolación.

Por ejemplo, la siguiente expresión interpola los valores de COS de acuerdo con el modelo definido en g1 y la plantilla de raster definida en stars.rrr:

stars.rrr = st_as_stars(rrr)
## [inverse distance weighted interpolation]
z1= predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

El objeto resultante de la interpolación es z1

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  11.45828 20.41773 24.93884 26.78872 32.35245 106.7079     0
## var1.var         NA       NA       NA      NaN       NA       NA 34026
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 159  -75.8  0.008736 +proj=longlat +datum=WGS8... [x]
## y    1 214  10.15 -0.008736 +proj=longlat +datum=WGS8... [y]
z1 = z1["var1.pred",,]
names(z1) = "soc"

Los dos atributos incluidos en el objeto z1 son:

var1.pred: Representa los valores predichos de COS por el método IDW. Muestra estadísticas como el mínimo, primer cuartil, mediana, media, tercer cuartil y máximo de los valores interpolados.

var1.var: Esta es la varianza de la predicción. Para la interpolación IDW, este atributo generalmente no se calcula (o es NA), ya que IDW es un método determinístico y no proporciona una medida de la incertidumbre de la predicción como lo hacen los métodos geoestadísticos.

Luego, se define una paleta de colores para la visualización del mapa de COS interpolado. La figura resultante muestra la superficie continua de COS interpolada utilizando el método IDW.

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m  # Print the map

5.3 OK interpolation

Los métodos de Kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos y asignar pesos en consecuencia al hacer predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función variogram. Esta función requiere dos argumentos: formula (especifica la variable dependiente y las covariables) y data (la capa de puntos con la variable dependiente y las covariables como atributos de punto). Por ejemplo, la siguiente expresión calcula el variograma empírico de muestras, sin covariables.

v_emp_ok = variogram(soc_igh_15_30_sucre ~ 1, data=nmuestras)

Luego, se grafica el variograma empírico para visualizar la autocorrelación espacial.

plot(v_emp_ok)

Existen varias formas de ajustar un modelo de variograma a un variograma empírico. Aquí se utiliza el ajuste automático mediante la función autofitVariogram del paquete automap. Esta función elige el mejor tipo de modelo y ajusta sus parámetros. Es importante tener en cuenta que autofitVariogram no funciona con objetos sf directamente, por lo que el objeto se convierte a SpatialPointsDataFrame (del paquete sp). El modelo ajustado puede ser visualizado con plot.

v_mod_ok = autofitVariogram(soc_igh_15_30_sucre ~ 1, as(nmuestras, "Spatial"))
plot(v_mod_ok)

El objeto resultante es una lista con varios componentes, incluyendo el variograma empírico y el modelo de variograma ajustado. El componente $var_model del objeto resultante contiene el modelo real:

v_mod_ok$var_model
##   model       psill    range kappa
## 1   Nug    45.11466     0.00   0.0
## 2   Ste 32809.58848 88417.12   0.4

Ahora, el modelo de variograma puede pasarse a la función gstat, y podemos continuar con la interpolación de Kriging Ordinario.

## [using ordinary kriging]
g2 = gstat(formula = soc_igh_15_30_sucre ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
z2 = z2["var1.pred",,]
names(z2) = "soc"
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "OK SOC interpolation [%]"
    )
m  # Print the map
colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.soc, opacity = 0.8, colorOptions = colores, group="RealWorld") %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Soil organic carbon [%]"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m  # Print the map

6 Evaluación de resultados

6.1 Evaluacion cualitativa

Se presenta una vista comparativa de las tres salidas de interpolación: los datos del “mundo real” (los datos de referencia), la interpolación IDW y la interpolación OK. Esta comparación visual permite una evaluación cualitativa de qué tan bien cada método reproduce la distribución espacial del COS en relación con los datos originales.

### run the following code from the console -- line-by-line
cv1 = gstat.cv(g1)
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
cv2 = gstat.cv(g2)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
cv1 = na.omit(cv1)
cv1
## class       : SpatialPointsDataFrame 
## features    : 383 
## extent      : -75.79385, -74.44849, 8.278027, 10.1432  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 12.9077961428334,       NA, 11.2577524185181, -26.3761122459709,     NA,    1 
## max values  : 96.4049971336662,       NA, 113.313278198242,  68.6534320168222,     NA,  383

Luego, se convierte el objeto cv1 a un formato adecuado para el análisis y visualización de los residuales. Se grafica la distribución de los residuales para evaluar visualmente el rendimiento del modelo.

cv1 = st_as_sf(cv1)
sp::bubble(as(cv1[, "residual"], "Spatial"))

Finalmente, se calculan los índices de precisión de la predicción, como el Error Cuadrático Medio (RMSE). Para IDW, el RMSE es:

# This is the RMSE value for the IDW interpolation
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 9.18764

Se repite el proceso de validación cruzada para los resultados de Kriging Ordinario y se convierte el objeto cv2. Se calcula el RMSE para los resultados de OK:

cv2 = st_as_sf(cv2)
# This is the RMSE value for the OK interpolation
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 8.953999

7 Resumen de resultados y precisión de la técnica de interpolación:

Según los resultados obtenidos, el Error Cuadrático Medio (RMSE) para la interpolación IDW es de 9.18764 mientras que para la interpolación de Kriging Ordinario (OK) es de 8.953999.

El RMSE es una métrica comúnmente utilizada para evaluar la precisión de los modelos de predicción, donde un valor de RMSE más bajo indica un mejor ajuste del modelo a los datos observados. En este caso, el RMSE del Kriging Ordinario es menor que el RMSE del IDW.

Por lo tanto, la técnica de interpolación de Kriging Ordinario parece ser más precisa para estimar el carbono orgánico del suelo en el departamento de Sucre, Colombia. Esto se debe a que Kriging, al ser un método geoestadístico probabilístico, no solo considera la distancia entre los puntos de muestreo como lo hace IDW, sino que también tiene en cuenta la autocorrelación espacial de los datos a través de un modelo de variograma. Esta capacidad de modelar la estructura espacial de la variabilidad de los datos le permite generar predicciones más robustas y, en este caso, más precisas.

8 Referencia

Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. Available at https://rpubs.com/ials2un/soc_interp.

sessionInfo()
## R version 4.5.1 (2025-06-13 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] curl_6.4.0     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4   leaflet_2.2.2 
##  [6] automap_1.1-20 gstat_2.1-4    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4     sass_0.4.10        class_7.3-23       KernSmooth_2.23-26
##  [5] lattice_0.22-7     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.4    
##  [9] grid_4.5.1         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.5.1        parallel_4.5.1    
## [33] spacetime_1.3-3    png_0.1-8          vctrs_0.6.5        R6_2.6.1          
## [37] zoo_1.8-14         proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.11.0      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.1.0         tidyselect_1.2.1  
## [49] tibble_3.3.0       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.5.1