rm(list=ls())
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(stars)
## Cargando paquete requerido: abind
library(leaflet)
library(gstat)
library(automap)
library(raster)
## Cargando paquete requerido: sp
library(RColorBrewer)
getwd()
## [1] "C:/Users/Home/Downloads/TALLER_PARCIAL2_GEOMATICA/drive-download-20250225T221640Z-001"
setwd("C:/Users/Home/Downloads/TALLER_PARCIAL2_GEOMATICA/drive-download-20250225T221640Z-001")
ruta_aoi="./municipios_muestreo.shp"
ruta_puntos="./puntos_muestreo.shp"
ruta_raster="./dem_srtm_9377.tif"
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\Home\Downloads\TALLER_PARCIAL2_GEOMATICA\drive-download-20250225T221640Z-001\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]]
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)
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]]
var="pH"
puntos[var]
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))
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",,]
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)
interpolacion_idw="./ph_idw.tif"
write_stars(z, dsn = interpolacion_idw)
v_emp_ok = variogram(pH ~ 1, puntos_reproj)
plot(v_emp_ok)
v_mod_ok = autofitVariogram(pH ~ 1, as(puntos_reproj, "Spatial"))
plot(v_mod_ok)
g = gstat(formula = pH ~ 1, model = v_mod_ok$var_model, data = puntos)
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]
prediccion = z["var1.pred",,]
names(prediccion) = "pH"
zona_optima = prediccion
zona_optima[zona_optima < 6.0 | zona_optima > 6.8] <- NA
plot(zona_optima, col = hcl.colors(10, "Blues"), main = "Zonas óptimas para el jitomate")
b_predict = seq(4, 7, 0.1)
plot(prediccion, breaks = b_predict, col = hcl.colors(length(b_predict)-1, "Spectral"), reset = FALSE)
varianza = z["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)
cv2 = 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]
cv2 = st_as_sf(cv2)
sp::bubble(as(cv2[, "residual"], "Spatial"))
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 0.7066087
cv3 = gstat.cv(g)
## [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]
cv3 = st_as_sf(cv3)
sp::bubble(as(cv3[, "residual"], "Spatial"))
sqrt(sum((cv3$var1.pred - cv3$observed)^2) / nrow(cv3))
## [1] 0.7066087
write_stars(prediccion, dsn = "./ph_kriging.tif")
write_stars(varianza, dsn = "./ph_varianza.tif")
¿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”
¿En qué sistema de coordenadas se encuentra el archivo de muestreo
de suelos?
MAGNA-SIRGAS 2018 / Origen-Nacional
¿Cuáles son las zonas más adecuadas en términos de ph del suelo para
el desarrollo del jitomate?
Las zonas adecuadas para jitomate son las que presentan pH 6,5.