rm(list=ls())
library(XML)
library(sf)
library(terra)
library(dplyr)
library(leaflet)
library(RColorBrewer)
library(gdalUtilities)
url = "https://files.isric.org/soilgrids/latest/data/"
## [1] "https://files.isric.org/soilgrids/latest/data/soc"
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"
csnr = st_read("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/CASANARE.gpkg")
## Reading layer `municipios' from data source 
##   `C:\Users\Janus\Desktop\GB2\proyecto3\datos\departamentos\CASANARE\CASANARE.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 19 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
## Geodetic CRS:  MAGNA-SIRGAS
igh= "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs"
csnr_igh  <- st_transform(csnr, igh)
(bbox <- st_bbox(csnr_igh))
##       xmin       ymin       xmax       ymax 
## -8145533.4   477279.6 -7792746.2   706445.8
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
##       xmin       ymax       xmax       ymin 
## -8145533.4   706445.8 -7792746.2   477279.6
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)
(csnr_soc <- terra::rast(file)/10)
## class       : SpatRaster 
## size        : 917, 1413, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8145750, -7792500, 477250, 706500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source(s)   : memory
## varname     : soc_igh_15_30 
## name        : soc_igh_15_30 
## min value   :           6.6 
## max value   :         118.9
terra::hist(csnr_soc)
## Warning: [hist] a sample of 77% of the cells was used (of which 1% was NA)

summary(csnr_soc)
## Warning: [summary] used a sample
##  soc_igh_15_30   
##  Min.   :  6.70  
##  1st Qu.: 13.50  
##  Median : 16.10  
##  Mean   : 20.73  
##  3rd Qu.: 20.80  
##  Max.   :114.80  
##  NA's   :1179
(names(csnr_soc) <- "soc")
## [1] "soc"
valores <- values(csnr_soc, na.rm=TRUE)
orangecyan <- colorNumeric(c("orange","yellow2", "darkseagreen", "cyan" ), valores,
  na.color = "transparent")
leaflet::leaflet(csnr) %>% 
  addTiles() %>% 
  setView(-72.0, 5.3, 9) %>% 
   addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.5, fillOpacity = 0.10,
    popup = paste("Municipio: ", csnr$MPIO_CNMBR)) %>%
  addRasterImage(csnr_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)
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(csnr_soc, geog))
## class       : SpatRaster 
## size        : 934, 1489, 1  (nrow, ncol, nlyr)
## resolution  : 0.002205488, 0.002205488  (x, y)
## extent      : -73.09924, -69.81527, 4.286671, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :        soc 
## min value   :   6.961041 
## max value   : 116.768845
(samples <- terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2000, 1  (geometries, attributes)
##  extent      : -73.09594, -69.81638, 4.287774, 6.345495  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      : 10.55
##                23.49
##                15.01
(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. 
##   8.372  13.600  16.117  20.679  20.409  93.331
length(soc)
## [1] 1925
id <- seq(1,1925,1)
(sitios <- data.frame(id, longit, latit, soc))
m <- leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m
sf::st_write(nmuestras, "C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/soc_csnr.gpkg", driver = "GPKG",
             append=FALSE)
## Deleting layer `soc_csnr' using driver `GPKG'
## Writing layer `soc_csnr' to data source 
##   `C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/soc_csnr.gpkg' using driver `GPKG'
## Writing 1925 features with 1 fields and geometry type Point.
sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## 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] gdalUtilities_1.2.5 RColorBrewer_1.1-3  leaflet_2.2.2      
## [4] dplyr_1.1.4         terra_1.8-54        sf_1.0-21          
## [7] XML_3.99-0.18      
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_2.0.0     compiler_4.5.0     tidyselect_1.2.1   Rcpp_1.0.14       
##  [5] jquerylib_0.1.4    png_0.1-8          scales_1.4.0       yaml_2.3.10       
##  [9] fastmap_1.2.0      R6_2.6.1           generics_0.1.4     classInt_0.4-11   
## [13] knitr_1.50         htmlwidgets_1.6.4  tibble_3.2.1       units_0.8-7       
## [17] DBI_1.2.3          bslib_0.9.0        pillar_1.10.2      rlang_1.1.6       
## [21] cachem_1.1.0       xfun_0.52          sass_0.4.10        cli_3.6.5         
## [25] magrittr_2.0.3     class_7.3-23       crosstalk_1.2.1    digest_0.6.37     
## [29] grid_4.5.0         rstudioapi_0.17.1  lifecycle_1.0.4    vctrs_0.6.5       
## [33] KernSmooth_2.23-26 proxy_0.4-27       evaluate_1.0.3     glue_1.8.0        
## [37] farver_2.1.2       codetools_0.2-20   e1071_1.7-16       rmarkdown_2.29    
## [41] tools_4.5.0        pkgconfig_2.0.3    htmltools_0.5.8.1

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.