Interpolacion Espacial
Stefany Torres
2025-02-23
options(repos = c(CRAN = "https://cloud.r-project.org/"))En este script utilizaremos dos tipos de datos
Datos de muestreos de propiedades químicas y físicas de los suelos en sistemas productivos de tomate en el departamento de Santander en formato shapefile. Los datos completos en formato csv se encuentran disponibles en https://www.sciencedirect.com/science/article/pii/S2352340919301957
Archivo ráster con Modelo Digital de Terreno SRTM. Descargado de Google Earth Engine
Realizaremos las siguientes actividades:
- Instalar las librerías/paquetes con contienen las funciones que vamos a utilizar
- Traer las librerías/paquetes con contienen las funciones que vamos a utilizar
- Ajustar el directorio de trabajo
- Declarar las variables que vamos a utilizar
- Leer y visualizar los archivos de datos (shp y tif) que se van a usar
- Interpolar valores de pH del suelo usando el método IDW
- Interpolar valores de pH del suelo usando el método Kriging
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
library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(stars)## Warning: package 'stars' was built under R version 4.3.3## Loading required package: abind## Warning: package 'abind' was built under R version 4.3.3library(leaflet)## Warning: package 'leaflet' was built under R version 4.3.3library(gstat)## Warning: package 'gstat' was built under R version 4.3.3library(automap)## Warning: package 'automap' was built under R version 4.3.3library(raster)## Warning: package 'raster' was built under R version 4.3.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.3.3library(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/EDGAR/Desktop/geomatica/interpolacion"#Si no corresponde al directorio dentro del cual se encuentra almacenado este archivo
#vaya a Session-Set working directory- to source file location
#Revise nuevamente cual es su directorio de trabajo4. 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 estudio5. 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 alarchivo 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\EDGAR\Desktop\geomatica\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-Nacionalcrs_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]]]]reproyectar
puntos_reproj <- st_transform(puntos, crs_raster) #Reproyectar vector puntos a 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)crs_puntos_reproj= st_crs(puntos_reproj)
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]]¿Qué atributos tienen los puntos de muestreo de suelos?l
- Departamento,
- Municipio,
- Latitud,
- Altitud en metros,
- pH,
- Conductividad eléctrica del suelo, medida en decisiemens por metro (dS/m),
- Concentración de amonio (NH₄⁺) en partes por millón (ppm
- Concentración de nitrato (NO₃⁻) en partes por millón (ppm)
- Concentración de óxido de potasio (K₂O) en partes por millón
- Concentración de pentóxido de fósforo (P₂O₅) en partes por millón (ppm)
- Contenido de carbono orgánico del suelo (SOC) en porcentaje,
- Porcentaje de arena en el suelo,
- Porcentaje de limo en el suelo,
- Porcentaje de arcilla en el suelo,
- Porcentaje de pendiente del terreno,
- Geometría del punto de muestreo
¿En qué sistema de coordenadas se encuentra el archivo de muestreo de suelos? MAGNA-SIRGAS_Origen-Nacional
#código para revisar el Sistema de Coordenadas
# archivo shapefile
soil_data = st_read(ruta_aoi)## Reading layer `municipios_muestreo' from data source
## `C:\Users\EDGAR\Desktop\geomatica\interpolacion\municipios_muestreo.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 90 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 4966256 ymin: 2254942 xmax: 5018458 ymax: 2297432
## Projected CRS: MAGNA-SIRGAS_Origen-Nacional# Verificar el sistema de coordenadas
crs = st_crs(soil_data)
print(crs)## 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]]]]6. Interpolar el valor de pH usando IDW (Inverse Distance Weighting)
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
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=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]]g_idw = gstat(formula = pH ~ 1, data = puntos_reproj,set=list(idp=2)) - 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.
6.3 almacenar los valores de la prediccion Notará que el raster que contiene la interpolación crea dos atributos: var1.pred que contiene la suiperficie con los valores de las predicciones y var1.var representa la varianza que se genera únicamente cuando se usa kriging
z = predict(g_idw, 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]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",,] #guardar solo los datos de la prediccion6.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) 6.5 gráfico de la predicción
plot_iwf = plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
#Añadir al gráfico existente los puntos de muestreo- 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
# Crear el gráfico base
plot_iwf <- plot(z, breaks = b, col = hcl.colors(length(b) - 1, "Spectral"), reset = FALSE)
# Añadir los puntos de muestreo al gráfico existente
points(st_geometry(puntos_reproj), pch = 3)
# Añadir curvas de nivel al gráfico existente
contour(z, breaks = b, add = TRUE)
- 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.
6.6. Almacenar raster en archivo geotiff en su computador Definir variable con la ruta del archivo tif que contendrá z
interpolacion_idw="./ph_idwStefa.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?
Las zonas más adecuadas para el cultivo de jitomate son aquellas representadas de color azul claro y agua marian, ya que estos colores indican valores de pH más cercanos al rango óptimo de 6.0-6.5.
https://www.yara.com.mx/nutricion-vegetal/tomate/principios-agronomicos-en-tomate/
7. Interpolar el valor de pH usando Kriging
7.1. Crear el variograma empírico
v_emp_ok = variogram(pH ~ 1, puntos_reproj)
# Imprimir el variograma empírico
plot(v_emp_ok)
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_reproj: variable que almacena los puntos de muestreo con los valores de pH.
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):
## 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 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)
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)
7.7 Exporte los raster de la predicción y de varianza a archivos geotiff
# Definir variable con la ruta del archivo tif que contendrá z
interpolacion_krig="./ph_kri_predic.tif"
# Guardar el raster en la ruta almacenada en interpolacion_krig (formato geotiff)
write_stars(prediccion, dsn = interpolacion_krig) #guarda el texto de la ruta
# Definir variable con la ruta del archivo tif que contendrá z
interpolacion_krig_var="./ph_kri_var.tif"
# Guardar el raster en la ruta almacenada en interpolacion_krig_var (formato geotiff)
write_stars(varianza, dsn = interpolacion_krig_var) #guarda el texto de la ruta8. Validación Cruzada de Leave-One-Out
8.1 Validación Cruzada de idw Validación Cruzada con gstat.cv
cv_idw = gstat.cv(g_idw)## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]gstat. cv : Aplica validación cruzada al modelo de interpolación y Evalúa el desempeño del modelo dejando un punto fuera en cada iteración y prediciendo su valor con el resto de los datos.
cv_idw Guarda los resultados de la validación cruzada.
Conversión del Objeto cv_idw
cv_idw_sf = st_as_sf(cv_idw)st_as_sf : Convierte cv_idw en un objeto sf para poder trabajar con él en análisis y visualización espacial.
sf : es una librería que maneja datos espaciales de forma eficiente en R.
Graficar los Residuos
sp::bubble(as(cv_idw_sf[, "residual"], "Spatial"))
cv_idw_sf[, “residual”]: Extrae solo la columna residual del objeto cv_idw_sf.
Convierte el objeto sf (cv_idw_sf) a un objeto de tipo SpatialPointsDataFrame del paquete {sp}.
Cálculo del rmse error del cuadratico medio
rmse_idw = sqrt(sum((cv_idw_sf$var1.pred - cv_idw_sf$observed)^2) / nrow(cv_idw_sf))cv_idw_sf$var.pred : Contiene los valores predichos por la interpolación IDW.
cv_idw_sf$observed : Contiene los valores observados (datos reales).
(cv_idw_sf\(var1.pred - cv_idw_sf\)observed)^2): Calcula los errores cuadráticos para cada punto.
sum: Suma todos los errores cuadráticos.
/ nrow(cv_idw_sf): Divide la suma por el número total de observaciones (promedio del error cuadrático).
sqrt: Aplica la raíz cuadrada para obtener el RMSE.
8.2 Validación Cruzada de Kriging
Validación Cruzada con gstat.cv
cv_OK = gstat.cv(g)## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]Conversión del Objeto cv_OK
cv_OK_sf = st_as_sf(cv_OK)
head(cv_OK_sf)## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 4968206 ymin: 2257635 xmax: 4976930 ymax: 2260816
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## var1.pred var1.var observed residual zscore fold
## 1 4.829246 0.5118545 4.8 -0.02924598 -0.04087831 1
## 2 4.821140 0.5213360 4.9 0.07886035 0.10921942 2
## 3 4.862078 0.4961376 4.8 -0.06207771 -0.08813220 3
## 4 4.845078 0.4650866 4.7 -0.14507779 -0.21273260 4
## 5 4.841381 0.4939118 4.8 -0.04138118 -0.05888141 5
## 6 4.825614 0.4370995 4.8 -0.02561391 -0.03874233 6
## 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)Graficar los Residuos
sp::bubble(as(cv_OK_sf[, "residual"], "Spatial"))
Cálculo del RMSE
rmse_OK = sqrt(sum((cv_OK_sf$var1.pred - cv_OK_sf$observed)^2) / nrow(cv_OK_sf))7.3 Comparar los Resultados
print(paste("RMSE para OK:", rmse_OK))## [1] "RMSE para OK: 0.706608746956368"print(paste("RMSE para IDW:", rmse_idw))## [1] "RMSE para IDW: 0.719633591538204"if (rmse_idw < rmse_OK) {
print("IDW es más preciso que OK.")
} else {
print("OK es más preciso que IDW.")
}## [1] "OK es más preciso que IDW."# "OK es más preciso que IDW."