Taller examen interpolación - Karen González

1. Primero vamos a instalar aquellos paquetes o librerías que no hemos instalado previamente.

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

2. Ahora vamos a cargar 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. Vamos a ajustar 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/Karen González/OneDrive/Escritorio/Taller examen geomática/Interpolacion"

4. Declarar las variables

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 Tiene los siguientes componentes:

read_stars: Es una función del paquete stars en R que se utiliza para leer archivos ráster. ruta_raster: CVariable 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\Karen González\OneDrive\Escritorio\Taller examen geomática\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

    crs_puntos=st_crs(puntos)
    crs_puntos

Coordinate Reference System:
  User input: MAGNA-SIRGAS_Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS_Origen-Nacional",
    BASEGEOGCRS["GCS_MAGNA-SIRGAS_2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ID["EPSG",1329]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

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.

Comparar 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)
    }

Los sistemas de referencia de coordenadas son diferentes.
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]]
Puntos:
Coordinate Reference System:
  User input: MAGNA-SIRGAS_Origen-Nacional 
  wkt:
PROJCRS["MAGNA-SIRGAS_Origen-Nacional",
    BASEGEOGCRS["GCS_MAGNA-SIRGAS_2018",
        DATUM["Marco Geocentrico Nacional de Referencia 2018",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ID["EPSG",1329]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

puntos_reproj <-  st_transform(puntos, crs_raster)
    
    (puntos_reproj) #Muestra la tabla de atributos de los puntos de muestreo

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 NH4.ppm NO3.ppm K2O.ppm P2O5.ppm
1   Santander   Confines   6.331580  -73.25827       1533 4.8   0.08    10.3     8.4   173.5     10.1
2   Santander   Confines   6.345463  -73.28759       1489 4.9   0.07    10.3     5.4   122.9      8.9
3   Santander   Confines   6.357339  -73.28138       1482 4.8   0.08    18.2     0.9   119.3      6.2
4   Santander   Confines   6.360368  -73.24903       1581 4.7   0.07     6.7     2.3   181.9     12.1
5   Santander   Confines   6.340922  -73.20868       1878 4.8   0.15    24.0    12.3   130.1      7.3
6   Santander   Confines   6.357226  -73.23005       1742 4.8   0.09    11.7     1.2   131.3      4.8
7   Santander   Confines   6.339062  -73.23301       1752 4.9   0.09     9.4     2.0    73.5     11.5
8   Santander   Confines   6.402768  -73.22524       1715 4.8   0.08     9.4     5.9   194.0      6.6
9   Santander   Confines   6.365608  -73.21098       1878 4.8   0.10    19.3     0.3   273.5      3.4
10  Santander     Paramo   6.384360  -73.19552       1912 4.9   0.12    19.1     0.1   128.9     18.8
   SOC.pct Sand.pct Silt.pct Clay.pct Slope.pct                geometry
1     2.46     40.7     32.6     26.6      10.5 POINT (4971447 2257635)
2     1.62     51.2     26.1     22.6       9.2 POINT (4968206 2259171)
3     2.95     51.9     35.4     12.7       7.4 POINT (4968893 2260483)
4     1.02     37.3     36.1     26.5      12.6 POINT (4972470 2260816)
5     2.48     65.6     30.3      4.1       8.5 POINT (4976930 2258665)
6     1.60     40.8     46.5     12.6       7.5 POINT (4974568 2260467)
7     3.65     55.0     42.5      2.5       5.6 POINT (4974240 2258460)
8     2.12     33.7     45.4     20.9      13.7 POINT (4975101 2265500)
9     3.81     64.4     24.9     10.6       1.5 POINT (4976677 2261393)
10    3.14     39.1     50.7     10.3      15.8 POINT (4978386 2263464)

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

Atributos: Department, Municipali, Latitude.D, Longitude., Altitude.m, pH, EC.dSm, NH4.ppm, NO3.ppm, K2O.ppm, P2O5.ppm, SOC.pct.

¿En qué sistema de coordenadas se encuentra el archivo de muestreo de suelos? MAGNA-SIRGAS 2018 / Origen-Nacional

Escriba el código para revisar el Sistema de Coordenadas abajo:

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]]

Se reproyectó a “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

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_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.

z = predict(g, raster_dem)

[inverse distance weighted interpolation]

print(z) #  Muestra los resultados almacenados en 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.  NA's
var1.pred  5.04681 5.130046 5.162891 5.174822 5.238671 5.283381 0e+00
var1.var        NA       NA       NA      NaN       NA       NA 1e+05
dimension(s):

    #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 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 datos.

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.

Añadir 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.

Añadir 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.

Definir variable con la ruta del archivo tif que contendrá z

Guardar el raster en la ruta almacenada en interpolacion_idw (formato geotiff)

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

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

Hay dos zonas aptas para el desarrollo del culivo de tomate, que son las del pH 6 y 6.5 (color azúl oscuro), esta hortaliza difícilmente tolera pH más ácidos que estos.

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: variable que almacena los puntos de muestreo con los valores de pH.

v_emp_ok = variogram(pH ~ 1, puntos_reproj)
  # Imprimir 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"))
#Imprimir el variograma ajustado
    plot(v_mod_ok)

7.3. Definir modelo de interpolación usando el variograma

g = gstat(formula = pH ~ 1, model = v_mod_ok$var_model, data = puntos_reproj)
  

7.4. Llevar a cabo la predicción

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)

 z = predict(g, raster_dem)

[using ordinary kriging]

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 visualizamos 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

NA

Exporte los raster de la predicción y de varianza a archivos geotiff

Escriba acá su código:

Exportar geotiff para la predicción

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

Exportar geotiff para la varianza

interpolacion_k_varianza="./ph_kriging_var.tif"
write_stars(varianza, dsn = interpolacion_k_varianza)