Este cuaderno presenta el proceso de descarga de datos de suelos proporcionados por SoilGrids, una plataforma de modelado digital de suelos a nivel mundial. Los datos utilizados serán de carbono orgánico del suelo.
Antes de comenzar, se debe limpiar el espacio de trabajo y asegurarse de que las librerías necesarias estén instaladas y cargadas:
rm(list=ls())
library(XML)
library(sf)
library(gdalUtilities)
library(terra)
library(dplyr)
library(leaflet)
library(RColorBrewer)
SoilGrids proporciona acceso a datos a través de varios servicios. Utilizaremos WebDAV para descargar los datos:
url = "https://files.isric.org/soilgrids/latest/data/"
voi = "soc"
depth = "15-30cm"
quantile = "mean"
(variable = paste(url, voi, sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc"
(layer = paste(variable,depth,quantile, sep="_"))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"
(vrt_layer = paste(layer, '.vrt', sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"
Usaremos un archivo GeoPackage con los municipios de Santander (Colombia):
Cargar los datos de la región de interés, definir el sistema de coordenadas proyectadas, las coordenadas de los limites, obtener la caja delimitadora.
(stder = st_read("datos/dpto_santander.gpkg"))
## Reading layer `dpto_santandercorr' from data source
## `C:\Users\lausi\OneDrive\Documents\2024-2\GEOMATICA\GB02\Proyecto R\datos\dpto_santander.gpkg'
## using driver `GPKG'
## Simple feature collection with 87 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 87 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp mpio_cnmbr area_km2
## 1 68 001 68001 BUCARAMANGA 153.74283
## 2 68 013 68013 AGUADA 75.22441
## 3 68 020 68020 ALBANIA 166.21369
## 4 68 051 68051 ARATOCA 169.73316
## 5 68 077 68077 BARBOSA 46.71246
## 6 68 079 68079 BARICHARA 137.24491
## 7 68 081 68081 BARRANCABERMEJA 1326.79871
## 8 68 092 68092 BETULIA 431.18573
## 9 68 101 68101 BOLÍVAR 1010.10956
## 10 68 121 68121 CABRERA 65.56073
## geom
## 1 MULTIPOLYGON (((-73.15042 7...
## 2 MULTIPOLYGON (((-73.56261 6...
## 3 MULTIPOLYGON (((-73.73616 5...
## 4 MULTIPOLYGON (((-72.98158 6...
## 5 MULTIPOLYGON (((-73.58988 5...
## 6 MULTIPOLYGON (((-73.22126 6...
## 7 MULTIPOLYGON (((-73.6939 7....
## 8 MULTIPOLYGON (((-73.53993 7...
## 9 MULTIPOLYGON (((-74.50132 6...
## 10 MULTIPOLYGON (((-73.25696 6...
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
stder_igh <- st_transform(stder, igh)
(bbox <- st_bbox(stder_igh))
## xmin ymin xmax ymax
## -8313449.6 635360.0 -8089703.8 906698.4
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
## xmin ymax xmax ymin
## -8313449.6 906698.4 -8089703.8 635360.0
Definir el archivo de salida y la URL de los datos.
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"
Descargar los datos utilizando GDAL.
gdal_translate(paste0(sg_url,datos), file ,
tr=c(250,250),
projwin=bb,
projwin_srs =igh)
## Warning in CPL_gdaltranslate(source, destination, options, oo, config_options,
## : GDAL Error 3: Cannot read 5933045 bytes
Cargar el raster descargado y convertir unidades.
(stder_soc <- terra::rast(file)/10)
## class : SpatRaster
## dimensions : 1085, 895, 1 (nrow, ncol, nlyr)
## resolution : 250, 250 (x, y)
## extent : -8313500, -8089750, 635500, 906750 (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s) : memory
## varname : soc_igh_15_30
## name : soc_igh_15_30
## min value : 8.7
## max value : 130.3
terra::hist(stder_soc)
summary(stder_soc)
## Warning: [summary] used a sample
## soc_igh_15_30
## Min. : 9.50
## 1st Qu.: 25.10
## Median : 33.20
## Mean : 35.81
## 3rd Qu.: 44.40
## Max. :130.20
## NA's :1588
(names(stder_soc) <- "soc")
## [1] "soc"
valores <- values(stder_soc, na.rm=TRUE)
Definir una paleta de colores para la visualización.
orangecyan <- colorNumeric(c("orange","yellow2", "darkseagreen", "cyan" ), valores,
na.color = "transparent")
Crear un mapa interactivo con Leaflet
leaflet::leaflet(stder) %>%
addTiles() %>%
setView(-74, 6, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.5, fillOpacity = 0.10,
popup = paste("Municipio: ", stder$MPIO_CNMBR)) %>%
addRasterImage(stder_soc, colors ="Spectral", opacity = 0.8) %>%
addLegend(pal = orangecyan, values = valores, title = "Soil Organic Carbon (SOC) [g/kg]")
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'