TALLER/PARCIAL INTERPOLACIÓN
Daniela Guerrero Poveda
2025-02-25
- 1.Primero vamos a instalar aquellos paquetes o librerías que no hemos instalado
- 2. Ahora vamos a cargar las librerías que requerimos usar
- 3. Vamos a ajustar el directorio de trabajo
- 4. Declarar las variables
- 5. Leer y visualizar los archivos de datos (shp y tif) que se van a usar
- 6. Interpolar el valor de
pH usando IDW
- 6.1. Primero échele un vistazo a los valores de la variable pH
- 6.2. La siguiente línea define el modelo de interpolación
- 6.3. Luego de definir los parámetros de nuestro modelo de interpolación IDS
- 6.4. Visualice en R es resultado de la interpolación
- 6.5. Almacenar raster en archivo geotiff en su computador
- ¿Cuáles son las zonas más adecuadas en términos de ph del suelo para el desarrollo del jitomate?
- 7. Interpolar el valor
de pH usando Kriging
- 7.1. Crear el variograma empírico
- 7.2. Ajustar el variograma empírico a una función
- 7.3. Definir modelo de interpolación usando el variograma
- 7.4. Llevar a cabo la predicción
- 7.5. Ahora visalizamos la superficie de la prediccion
- 7.6. Ahora visalizamos la superficie de la varianza
- 7.7. RMSE para Krigging
- 7.7. RMSE para IDW
- Gracias por su atención! Fin del Script
1.Primero vamos a instalar aquellos paquetes o librerías que no hemos instalado
install.packages( c( "sf", "stars","leaflet", "gstat","automap","raster", "RColorBrewer"))2. Ahora vamos a cargar las librerías que requerimos usar
library(sf)
library(stars)
library(leaflet)
library(gstat)
library(automap)
library(raster)
library(RColorBrewer)3. Vamos a ajustar el directorio de trabajo
getwd()## [1] "C:/Users/dguer/Desktop/GEOMÁTICA/INTERPOLACIÒN/Interpolacion-20250220T160855Z-001/Interpolacion"4. Declarar las variables
ruta_aoi="./municipios_muestreo.shp"
ruta_puntos="./puntos_muestreo.shp"
ruta_raster="./dem_srtm_9377.tif"5. Leer y visualizar los archivos de datos (shp y tif) que se van a usar
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]] puntos = st_read(ruta_puntos)## Reading layer `puntos_muestreo' from data source
## `C:\Users\dguer\Desktop\GEOMÁTICA\INTERPOLACIÒN\Interpolacion-20250220T160855Z-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 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]]]]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)## 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)¿En qué sistema de coordenadas se encuentra el archivo de muestreo de suelos?
Este código nos ayudará a verificar el sistema coordenado, gracias a él podemos confirmar que su sistema coordenado es MAGNA-SIRGAS 2018 / Origen-Nacional
Adicionalmente, los atributos que nos muestran los archivos incluyen coordenadas de latitud, longitud, pH del suelo, y concentración de nutrientes como NH4, NO3, K2O, P2O5 y SOC
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]]6. Interpolar el valor de pH usando IDW
6.1. Primero échele un vistazo a los valores de la variable pH
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_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
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
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. 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] z = z["var1.pred",,] 6.4. Visualice en R es resultado de la interpolación
names(z) = "pH"
b = seq(4, 7.5, 0.1)
plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
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
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?
Según el gráfico anterior, podemos afirmar que las zonas más adecuadas para establecer un cultivo de jitomate son aquelllas que se visualizan en color azul, ya que su pH facilitará el rendimiento adecuado del cultivo.
7. Interpolar el valor de pH usando Kriging
7.1. Crear el variograma empírico
v_emp_ok = variogram(pH ~ 1, puntos_reproj)
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"))
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_reproj)7.4. Llevar a cabo la predicción
z2 = predict(g2, raster_dem)## [using ordinary kriging] print(z2)## 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 = z2["var1.pred",,]
names(prediccion) = "pH"
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 = z2 ["var1.var",,]
names(varianza) = "varianza pH"
b_var = seq(0.1, 1.4, 0.1)
plot(varianza, breaks = b_var, col = hcl.colors(length(b_var)-1, "Spectral"), reset = FALSE)
plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)
7.7. RMSE para Krigging
cv2 = gstat.cv(g2)
cv2 = st_as_sf(cv2) sp::bubble(as(cv2[, "residual"], "Spatial"))
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))## [1] 0.70660877.7. RMSE para IDW
cv3 = gstat.cv(g)
cv3 = st_as_sf(cv3) sp::bubble(as(cv3[, "residual"], "Spatial"))
sqrt(sum((cv3$var1.pred - cv3$observed)^2) / nrow(cv3))## [1] 0.7196336