interpolacion

2) Se leen los datos

Trabajaremos con los datos obtenidos en un cuaderno pasado de variables agronomicas presentadas por soilgrids

archivo <- ("phh2o_igh_15_30.tif")
(Ph <- rast(archivo))# se transforma a un raster

## class       : SpatRaster 
## size        : 1693, 1107, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8517000, -8240250, 779250, 1202500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : phh2o_igh_15_30.tif 
## name        : phh2o_igh_15_30

Pasamos los datos transformados por soils grid y los dejamos en unidades reales para el ph se presenta el siguiente factor

Phr<-Ph/10

Volvemos los datos en coordenadas reales

geog ="+proj=longlat +datum=WGS84"
(geog.ph = project(Phr, geog))

## class       : SpatRaster 
## size        : 1739, 1271, 1  (nrow, ncol, nlyr)
## resolution  : 0.002185969, 0.002185969  (x, y)
## extent      : -76.3331, -73.55473, 7.000841, 10.80224  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : phh2o_igh_15_30 
## min value   :             4.6 
## max value   :             7.3

Convertimos el objeto a un objeto stars

stars.ph = st_as_stars(geog.ph)

Graficamos para verificar si quedo bien la conversion

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.ph,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 3:8)) 
m

####3)Muestreo

con la confirmacion de que quedo correcta la transformacion vamos a tomar 500 muestras aleatoriamente de sitios reales

set.seed(123456)#da numeros aleatorios


(samples <- spatSample(geog.ph, 1000, "random", as.points=TRUE))#se obtienen 500 muestras aleatorias

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1000, 1  (geometries, attributes)
##  extent      : -76.332, -73.55582, 7.001934, 10.80115  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : phh2o_igh_15_30
##  type        :           <num>
##  values      :             5.4
##                          5.701
##                            5.3

Presenta 1000 puntos tomados al azar, con coordenadas wgs84 y una extension entre los -76 a los -73 y de los 7 a los 10 en sus respectivos ejes y valores que no van mas alla de 7 Ph

Ahora se pasara el objeto spat a un objeto mas simple dandonos la extension de los datos

(muestras <- sf::st_as_sf(samples))

nmuestras <- na.omit(muestras)# quitamos las datos NA presentes en la toma de los datos

Se configura para observar el muestreo hecho

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
ph <- muestras$phh2o_igh_15_30

id <- seq(1,1000,1)

(sitios <- data.frame(id, longit, latit, ph))

sitios <- na.omit(sitios)# quitamos los valores NA

head(sitios)#verificamos la organizacion del data frame

Y visualizamos los datos

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.ph,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 3:8)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$Phr, clusterOptions = markerClusterOptions())
m  # Print the map

Quitando los NA obtenemos un muestreo de 779 puntos aleatorios en el area selecionada, con esta informacion estamos listos para desarrollar nuestro ejercicio de interpolación

####4)Interpolación

Para realizar una interpolación es necesario crear un objeto gstat el cual se trabaja con un librerya del mismo nombre estos contiene la informacion necesaria para realizar la interpolacion por cualquiera de los dos metos que vamos a trabajar estos son la definicion de un modelo y la calibracion de los datos.

Dependiendo de los argumentos dados a gstat estaremos trabajando en IDW o en OK para esto definimos tres parametros la formula, el modelo y los datos.

#####5)Interpolación IDW

Definimos el modelo con el que vamos a trabajar

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

con el modelo realizado ya podemos hacer predicciones con la funcion predict esta acepta datos de un objeto raster-stars como un dem o un modelo-gstat como un g1 con esto podemos realizar de manera especifica la localizacion donde se hace la predicción y se puede saber los coeficientes de variación pero este ultimo unicamente en OK.

Vamos a crear un objeto llamado rrr el cual es la disminución de la resolución del objeto raster geog.ph

(rrr = aggregate(geog.ph, 4))

## class       : SpatRaster 
## size        : 435, 318, 1  (nrow, ncol, nlyr)
## resolution  : 0.008743876, 0.008743876  (x, y)
## extent      : -76.3331, -73.55254, 6.998655, 10.80224  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : phh2o_igh_15_30 
## min value   :        4.637860 
## max value   :        7.017218

Ahora lo transformamos en un objeto con unicamente valores de 1

values(rrr) <-1

names(rrr) <- "valor"

rrr

## class       : SpatRaster 
## size        : 435, 318, 1  (nrow, ncol, nlyr)
## resolution  : 0.008743876, 0.008743876  (x, y)
## extent      : -76.3331, -73.55254, 6.998655, 10.80224  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Convertimos el obejto de spat raster a star

stars.rrr = st_as_stars(rrr)

Y con el objeto en star podemos realizar las predicciones,generamos un obejeto llamado z1 el cual va dar numericamente el valor de los ph por celdas de rrr

z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean 3rd Qu.     Max.   NA's
## var1.pred  4.665118 5.240374 5.644352 5.590594 5.88517 6.647747      0
## var1.var         NA       NA       NA      NaN      NA       NA 138330
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 318 -76.33  0.008744 +proj=longlat +datum=WGS8... [x]
## y    1 435   10.8 -0.008744 +proj=longlat +datum=WGS8... [y]

El valor predicho 1 sera el Ph por lo tanto podemos nombrarlo como tal

z1 = z1["var1.pred",,]
names(z1) = "ph"

contruimos una paleta de color con dominio en el alcance de los datos

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 3:8, na.color = "transparent")

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 3:8)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$ph, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$ph,
    title = "IDW ph interpolation"
    )
m  # Print the map

###5.1) Interpolacion Ok La interpolación por el metodo de OK es necesario que halla un variograma en el modelo por lo tanto es importante calcularlo dado que al hacer uso de este se llega a cuantificar el patron que hay en la autocorrelación de los dato de manera que se le da peso a las predicciones halladas.

Sabiendo la importancia del variograma lo calcularemos de manera empirica con la funcion varigrama que depende del modelo y los datos

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

La graficamos para observar el comportamiento que tiene

plot(v_emp_ok)

Ahora ajustaremos este variograma al modelo usando la funcion autofit variogram, cabe aclarar que ademas de este hay distintos metodos para ajustar el variograma al modelo , ademas para usar esta funcion los objetos sf no funcionan

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

Lo graficamos y vemos como los puntos tinen tendencia a ajustarse a una curva

plot(v_mod_ok)

Al variograma ajustado le uniremos el variograma teorico

v_mod_ok$var_model

Con el variograma ya podemos realizar la interpolación

g2 = gstat(formula = phh2o_igh_15_30 ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)#toma tiempo la interpolación por tanto es recomendable hacerlo una vez y poner como comentario esta linea de codigo

## [using ordinary kriging]

Se cambia el nombre de la variable

z2 = z2["var1.pred",,]
names(z2) = "ph"

Se grafica para mostrar las predicciones en los puntos

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(z2,opacity = 0.7,colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),domain = 4:8)) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$ph, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$ph,
    title = "OK Ph interpolation" )
m

6) Evaluación de los datos

Se pueden observar las dos interpolaciones y los datos del mundo real para evaluar cual de las interpolaciones estuvo mas acertada

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 4:8, na.color = "transparent")

m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.ph, 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$ph,
    title = "Ph")
m

Se observa que krigin ordinario se ajusta visualmente mas cerca con los datos reales.

Pese a haber observado la grafica y dar una idea se necesita un metodo cuantitativo para definir que interpolación presenta mejor interpolación por lo tanto se usa el metodo de validación cruzada

En la Validación Cruzada Leave-One-Out se sigue una serie de pasos siendo estos: -Sacamos un punto de los datos de calibración -Hacemos una predicción para ese punto -Repetimos para todos los puntos Al final, lo que obtenemos es una tabla con un valor observado y un valor predicho para todos los puntos.

cv1 = na.omit(cv1)

cv1

## class       : SpatialPointsDataFrame 
## features    : 779 
## extent      : -76.30796, -73.61484, 7.001934, 10.78803  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,           residual, zscore, fold 
## min values  : 4.80710830317459,       NA, 4.65010833740234, -0.443239058503761,     NA,    1 
## max values  : 6.42764837612551,       NA, 6.65222978591919,  0.490665944373585,     NA,  779

Se pasa sf para poder graficarlo

cv1 = st_as_sf(cv1)

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

La distribucion de los residuos es muy cercana a cero Calculamos las raices del cuadrado medio del error

sqrt(sum((cv1$var1.pred - cv1$observed)^2)/nrow(cv1)) #cuadrado medio del error de IDW

## [1] 0.1498236

cv2 = st_as_sf(cv2)

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

## [1] 0.1332452

Ambos modelos presenta un gran ajuste a ser cercanos a cero por lo tanto fue un buen modelo predictivo con ambos metdos sin embargo mediante el metodo OK se presento menor RSME ajustandose a lo visto en los graficos el modelo OK presento un mejor ajuste a los datos del mundo real, siendo así para los datos que se obtuvieron es ideal usar el modelo de krigin ordinario.

####7) Reference Cite this work as: Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. Available at https://rpubs.com/ials2un/soc_interp.

interpolacion

Cristian David Eljadue Silva

2025-07-23

1) Introduccion

2) Se leen los datos

6) Evaluación de los datos