Parcial Interpolación R
Alejandra Serrano
2025-02-25
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)## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(stars)## Cargando paquete requerido: abindlibrary(leaflet)
library(gstat)
library(automap)
library(raster)## Cargando paquete requerido: splibrary(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/aleja/Downloads/PARCIAL R"setwd("C:/Users/aleja/OneDrive/Escritorio/4TO SEMESTRE/GEOMATICA/PARCIAL R")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):
## from to offset delta refsys point x/y
## x 1 1754 4966240 29.78 MAGNA-SIRGAS 2018 / Orige... FALSE [x]
## y 1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... FALSE [y] plot(raster_dem)
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\aleja\Downloads\PARCIAL R\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. 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)
}## 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]] 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
## 1 Santander Confines 6.331580 -73.25827 1533 4.8 0.08 10.3
## 2 Santander Confines 6.345463 -73.28759 1489 4.9 0.07 10.3
## 3 Santander Confines 6.357339 -73.28138 1482 4.8 0.08 18.2
## 4 Santander Confines 6.360368 -73.24903 1581 4.7 0.07 6.7
## 5 Santander Confines 6.340922 -73.20868 1878 4.8 0.15 24.0
## 6 Santander Confines 6.357226 -73.23005 1742 4.8 0.09 11.7
## 7 Santander Confines 6.339062 -73.23301 1752 4.9 0.09 9.4
## 8 Santander Confines 6.402768 -73.22524 1715 4.8 0.08 9.4
## 9 Santander Confines 6.365608 -73.21098 1878 4.8 0.10 19.3
## 10 Santander Paramo 6.384360 -73.19552 1912 4.9 0.12 19.1
## NO3.ppm K2O.ppm P2O5.ppm SOC.pct Sand.pct Silt.pct Clay.pct Slope.pct
## 1 8.4 173.5 10.1 2.46 40.7 32.6 26.6 10.5
## 2 5.4 122.9 8.9 1.62 51.2 26.1 22.6 9.2
## 3 0.9 119.3 6.2 2.95 51.9 35.4 12.7 7.4
## 4 2.3 181.9 12.1 1.02 37.3 36.1 26.5 12.6
## 5 12.3 130.1 7.3 2.48 65.6 30.3 4.1 8.5
## 6 1.2 131.3 4.8 1.60 40.8 46.5 12.6 7.5
## 7 2.0 73.5 11.5 3.65 55.0 42.5 2.5 5.6
## 8 5.9 194.0 6.6 2.12 33.7 45.4 20.9 13.7
## 9 0.3 273.5 3.4 3.81 64.4 24.9 10.6 1.5
## 10 0.1 128.9 18.8 3.14 39.1 50.7 10.3 15.8
## geometry
## 1 POINT (4971447 2257635)
## 2 POINT (4968206 2259171)
## 3 POINT (4968893 2260483)
## 4 POINT (4972470 2260816)
## 5 POINT (4976930 2258665)
## 6 POINT (4974568 2260467)
## 7 POINT (4974240 2258460)
## 8 POINT (4975101 2265500)
## 9 POINT (4976677 2261393)
## 10 POINT (4978386 2263464)¿Qué atributos tienen los puntos de muestreo de suelos? Department Municipali Latitude.D Longitude. Altitude.m pH EC.dSm NH4.ppm NO3.ppm K2O.ppm P2O5.ppm SOC.pct Sand.pct Silt.pct Clay.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)## 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]]Se reproyecta a ESPG-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 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.
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):
## from to offset delta refsys x/y
## x 1 1754 4966240 29.78 MAGNA-SIRGAS 2018 / Orige... [x]
## y 1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... [y]** 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
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.
plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
# Añadir al gráfico existente los puntos de muestreo
plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)
#Añadir al gráfico existente curvas de nivel basadas en el ráster con la predicción
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** interpolacion_idw="./ph_idw.tif"**Guardar el raster en la ruta almacenada en interpolacion_idw (formato geotiff)**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? Laz zonas aptas son las que tienen un pH entre 6 y 6,8.
#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 jitomate")
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
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):
## from to offset delta refsys x/y
## x 1 1754 4966240 29.78 MAGNA-SIRGAS 2018 / Orige... [x]
## y 1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... [y]# 7.5. Ahora visalizamos la superficie de la predicción:
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)
plot(prediccion, breaks = b_predict, col = hcl.colors(length(b_predict)-1, "Spectral"), reset = FALSE)
7.6. Ahora visalizamos 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)
plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)
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]cv2 = st_as_sf(cv2)
sp::bubble(as(cv2[, "residual"], "Spatial"))
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))## [1] 0.7066087cv3 = gstat.cv(g)## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [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.7196336Exporte los raster de la predicción y de varianza a archivos geotiff Escriba acá su código:
interpolacion_kriging="./ph_kriging.tif"
write_stars(prediccion, dsn = interpolacion_kriging)
interpolacion_k_varianza="./ph_kriging_var.tif"
write_stars(varianza, dsn = interpolacion_k_varianza) Conclusión del calculo de RMSE: Teniendo en cuenta los resultados anteriores donde el RMSE por Kriking es de 0,7066087 y el RMSE por IDW es de 0,7196336, podemos concluir que el método por Kriking es el mas apropiado para usar debido que su valor de RMSE es menor.