Este cuaderno explica cómo descargar datos raster de SoilGrids para el departamento de Cundinamarca, el proceso se ilustra utilizando la variable de carbono orgánico del suelo. SoilGrids es un sistema para la cartografía digital global del suelo que hace uso de la información global del perfil del suelo y los datos de covarianza para modelar la distribución espacial de las propiedades del suelo en todo el mundo.
limpiamos la consola
rm(list=ls())
Y luego instalamos y llamamos las librerias necesarias:
library(XML)
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(gdalUtilities)
##
## Adjuntando el paquete: 'gdalUtilities'
## The following object is masked from 'package:sf':
##
## gdal_rasterize
library(terra)
## terra 1.8.15
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(leaflet)
library(RColorBrewer)
url = "https://files.isric.org/soilgrids/latest/data/"
# basic strings
voi = "soc" # soil organic carbon
depth = "15-30cm"
quantile = "mean"
# concatenate the strings
(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"
Leemos nuestro GeoPackage de los datos de municipios de Cundinamarca previamente descargados del DANE.
(cun = sf:: st_read("Cundinamarca/cundinamarca mun_Cun.gpkg"))
## Reading layer `mgn_adm_mpio_grafico' from data source
## `C:\Users\diazo\OneDrive\Documentos\GEOMATICA\PROYECTO_5\cundinamarca\cundinamarca mun_Cun.gpkg'
## using driver `GPKG'
## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 25 001 25001 CUNDINAMARCA AGUA DE DIOS
## 2 25 019 25019 CUNDINAMARCA ALBÁN
## 3 25 035 25035 CUNDINAMARCA ANAPOIMA
## 4 25 040 25040 CUNDINAMARCA ANOLAIMA
## 5 25 053 25053 CUNDINAMARCA ARBELÁEZ
## 6 25 086 25086 CUNDINAMARCA BELTRÁN
## 7 25 095 25095 CUNDINAMARCA BITUIMA
## 8 25 099 25099 CUNDINAMARCA BOJACÁ
## 9 25 120 25120 CUNDINAMARCA CABRERA
## 10 25 123 25123 CUNDINAMARCA CACHIPAY
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 78 Noviembre 29 de 1963 MUNICIPIO 85.95871 2023
## 2 Ordenanza 19 del 22 de Mayo de 1903 MUNICIPIO 50.83510 2023
## 3 Dispocisión de Septiembre 1 de 1864 MUNICIPIO 123.95549 2023
## 4 1882 MUNICIPIO 120.95146 2023
## 5 Decreto 32 de Enero 16 de 1886 MUNICIPIO 142.45716 2023
## 6 Ordenanza 5 del 12 de Noviembre de 1853 MUNICIPIO 176.85320 2023
## 7 1772 MUNICIPIO 61.04708 2023
## 8 1700 MUNICIPIO 102.31988 2023
## 9 Ordenanza 79 Noviembre 29 de 1963 MUNICIPIO 421.83959 2023
## 10 Ordenanza 6 Noviembre 26 de 1982 MUNICIPIO 53.45262 2023
## shape_Leng shape_Area geom
## 1 0.4963557 0.007002449 MULTIPOLYGON (((-74.67275 4...
## 2 0.5162142 0.004144411 MULTIPOLYGON (((-74.47852 4...
## 3 0.6705565 0.010100758 MULTIPOLYGON (((-74.55201 4...
## 4 0.5280079 0.009859268 MULTIPOLYGON (((-74.43043 4...
## 5 0.8065843 0.011603784 MULTIPOLYGON (((-74.42848 4...
## 6 0.9093670 0.014413273 MULTIPOLYGON (((-74.76139 4...
## 7 0.4907957 0.004976541 MULTIPOLYGON (((-74.51062 4...
## 8 0.4979169 0.008339597 MULTIPOLYGON (((-74.32819 4...
## 9 1.0188051 0.034347363 MULTIPOLYGON (((-74.46463 4...
## 10 0.4211104 0.004356744 MULTIPOLYGON (((-74.40513 4...
Reproyectamos nuestras capas
igh= '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
cun_igh <- sf::st_transform(cun, igh)
Conseguimos el cuadro delimitador del area que nos interesa
(bbox <- sf::st_bbox(cun_igh))
## xmin ymin xmax ymax
## -8344613.5 415236.1 -8142417.0 649800.5
Utilizamos los datos bbox para definir los límites de la caja límite tal y como los utiliza la librería GDAL.
ulx = bbox$xmin
uly = bbox$ymax
lrx = bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
## xmin ymax xmax ymin
## -8344613.5 649800.5 -8142417.0 415236.1
Usamos la función gdal_translate para descargar la capa vrt
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "Soilsgrids.org.tif"
Dividimos los valores de las predicciones (en unidades cartografiadas) por los valores de la columna Factor de conversión
(cun_soc <- terra::rast(file)/10)
## class : SpatRaster
## dimensions : 837, 885, 1 (nrow, ncol, nlyr)
## resolution : 0.002259887, 0.002389486 (x, y)
## extent : -74.9136, -72.9136, 3.8256, 5.8256 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : Soilsgrids.org
## name : Soilsgrids.org
## min value : 0.0
## max value : 92.5
Realizamos un histograma
terra::hist(cun_soc)
summary(cun_soc)
## Warning: [summary] used a sample
## Soilsgrids.org
## Min. : 0.00
## 1st Qu.:17.20
## Median :22.30
## Mean :21.12
## 3rd Qu.:25.80
## Max. :84.10
Cambiamos el nombre del atributo SOC
(names(cun_soc) <- "soc")
## [1] "soc"
Se ponen los nombres valores de SOC en una variable
valores <- values(cun_soc, na.rm=TRUE)
Creamos paleta de colores
orangecyan <- colorNumeric(c("orange", "yellow2", "darkseagreen", "cyan"), valores,
na.color = "transparent")
y hacemos la parcela
leaflet::leaflet(cun) %>%
addTiles() %>%
setView(-74, 6, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.5, fillOpacity = 0.10,
popup = paste("Municipio: ",cun$mpio_cnmbr)) %>%
addRasterImage(cun_soc, colors = "Spectral", opacity = 0.8) %>%
addLegend(pal = orangecyan, values = valores, title = "Contenido organico del suelo (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'