rm(list=ls())
##it is assummed that you already installed terra, sf, and leaflet
##paquetes= c('XML','rgdal', 'dplyr','devtools', 'stars', 'mapview', 'RColorBrewer')
##install.packages(paquetes)
##devtools::install_github("JoshOBrien/gdalUtilities")
##devtools::install_github('gearslaboratory/gdalUtils')
library(XML)
## Warning: package 'XML' was built under R version 4.4.3
library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(gdalUtilities)
## Warning: package 'gdalUtilities' was built under R version 4.4.3
##
## Adjuntando el paquete: 'gdalUtilities'
## The following object is masked from 'package:sf':
##
## gdal_rasterize
library(terra)
## Warning: package 'terra' was built under R version 4.4.3
## terra 1.8.54
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)
## Warning: package 'leaflet' was built under R version 4.4.3
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"
list.files("C:/Users/yorhl/Downloads/DATOS")
## [1] "aguas de guaviare.gpkg"
## [2] "Base agrícola 2023.csv"
## [3] "cities.gpkg"
## [4] "citis.gpkg"
## [5] "COL_AGUA"
## [6] "COL_ELEVACION"
## [7] "COL_msk_alt_tif"
## [8] "COL_rVIAS"
## [9] "colores poblacion.gpkg"
## [10] "DCD-POBLA.csv"
## [11] "dps-2.csv"
## [12] "Evaluaciones_Agropecuarias_Municipales_EVA.csv"
## [13] "guaviare"
## [14] "Guaviare.gpkg"
## [15] "Guaviare.tif"
## [16] "mapa de color.gpkg"
## [17] "MGN2023_DPTO_POLITICO"
## [18] "MGN2023_MPIO_POLITICO"
## [19] "mipo.gpkg"
## [20] "MPM"
## [21] "Mun-2023.csv"
## [22] "puntos_poblacion.gpkg"
## [23] "simplemaps_worldcities_basicv1.77"
## [24] "soc.tif"
## [25] "soil_GUAVIARE.Rmd"
## [26] "Vias.gpkg"
## [27] "worldcities.csv"
## [28] "ytree.gpkg"
(Guaviare = st_read("C:/Users/yorhl/Downloads/DATOS/Guaviare.gpkg"))
## Multiple layers are present in data source C:\Users\yorhl\Downloads\DATOS\Guaviare.gpkg, reading layer `mipio'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `mipio' from data source
## `C:\Users\yorhl\Downloads\DATOS\Guaviare.gpkg' using driver `GPKG'
## Simple feature collection with 4 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.66391 ymin: 0.6554126 xmax: -69.99511 ymax: 2.924987
## Geodetic CRS: MAGNA-SIRGAS
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
Guaviare_igh <- st_transform(Guaviare, igh)
(bbox <- st_bbox(Guaviare_igh))
## xmin ymin xmax ymax
## -8201803.20 72960.19 -7794310.65 325608.03
## ul means upper left
## lr means lower right
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
## xmin ymax xmax ymin
## -8201803.20 325608.03 -7794310.65 72960.19
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc.tif'
file = "soc.tif"
# UNCOMMENT FOR FIRST TIME USE
# THEN, COMMENT IT BACK
#gdal_translate(paste0(sg_url,datos), file ,
# tr=c(250,250),
# projwin=bb,
# projwin_srs =igh)
(Guaviare_soc <- terra::rast(file)/10)
## class : SpatRaster
## size : 837, 885, 1 (nrow, ncol, nlyr)
## resolution : 0.002259887, 0.002389486 (x, y)
## extent : -73.0659, -71.0659, 0.9126, 2.9126 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : soc
## name : soc
## min value : 0.0
## max value : 8.1
terra::hist(Guaviare_soc)
summary(Guaviare_soc)
## Warning: [summary] used a sample
## soc
## Min. :0.000
## 1st Qu.:4.100
## Median :4.400
## Mean :4.412
## 3rd Qu.:4.700
## Max. :7.900
(names(Guaviare_soc) <- "soc")
## [1] "soc"
valores <- values(Guaviare_soc, na.rm=TRUE)
orangecyan <- colorNumeric(c("orange","yellow2", "darkseagreen", "cyan" ), valores,
na.color = "transparent")
leaflet::leaflet(Guaviare) %>%
addTiles() %>%
setView(-74, 6, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.5, fillOpacity = 0.10,
popup = paste("Municipio: ", Guaviare$MPIO_CNMBR)) %>%
addRasterImage(Guaviare_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'
set.seed(123456)
# Conversion of SOC coordinate reference system
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(Guaviare_soc, geog))
## class : SpatRaster
## size : 861, 861, 1 (nrow, ncol, nlyr)
## resolution : 0.002321978, 0.002321978 (x, y)
## extent : -73.0659, -71.06668, 0.9133765, 2.9126 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : soc
## min value : 0.00000
## max value : 7.99161
(samples <- terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))
## class : SpatVector
## geometry : points
## dimensions : 2000, 1 (geometries, attributes)
## extent : -73.06242, -71.06784, 0.9145375, 2.911439 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : soc
## type : <num>
## values : 4.654
## 4.504
## 4.775
(muestras <- sf::st_as_sf(samples))
nmuestras <- na.omit(muestras)
longit <- st_coordinates(nmuestras)[,1]
latit <- st_coordinates(nmuestras)[,2]
soc <- nmuestras$soc
summary(soc)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 4.112 4.376 4.396 4.679 6.791
length(soc)
## [1] 2000
id <- seq(1,2000,1)
(sitios <- data.frame(id, longit, latit, soc))
m <- leaflet() %>%
addTiles() %>%
addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m # Print the map
# change path and filename to match your directories
# as well as your department name
sf::st_write(nmuestras, "C:/Users/yorhl/Downloads/DATOS/Guaviare.gpkg", driver = "GPKG",
append=FALSE)
## Deleting layer `Guaviare' using driver `GPKG'
## Writing layer `Guaviare' to data source
## `C:/Users/yorhl/Downloads/DATOS/Guaviare.gpkg' using driver `GPKG'
## Writing 2000 features with 1 fields and geometry type Point.
sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-3 leaflet_2.2.2 dplyr_1.1.4
## [4] terra_1.8-54 gdalUtilities_1.2.5 sf_1.0-21
## [7] XML_3.99-0.18
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.9 utf8_1.2.4 generics_0.1.3 class_7.3-22
## [5] KernSmooth_2.23-24 digest_0.6.37 magrittr_2.0.3 evaluate_1.0.1
## [9] grid_4.4.1 fastmap_1.2.0 jsonlite_1.8.9 e1071_1.7-16
## [13] DBI_1.2.3 fansi_1.0.6 crosstalk_1.2.1 scales_1.3.0
## [17] codetools_0.2-20 jquerylib_0.1.4 cli_3.6.3 rlang_1.1.4
## [21] units_0.8-5 munsell_0.5.1 cachem_1.1.0 yaml_2.3.10
## [25] tools_4.4.1 colorspace_2.1-1 vctrs_0.6.5 R6_2.5.1
## [29] png_0.1-8 proxy_0.4-27 lifecycle_1.0.4 classInt_0.4-10
## [33] htmlwidgets_1.6.4 pkgconfig_2.0.3 pillar_1.9.0 bslib_0.8.0
## [37] glue_1.8.0 Rcpp_1.0.13 xfun_0.49 tibble_3.2.1
## [41] tidyselect_1.2.1 rstudioapi_0.17.1 knitr_1.49 farver_2.1.2
## [45] htmltools_0.5.8.1 rmarkdown_2.29 compiler_4.4.1