rm(list = ls())
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="_")) # layer of interest
## [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"
(stder = st_read("D:/GB2_UNAL/Proyecto4/Datos/Mun_meta.gpkg"))
## Reading layer `municipios_colombia' from data source
## `D:\GB2_UNAL\Proyecto4\Datos\Mun_meta.gpkg' using driver `GPKG'
## Simple feature collection with 29 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.93335 ymin: 1.617732 xmax: -71.07753 ymax: 4.899101
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 29 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.93335 ymin: 1.617732 xmax: -71.07753 ymax: 4.899101
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 50 001 50001 META VILLAVICENCIO
## 2 50 006 50006 META ACACÍAS
## 3 50 110 50110 META BARRANCA DE UPÍA
## 4 50 124 50124 META CABUYARO
## 5 50 150 50150 META CASTILLA LA NUEVA
## 6 50 223 50223 META CUBARRAL
## 7 50 226 50226 META CUMARAL
## 8 50 245 50245 META EL CALVARIO
## 9 50 251 50251 META EL CASTILLO
## 10 50 270 50270 META EL DORADO
## mpio_crslc mpio_tipo
## 1 1840 MUNICIPIO
## 2 Ordenanza 23 de 1960 MUNICIPIO
## 3 Ordenanza 5 de Octubre 16 de 1990 MUNICIPIO
## 4 Decreto 47 de Septiembre 3 de 1912 MUNICIPIO
## 5 Ordenanza 08 de Julio 7 de 1961 MUNICIPIO
## 6 Ordenanza 23 de Noviembre 28 de 1960 MUNICIPIO
## 7 Decreto 2543 de 1955 MUNICIPIO
## 8 Decreto 2543 de 1955 MUNICIPIO
## 9 Ordenanza 25 de Febrero de 1976 MUNICIPIO
## 10 Ordenanza 44 Noviembre 24 de 1992. Decreto 2129 de Diciembr MUNICIPIO
## mpio_narea mpio_nano shape_Leng shape_Area geom
## 1 1285.8794 2023 2.0397750 0.104714061 MULTIPOLYGON (((-73.69288 4...
## 2 1123.7379 2023 2.0742072 0.091507447 MULTIPOLYGON (((-73.80268 4...
## 3 403.8228 2023 1.3120271 0.032894731 MULTIPOLYGON (((-73.01607 4...
## 4 913.6625 2023 2.0088484 0.074403262 MULTIPOLYGON (((-73.08865 4...
## 5 514.6860 2023 1.3493293 0.041899616 MULTIPOLYGON (((-73.7252 3....
## 6 1157.8652 2023 2.0750915 0.094269993 MULTIPOLYGON (((-74.20143 4...
## 7 625.0714 2023 1.8723272 0.050907259 MULTIPOLYGON (((-73.55407 4...
## 8 270.3071 2023 0.8271179 0.022020904 MULTIPOLYGON (((-73.74633 4...
## 9 568.2181 2023 1.1179992 0.046251583 MULTIPOLYGON (((-73.92004 3...
## 10 117.3411 2023 0.4364058 0.009552168 MULTIPOLYGON (((-73.83554 3...
stder %>% as_tibble()
## # A tibble: 29 × 12
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr mpio_crslc mpio_tipo
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 50 001 50001 META VILLAVICENC… 1840 MUNICIPIO
## 2 50 006 50006 META ACACÍAS Ordenanza… MUNICIPIO
## 3 50 110 50110 META BARRANCA DE… Ordenanza… MUNICIPIO
## 4 50 124 50124 META CABUYARO Decreto 4… MUNICIPIO
## 5 50 150 50150 META CASTILLA LA… Ordenanza… MUNICIPIO
## 6 50 223 50223 META CUBARRAL Ordenanza… MUNICIPIO
## 7 50 226 50226 META CUMARAL Decreto 2… MUNICIPIO
## 8 50 245 50245 META EL CALVARIO Decreto 2… MUNICIPIO
## 9 50 251 50251 META EL CASTILLO Ordenanza… MUNICIPIO
## 10 50 270 50270 META EL DORADO Ordenanza… MUNICIPIO
## # ℹ 19 more rows
## # ℹ 5 more variables: mpio_narea <dbl>, mpio_nano <int>, shape_Leng <dbl>,
## # shape_Area <dbl>, geom <MULTIPOLYGON [°]>
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
## -8345138.8 180085.1 -7916386.8 545365.5
ulx = bbox$xmin
uly = bbox$ymax
lrx = bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
## xmin ymax xmax ymin
## -8345138.8 545365.5 -7916386.8 180085.1
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"
#gdal_translate(paste0(sg_url,datos), file ,
# tr=c(250,250),
# projwin=bb,
# projwin_srs =igh)
(stder_soc <- terra::rast(file)/10)
## class : SpatRaster
## dimensions : 1461, 1715, 1 (nrow, ncol, nlyr)
## resolution : 250, 250 (x, y)
## extent : -8345250, -7916500, 180250, 545500 (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s) : memory
## varname : soc_igh_15_30
## name : soc_igh_15_30
## min value : 7.4
## max value : 116.5
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)

summary(stder_soc)
## Warning: [summary] used a sample
## soc_igh_15_30
## Min. : 7.70
## 1st Qu.: 12.10
## Median : 14.90
## Mean : 20.89
## 3rd Qu.: 21.20
## Max. :116.50
## NA's :795
(names(stder_soc) <- "soc")
## [1] "soc"
valores <- values(stder_soc, na.rm=TRUE)
orangecyan <- colorNumeric(c("Spectral"), valores,
na.color = "transparent")
# plot an interactive map
leaflet::leaflet(stder) %>%
addTiles() %>%
setView(-72, 3, 7) %>%
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 = orangecyan, 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'