Este cuaderno explica brevemente como descargar datos raster de SoilGrigds usando R. Este proceso se lleva a cabo con datos de carbono organico en el suelo como una variable de interes, claramente el codigo se puede replicar usando una variable distinta.
Posterior a la descarga de los datos procederemos a ilustrar dos metodos de interpolacion espacial: Inverse Dinstance Weighted (IDW) y Ordinary Kriging (OK). IDW es una tecnica deterministica y OK es una tecnica probabilistica. Ambas técnicas se utilizan aquí para obtener una superficie continua de carbono orgánico del suelo (SOC) a una profundidad de 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m.
limpiar el workspace
rm(list=ls())
Cargamos las librerias necesarias
Definimos el link que corresponde al sitio web de ISRIC
url = "https://files.isric.org/soilgrids/latest/data/"
Entonces, creamos algunos objetos indicando que necesitamos de ISRIC
#strings basicos
voi = "soc" #Soil organic carbon
depth = "15-30cm"
quantile = "mean"
#concatenar los strings
(variable = paste(url, voi, sep = ""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc"
Ahora definimos la propiedad del suelo que queremos descargar
(layer = paste(variable,depth,quantile, sep = "_")) #capa de interes
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"
asi, la verificación VRT esta completa
(vrt_layer = paste(layer, '.vrt', sep = "" ))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"
Lee un mapa vector con el area de interes previamente descargado
#Cambia la ruta y el nombre del archivo para que coincidan con tus datos.
(stder = st_read("C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/nariñito.gpkg"))
## Reading layer `municipios' from data source
## `C:\Users\vocav\Documents\GB2-2025 nicolin-20250723T070234Z-1-001\GB2-2025 nicolin\GB2-2025\GB2\PR4\nariñito.gpkg'
## using driver `GPKG'
## Simple feature collection with 64 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 64 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 52 001 52001 NARIÑO PASTO
## 2 52 019 52019 NARIÑO ALBÁN
## 3 52 022 52022 NARIÑO ALDANA
## 4 52 036 52036 NARIÑO ANCUYA
## 5 52 051 52051 NARIÑO ARBOLEDA
## 6 52 079 52079 NARIÑO BARBACOAS
## 7 52 083 52083 NARIÑO BELÉN
## 8 52 110 52110 NARIÑO BUESACO
## 9 52 203 52203 NARIÑO COLÓN
## 10 52 207 52207 NARIÑO CONSACÁ
## mpio_crslc mpio_tipo mpio_narea mpio_nano shape_Leng
## 1 1540 MUNICIPIO 1099.37825 2023 1.8748401
## 2 Ordenanza 41 de 1903 MUNICIPIO 38.61009 2023 0.3402693
## 3 Ordenanza 11 de 1911 MUNICIPIO 47.64055 2023 0.4520361
## 4 1544 MUNICIPIO 68.90636 2023 0.3736985
## 5 1859 MUNICIPIO 60.18576 2023 0.3754062
## 6 1916 MUNICIPIO 2732.93328 2023 3.6991487
## 7 Ordenanza 53 Noviembre 29 de 1985 MUNICIPIO 41.84541 2023 0.3732840
## 8 1899 MUNICIPIO 635.96083 2023 1.2292312
## 9 Ordenanza 37 de 1921 MUNICIPIO 61.75053 2023 0.4592866
## 10 1900 MUNICIPIO 119.35364 2023 0.4650471
## shape_Area geom
## 1 0.089062266 MULTIPOLYGON (((-77.31099 1...
## 2 0.003129126 MULTIPOLYGON (((-77.06531 1...
## 3 0.003855327 MULTIPOLYGON (((-77.69156 0...
## 4 0.005578855 MULTIPOLYGON (((-77.50924 1...
## 5 0.004877181 MULTIPOLYGON (((-77.13479 1...
## 6 0.220963473 MULTIPOLYGON (((-77.79456 1...
## 7 0.003391678 MULTIPOLYGON (((-77.07227 1...
## 8 0.051533090 MULTIPOLYGON (((-77.23516 1...
## 9 0.005005108 MULTIPOLYGON (((-77.04473 1...
## 10 0.009664908 MULTIPOLYGON (((-77.50158 1...
Nota el CRS de stder. Como las capas de ISRIC usan la proyección Homolosine, necesitamos reproyectar nuestra capa usando la biblioteca sf.
repro = '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
stder_reproj = st_transform(stder, repro)
Ahora definimos los limites de nuestra area de interes
(box = st_bbox(stder_reproj))
## xmin ymin xmax ymax
## -8796326.48 40225.08 -8554103.93 298770.17
Ahora, usemos los datos de bbox para definir los límites de la caja de delimitación tal como los usa la biblioteca GDAL. Por cierto, esta es una de las partes más complicadas de usar GDAL.
ulx = box$xmin
uly = box$ymax
lrx = box$xmax
lry = box$ymin
(bbx = c(ulx, uly, lrx,lry))
## xmin ymax xmax ymin
## -8796326.48 298770.17 -8554103.93 40225.08
Ahora, podemos usar la función gdal_translate para descargar la capa VRT. Primero, definamos dónde guardar los datos:
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/carbono_orgánico_suelo.tif"
Ten en cuenta que la siguiente tarea puede tardar varios minutos. ¡Ejecútala solo una vez!
Las unidades de la capa SOC no son las más comunes, para facilitar su uso vamos a usar un factor de conversión del cual obtenemos unidades familiares.
(stder_soc <- terra::rast(file)/10)
## class : SpatRaster
## size : 972, 963, 1 (nrow, ncol, nlyr)
## resolution : 0.002260748, 0.002389609 (x, y)
## extent : -79.0102, -76.8331, 0.3613, 2.684 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : carbono_orgánico_suelo
## name : carbono_orgánico_suelo
## min value : 0.0
## max value : 262.3
##4. Exploración de datos
Un histograma puede ayudar con esto
## añade un titulo significativo y las unidades de SOC
(media_soc = mean(values(stder_soc), na.rm = T))
## [1] 71.41871
par(bg = "white")
options(scipen = 10)
hist(values(stder_soc), main= "Contenido de carbono orgánico del suelo en la fracción de tierra fina",
xlab= "g/kg", xlim = c(0, 120), col="lightblue", breaks=20)
Un resumen también es útil.
summary(stder_soc)
## Warning: [summary] used a sample
## carbono_orgánico_suelo
## Min. : 0.0
## 1st Qu.: 35.6
## Median : 56.4
## Mean : 71.4
## 3rd Qu.: 98.3
## Max. :251.8
Ahora cambiamos el nombre del atributo SOC
(names(stder_soc) <- "soc")
## [1] "soc"
Es hora de graficar. ### Ponemos los valores de SOC en una variable
val = values(stder_soc, na.rm = T)
creamos una paleta de colores
orangecyan = colorNumeric(c("orange","yellow2", "darkseagreen", "cyan" ), val,
na.color = "transparent")
Por ultimo la grafica.
#Grafica un mapa interactivo
stder <- st_transform(stder, crs = 4326)
leaflet::leaflet(stder) %>%
addTiles() %>%
setView(-74, 6, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.5, fillOpacity = 0.10,
) %>%
addRasterImage(stder_soc, colors ="Spectral", opacity = 0.8) %>%
addLegend(pal = orangecyan, values = val, title = "Soil Organic Carbon (SOC) [g/kg]")