Examen interpolación

1.Primero se instalan aquellos paquetes o librerías que no hemos instalado previamente

install.packages( c( "sf", "stars","leaflet", "gstat","automap","raster", "RColorBrewer"))

Error in install.packages : Updating loaded packages

2. Ahora se cargan las librerías que requerimos usar las cuales ya debieron ser instaladas en el paso anterior

library(sf)
library(stars)
library(leaflet)
library(gstat)
library(automap)
library(raster)
library(RColorBrewer)

3.Se ajusta el directorio de trabajo con el fin que sea el mismo en el cual se encuentra almacenado el archivo R script

• El directorio de trabajo es aquel en el cual se guardan por defecto los archivos creados al trabajar con rutas relativas
• La siguiente línea nos muestra cual es el directorio de trabajo actual

getwd()

[1] "C:/Users/LENOVO/Desktop/INTREPOLACION_JOHA/interpolacion/Interpolacion-20250224T023201Z-001/Interpolacion"

    setwd("C:/Users/LENOVO/Desktop/INTREPOLACION_JOHA/interpolacion/Interpolacion-20250224T023201Z-001/Interpolacion")

4. Se declaran las variables a continuación

ruta_aoi="./municipios_muestreo.shp"
ruta_puntos="./puntos_muestreo.shp" #Ruta de archivo shp de puntos de muestreo
ruta_raster="./dem_srtm_9377.tif"   #Ruta de DEM (tif) de área de estudio

5. Leer y visualizar los archivos de datos (shp y tif) que se van a usar

• La siguiente línea lee el archivo del DEM como un elemento de R y lo almacena en la variable raster_dem, y contiene los siguientes componentes:
• read_stars: Es una función del paquete stars en R que se utiliza para leer archivos ráster, ruta_raster: Variable que contiene la ruta al archivo ráster que se leerá, RasterIO = Argumento que permite especificar opciones para la lectura del archivo.
• En este caso, RasterIO se usa para especificar que va a leer solamente la banda 1 del ráster (list(bands = 1)).
• El ráster del DEM solo tiene 1 banda pero es común trabajar con rásters multibanda, esta opción es muy útil para especificar la banda específica con la que se desea trabajar

 raster_dem=read_stars(ruta_raster, RasterIO = list(bands = 1))
    raster_dem

stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                   Min. 1st Qu. Median    Mean 3rd Qu. Max.
dem_srtm_9377.tif  -247       0      0 131.274       0 1673
dimension(s):

    plot(raster_dem)

downsample set to 1

    crs_raster=st_crs(raster_dem)
    crs_raster

Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

• La siguiente línea permite leer el archivo shapefile de los puntos de muestreo como un objeto de R

puntos = st_read(ruta_puntos)

Reading layer `puntos_muestreo' from data source 
  `C:\Users\LENOVO\Desktop\INTREPOLACION_JOHA\interpolacion\Interpolacion-20250224T023201Z-001\Interpolacion\puntos_muestreo.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 67 features and 16 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
Projected CRS: MAGNA-SIRGAS_Origen-Nacional

    gs_crs = st_crs(raster_dem)
    puntos = st_transform(puntos, crs = gs_crs)
    crs_puntos=st_crs(puntos)
    crs_puntos

Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

• Aunque aparentemente tienen el mismo sistema de referencia de coordenadas,existen diferencias en la representación del CRS, estas diferencias en la representación pueden hacer que R o los paquetes de manejo de CRS consideren los objetos como distintos, a pesar de que geométricamente podrían representar el mismo sistema de referencia.
• Se comparan los sistemas de referencia

if (identical(crs_raster, crs_puntos)) {
      cat("El sistema de referencia de coordenadas del raster y de los puntos es el mismo:\n")
      print(crs_raster)
    } else {
      cat("Los sistemas de referencia de coordenadas son diferentes.\n")
      cat("Raster:\n")
      print(crs_raster)
      cat("Puntos:\n")
      print(crs_puntos)
    }

El sistema de referencia de coordenadas del raster y de los puntos es el mismo:
Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

• A continuación se muestra la tabla de atributos de los puntos de muestreo

 puntos_reproj <-  st_transform(puntos, crs_raster)
    
    (puntos_reproj)

Simple feature collection with 67 features and 16 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
First 10 features:
   Department Municipali Latitude.D Longitude. Altitude.m  pH EC.dSm
1   Santander   Confines   6.331580  -73.25827       1533 4.8   0.08
2   Santander   Confines   6.345463  -73.28759       1489 4.9   0.07
3   Santander   Confines   6.357339  -73.28138       1482 4.8   0.08
4   Santander   Confines   6.360368  -73.24903       1581 4.7   0.07
5   Santander   Confines   6.340922  -73.20868       1878 4.8   0.15
6   Santander   Confines   6.357226  -73.23005       1742 4.8   0.09
7   Santander   Confines   6.339062  -73.23301       1752 4.9   0.09
8   Santander   Confines   6.402768  -73.22524       1715 4.8   0.08
9   Santander   Confines   6.365608  -73.21098       1878 4.8   0.10
10  Santander     Paramo   6.384360  -73.19552       1912 4.9   0.12
   NH4.ppm NO3.ppm K2O.ppm P2O5.ppm SOC.pct Sand.pct Silt.pct
1     10.3     8.4   173.5     10.1    2.46     40.7     32.6
2     10.3     5.4   122.9      8.9    1.62     51.2     26.1
3     18.2     0.9   119.3      6.2    2.95     51.9     35.4
4      6.7     2.3   181.9     12.1    1.02     37.3     36.1
5     24.0    12.3   130.1      7.3    2.48     65.6     30.3
6     11.7     1.2   131.3      4.8    1.60     40.8     46.5
7      9.4     2.0    73.5     11.5    3.65     55.0     42.5
8      9.4     5.9   194.0      6.6    2.12     33.7     45.4
9     19.3     0.3   273.5      3.4    3.81     64.4     24.9
10    19.1     0.1   128.9     18.8    3.14     39.1     50.7
   Clay.pct Slope.pct                geometry
1      26.6      10.5 POINT (4971447 2257635)
2      22.6       9.2 POINT (4968206 2259171)
3      12.7       7.4 POINT (4968893 2260483)
4      26.5      12.6 POINT (4972470 2260816)
5       4.1       8.5 POINT (4976930 2258665)
6      12.6       7.5 POINT (4974568 2260467)
7       2.5       5.6 POINT (4974240 2258460)
8      20.9      13.7 POINT (4975101 2265500)
9      10.6       1.5 POINT (4976677 2261393)
10     10.3      15.8 POINT (4978386 2263464)

Respondiendo las siguientes preguntas :

• ¿Qué atributos tienen los puntos de muestreo de suelos?

• ¿En qué sistema de coordenadas se encuentra el archivo de muestreo de suelos?escriba el código para revisar el Sistema de Coordenadas abajo

 st_crs(puntos)

Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

    st_crs(raster_dem)

Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

6. Interpolar el valor de pH usando IDW

6.1. Primero échele un vistazo a los valores de la variable pH

• var=“pH” : Define que la variable var almacenará la cadena de texto pH
• puntos[var] : Extrae solo la columna pH de la tabla de atributos de los puntos de muestreo

var="pH"
 puntos[var]

Simple feature collection with 67 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
First 10 features:
    pH                geometry
1  4.8 POINT (4971447 2257635)
2  4.9 POINT (4968206 2259171)
3  4.8 POINT (4968893 2260483)
4  4.7 POINT (4972470 2260816)
5  4.8 POINT (4976930 2258665)
6  4.8 POINT (4974568 2260467)
7  4.9 POINT (4974240 2258460)
8  4.8 POINT (4975101 2265500)
9  4.8 POINT (4976677 2261393)
10 4.9 POINT (4978386 2263464)

6.2. La siguiente línea define el modelo de interpolación

En la función gstat del paquete gstat, idw es el método de interpolación por defecto por este motivo no ve el término idw en ninguna parte.

• g : Es la variable en el cual se almacenará el resultado de la función gstat.

• gstat(): Es la función principal a la cual se le ingresan los siguientes argumentos: formula = pH ~ 1 : Fórmula en la cual la variable pH es la variable de respuesta (la que vamos a predecir), y ~ 1: indica que se está modelando pH como una constante (sin otras variables predictoras).

• El argumento data especifica la capa de puntos con valores de pH que se usará para generar la interpolación. Mediante el argumento set se pasan opciones adicionales a la función gstat.

• En este caso, idp=2 se refiere a la potencia a la que se eleva la distancia en el método idw.

st_crs(puntos_reproj)

Coordinate Reference System:
  User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
    BASEGEOGCRS["MAGNA-SIRGAS 2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",20046]],
    CONVERSION["Colombia Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-73,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9992,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, topographic mapping."],
        AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
        BBOX[-4.23,-84.77,15.51,-66.87]],
    ID["EPSG",9377]]

    crs(raster_dem)

[1] NA

g = gstat(formula = pH ~ 1, data = puntos_reproj,set=list(idp=2))

6.3. Luego de definir los parámetros de nuestro modelo de interpolación IDS creamos la superficie de predicción.

• predict : Es una función para realizar predicciones sobre un conjunto de datos espaciales
• g : Es la variable de tipo gstat que contiene la configuración del modelo idw
• raster_dem: Proporciona las coordenadas espaciales sobre las cuales se quiere hacer la interpolación, es decir, raster_dem define la “malla” o el raster plantilla sobre la que se calcularán las predicciones.
• z :Es la variable que almacenará los resultados de la interpolación, es decir que z contendrá un raster con los valores predichos para cada celda del ráster almacenado en raster_dem.
• print(z) : Muestra los resultados almacenados en z.

z = predict(g, raster_dem)

[inverse distance weighted interpolation]

    
    print(z)

stars object with 2 dimensions and 2 attributes
attribute(s), summary of first 1e+05 cells:
              Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
var1.pred  5.04681 5.130046 5.162891 5.174822 5.238671 5.283381
var1.var        NA       NA       NA      NaN       NA       NA
            NA's
var1.pred  0e+00
var1.var   1e+05
dimension(s):

• Notará que el raster que contiene la interpolación crea dos atributos:

  var1.pred : Contiene la suiperficie con los valores de las predicciones

  var1.var : Representa la varianza que se genera únicamente cuando se usa kriging

• Como para IDW no generamos un atributo de varianza,solo traemos el primer atributo de z y lo almacenamos en la variable z, es decir que reemplazamos el valor que tenía z con 2 atributos por el atributo var1.pred.

z = z["var1.pred",,]

6.4. Visualice en R es resultado de la interpolación

• Cambiamos el nombre del atributo de var1.pred a “pH” así el raster z tendrá un único atributo llamado “pH” que contiene los valores interpolados.

names(z) = "pH"

-La siguiente línea crea una secuencia de números desde 4 hasta 7.5 con un incremento de 0.1 y la almacena en la variable b.

b = seq(4, 7.5, 0.1)

• La siguiente línea genera un gráfico de la predicción.

  breaks = b: Define los intervalos de los datos usando la secuencia generada en el paso anterior. col = hcl.colors(length(b)-1, “Spectral”): Define una paleta de colores “Spectral” para visualizar los datosas.

• La cantidad de colores se determina como length(b)-1, es decir, un color para cada intervalo definido en b. reset = FALSE: Este argumento evita que los límites del gráfico se reseteen después de este plot, es útil porque se van a agregar más elementos al gráfico en los siguientes pasos.

• Se añaden al gráfico existente los puntos de muestreo.

  El argumento pch = 3 define el tipo de símbolo (en este caso, una pequeño cruz), add = TRUE indica que se deben añadir estos puntos al gráfico existente en lugar de crear un nuevo gráfico.

• Se añade al gráfico existente curvas de nivel basadas en el ráster con la predicción.

 plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)

downsample set to 1

  plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)

  contour(z, breaks = b, add = TRUE)

6.5. Almacenar raster en archivo geotiff en su computador

• Se define la variable con la ruta del archivo tif que contendrá z
• Se guarda el raster en la ruta almacenada en interpolacion_idw (formato geotiff)

interpolacion_idw="./ph_idw.tif"
write_stars(z, dsn = interpolacion_idw)    

Respondiendo la siguiente pregunta

• ¿Cuáles son las zonas más adecuadas en términos de ph del suelo para el desarrollo del jitomate?

 zona_optima = prediccion 
    zona_optima[zona_optima < 6.0 | zona_optima > 6.8] <- NA  
    plot(zona_optima, col = hcl.colors(10, "Blues"), main = "Zonas óptimas para el jitomate")

downsample set to 1

7. Interpolar el valor de pH usando Kriging

7.1. Crear el variograma empírico

• variogram es una función del paquete gstat que calcula el variograma empírico para los puntos de muestreo de suelos.

• pH ~ 1 correspone al modelo de fórmula que se está usando para calcular el variograma.

En este caso, se está interesado en el variograma del pH

• y ~ 1 Indica que no hay ninguna variable explicativa o covariable (solo pH)

• puntos Es una variable que almacena los puntos de muestreo con los valores de pH.

v_emp_ok = variogram(pH ~ 1, puntos_reproj)

• Se imprime el variograma empírico

plot(v_emp_ok)

7.2. Ajustar el variograma empírico a una función

v_mod_ok = autofitVariogram(pH ~ 1, as(puntos_reproj, "Spatial"))

• Se imprime el variograma ajustado

plot(v_mod_ok)

7.3. Definir modelo de interpolación usando el variograma

g2 = gstat(formula = pH ~ 1, model = v_mod_ok$var_model, data = puntos)

7.4. Llevar a cabo la predicción

 z = predict(g2, raster_dem)

[using ordinary kriging]

• Notará que el raster que contiene la interpolación crea dos atributos:

var1.pred que contiene la suiperficie con los valores de las predicciones

var1.var representa la varianza (desviación estándard al cuadrado)

  print(z)

stars object with 2 dimensions and 2 attributes
attribute(s), summary of first 1e+05 cells:
               Min.   1st Qu.    Median      Mean 3rd Qu.     Max.
var1.pred  4.973535 5.0923192 5.1475910 5.1348409 5.18763 5.238132
var1.var   0.604284 0.7349258 0.8883815 0.9006311 1.04760 1.302663
dimension(s):

7.5. Ahora visalizamos la superficie de la prediccion

prediccion = z["var1.pred",,]
  names(prediccion) = "pH"

• La siguiente línea crea una secuencia de números desde 4 hasta 7 con un incremento de 0.1 y la almacena en la variable b_predict.

b_predict = seq(4, 7, 0.1)  

• La siguiente línea muestra un gráfico de la predicción.

plot(prediccion, breaks = b_predict, col = hcl.colors(length(b_predict)-1, "Spectral"), reset = FALSE)

downsample set to 1

7.6. Ahora visualizamos la superficie de la varianza

varianza = z["var1.var",,]
names(varianza) = "varianza pH"

b_var = seq(0.1, 1.4, 0.1)

• La siguiente línea muestra la distribución de la varianza

 plot(varianza, breaks = b_var, col = hcl.colors(length(b_var)-1, "Spectral"), reset = FALSE)

downsample set to 1

  plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)

cv2 = gstat.cv(g2)

  |                                                                 
  |                                                           |   0%[using ordinary kriging]

  |                                                                 
  |=                                                          |   2%[using ordinary kriging]

  |                                                                 
  |==                                                         |   3%[using ordinary kriging]

  |                                                                 
  |===                                                        |   5%[using ordinary kriging]

  |                                                                 
  |====                                                       |   6%[using ordinary kriging]

  |                                                                 
  |====                                                       |   8%[using ordinary kriging]

  |                                                                 
  |=====                                                      |   9%[using ordinary kriging]

  |                                                                 
  |======                                                     |  11%[using ordinary kriging]

  |                                                                 
  |=======                                                    |  12%[using ordinary kriging]

  |                                                                 
  |========                                                   |  14%[using ordinary kriging]

  |                                                                 
  |=========                                                  |  15%[using ordinary kriging]

  |                                                                 
  |==========                                                 |  17%[using ordinary kriging]

  |                                                                 
  |===========                                                |  18%[using ordinary kriging]

  |                                                                 
  |============                                               |  20%[using ordinary kriging]

  |                                                                 
  |=============                                              |  21%[using ordinary kriging]

  |                                                                 
  |=============                                              |  23%[using ordinary kriging]

  |                                                                 
  |==============                                             |  24%[using ordinary kriging]

  |                                                                 
  |===============                                            |  26%[using ordinary kriging]

  |                                                                 
  |================                                           |  27%[using ordinary kriging]

  |                                                                 
  |=================                                          |  29%[using ordinary kriging]

  |                                                                 
  |==================                                         |  30%[using ordinary kriging]

  |                                                                 
  |===================                                        |  32%[using ordinary kriging]

  |                                                                 
  |====================                                       |  33%[using ordinary kriging]

  |                                                                 
  |=====================                                      |  35%[using ordinary kriging]

  |                                                                 
  |=====================                                      |  36%[using ordinary kriging]

  |                                                                 
  |======================                                     |  38%[using ordinary kriging]

  |                                                                 
  |=======================                                    |  39%[using ordinary kriging]

  |                                                                 
  |========================                                   |  41%[using ordinary kriging]

  |                                                                 
  |=========================                                  |  42%[using ordinary kriging]

  |                                                                 
  |==========================                                 |  44%[using ordinary kriging]

  |                                                                 
  |===========================                                |  45%[using ordinary kriging]

  |                                                                 
  |============================                               |  47%[using ordinary kriging]

  |                                                                 
  |=============================                              |  48%[using ordinary kriging]

  |                                                                 
  |==============================                             |  50%[using ordinary kriging]

  |                                                                 
  |==============================                             |  52%[using ordinary kriging]

  |                                                                 
  |===============================                            |  53%[using ordinary kriging]

  |                                                                 
  |================================                           |  55%[using ordinary kriging]

  |                                                                 
  |=================================                          |  56%[using ordinary kriging]

  |                                                                 
  |==================================                         |  58%[using ordinary kriging]

  |                                                                 
  |===================================                        |  59%[using ordinary kriging]

  |                                                                 
  |====================================                       |  61%[using ordinary kriging]

  |                                                                 
  |=====================================                      |  62%[using ordinary kriging]

  |                                                                 
  |======================================                     |  64%[using ordinary kriging]

  |                                                                 
  |======================================                     |  65%[using ordinary kriging]

  |                                                                 
  |=======================================                    |  67%[using ordinary kriging]

  |                                                                 
  |========================================                   |  68%[using ordinary kriging]

  |                                                                 
  |=========================================                  |  70%[using ordinary kriging]

  |                                                                 
  |==========================================                 |  71%[using ordinary kriging]

  |                                                                 
  |===========================================                |  73%[using ordinary kriging]

  |                                                                 
  |============================================               |  74%[using ordinary kriging]

  |                                                                 
  |=============================================              |  76%[using ordinary kriging]

  |                                                                 
  |==============================================             |  77%[using ordinary kriging]

  |                                                                 
  |==============================================             |  79%[using ordinary kriging]

  |                                                                 
  |===============================================            |  80%[using ordinary kriging]

  |                                                                 
  |================================================           |  82%[using ordinary kriging]

  |                                                                 
  |=================================================          |  83%[using ordinary kriging]

  |                                                                 
  |==================================================         |  85%[using ordinary kriging]

  |                                                                 
  |===================================================        |  86%[using ordinary kriging]

  |                                                                 
  |====================================================       |  88%[using ordinary kriging]

  |                                                                 
  |=====================================================      |  89%[using ordinary kriging]

  |                                                                 
  |======================================================     |  91%[using ordinary kriging]

  |                                                                 
  |=======================================================    |  92%[using ordinary kriging]

  |                                                                 
  |=======================================================    |  94%[using ordinary kriging]

  |                                                                 
  |========================================================   |  95%[using ordinary kriging]

  |                                                                 
  |=========================================================  |  97%[using ordinary kriging]

  |                                                                 
  |========================================================== |  98%[using ordinary kriging]

  |                                                                 
  |===========================================================| 100%[using ordinary kriging]

cv2 = st_as_sf(cv2)

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

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

[1] 0.7066087

 cv3 = gstat.cv(g)

  |                                                                 
  |                                                           |   0%[inverse distance weighted interpolation]

  |                                                                 
  |=                                                          |   2%[inverse distance weighted interpolation]

  |                                                                 
  |==                                                         |   3%[inverse distance weighted interpolation]

  |                                                                 
  |===                                                        |   5%[inverse distance weighted interpolation]

  |                                                                 
  |====                                                       |   6%[inverse distance weighted interpolation]

  |                                                                 
  |====                                                       |   8%[inverse distance weighted interpolation]

  |                                                                 
  |=====                                                      |   9%[inverse distance weighted interpolation]

  |                                                                 
  |======                                                     |  11%[inverse distance weighted interpolation]

  |                                                                 
  |=======                                                    |  12%[inverse distance weighted interpolation]

  |                                                                 
  |========                                                   |  14%[inverse distance weighted interpolation]

  |                                                                 
  |=========                                                  |  15%[inverse distance weighted interpolation]

  |                                                                 
  |==========                                                 |  17%[inverse distance weighted interpolation]

  |                                                                 
  |===========                                                |  18%[inverse distance weighted interpolation]

  |                                                                 
  |============                                               |  20%[inverse distance weighted interpolation]

  |                                                                 
  |=============                                              |  21%[inverse distance weighted interpolation]

  |                                                                 
  |=============                                              |  23%[inverse distance weighted interpolation]

  |                                                                 
  |==============                                             |  24%[inverse distance weighted interpolation]

  |                                                                 
  |===============                                            |  26%[inverse distance weighted interpolation]

  |                                                                 
  |================                                           |  27%[inverse distance weighted interpolation]

  |                                                                 
  |=================                                          |  29%[inverse distance weighted interpolation]

  |                                                                 
  |==================                                         |  30%[inverse distance weighted interpolation]

  |                                                                 
  |===================                                        |  32%[inverse distance weighted interpolation]

  |                                                                 
  |====================                                       |  33%[inverse distance weighted interpolation]

  |                                                                 
  |=====================                                      |  35%[inverse distance weighted interpolation]

  |                                                                 
  |=====================                                      |  36%[inverse distance weighted interpolation]

  |                                                                 
  |======================                                     |  38%[inverse distance weighted interpolation]

  |                                                                 
  |=======================                                    |  39%[inverse distance weighted interpolation]

  |                                                                 
  |========================                                   |  41%[inverse distance weighted interpolation]

  |                                                                 
  |=========================                                  |  42%[inverse distance weighted interpolation]

  |                                                                 
  |==========================                                 |  44%[inverse distance weighted interpolation]

  |                                                                 
  |===========================                                |  45%[inverse distance weighted interpolation]

  |                                                                 
  |============================                               |  47%[inverse distance weighted interpolation]

  |                                                                 
  |=============================                              |  48%[inverse distance weighted interpolation]

  |                                                                 
  |==============================                             |  50%[inverse distance weighted interpolation]

  |                                                                 
  |==============================                             |  52%[inverse distance weighted interpolation]

  |                                                                 
  |===============================                            |  53%[inverse distance weighted interpolation]

  |                                                                 
  |================================                           |  55%[inverse distance weighted interpolation]

  |                                                                 
  |=================================                          |  56%[inverse distance weighted interpolation]

  |                                                                 
  |==================================                         |  58%[inverse distance weighted interpolation]

  |                                                                 
  |===================================                        |  59%[inverse distance weighted interpolation]

  |                                                                 
  |====================================                       |  61%[inverse distance weighted interpolation]

  |                                                                 
  |=====================================                      |  62%[inverse distance weighted interpolation]

  |                                                                 
  |======================================                     |  64%[inverse distance weighted interpolation]

  |                                                                 
  |======================================                     |  65%[inverse distance weighted interpolation]

  |                                                                 
  |=======================================                    |  67%[inverse distance weighted interpolation]

  |                                                                 
  |========================================                   |  68%[inverse distance weighted interpolation]

  |                                                                 
  |=========================================                  |  70%[inverse distance weighted interpolation]

  |                                                                 
  |==========================================                 |  71%[inverse distance weighted interpolation]

  |                                                                 
  |===========================================                |  73%[inverse distance weighted interpolation]

  |                                                                 
  |============================================               |  74%[inverse distance weighted interpolation]

  |                                                                 
  |=============================================              |  76%[inverse distance weighted interpolation]

  |                                                                 
  |==============================================             |  77%[inverse distance weighted interpolation]

  |                                                                 
  |==============================================             |  79%[inverse distance weighted interpolation]

  |                                                                 
  |===============================================            |  80%[inverse distance weighted interpolation]

  |                                                                 
  |================================================           |  82%[inverse distance weighted interpolation]

  |                                                                 
  |=================================================          |  83%[inverse distance weighted interpolation]

  |                                                                 
  |==================================================         |  85%[inverse distance weighted interpolation]

  |                                                                 
  |===================================================        |  86%[inverse distance weighted interpolation]

  |                                                                 
  |====================================================       |  88%[inverse distance weighted interpolation]

  |                                                                 
  |=====================================================      |  89%[inverse distance weighted interpolation]

  |                                                                 
  |======================================================     |  91%[inverse distance weighted interpolation]

  |                                                                 
  |=======================================================    |  92%[inverse distance weighted interpolation]

  |                                                                 
  |=======================================================    |  94%[inverse distance weighted interpolation]

  |                                                                 
  |========================================================   |  95%[inverse distance weighted interpolation]

  |                                                                 
  |=========================================================  |  97%[inverse distance weighted interpolation]

  |                                                                 
  |========================================================== |  98%[inverse distance weighted interpolation]

  |                                                                 
  |===========================================================| 100%[inverse distance weighted interpolation]

  cv3 = st_as_sf(cv3)

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

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

[1] 0.7196336

• Exportamos los raster de la predicción y de varianza a archivos geotiff
• Se escribe el código

write_stars(prediccion, dsn = "./ph_kriging.tif")
  write_stars(varianza, dsn = "./ph_varianza.tif")

Cálculo del RMSE

• El cálculo del RMSE en Kriging es de 0.7066087 y el RMSE de IDW es de 0.7196336 por lo tanto podríamos decir que el mejor es el de Kriging al tener un menor valor.