Interpolación espacial del carbono orgánico del suelo en el departamento de Bolivar
Danna Castañeda, Manuel Ignacio Vargas
2023-06-22
1. Introducción
En el desarrollo de este cuaderno se va a representar dos estrategias de interpolación Distancia ponderada inversa (IDW) y Kriging ordinario (OK). Donde IDW es una técnica determinista y OK es probabilístico para la obtención de una superficie continua del carbono orgánico a 15-30cm, cabe resaltar que los datos aquí usados se extrajeron de Soilgrids
#Las librerias que se usaron fueron las siguientes
library(sp)## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0library(terra)## terra 1.7.37library(sf)## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(stars)## Loading required package: abindlibrary(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)##
## Attaching package: '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, unionlibrary(curl)## Using libcurl 7.84.0 with Schannel#rm(list=ls())
#Controlador para la subida del TIFF
h <- new_handle()
handle_setopt(h, http_version = 2)2. Leer los datos de entrada
A continuación se va a leer el archivo TIFF del carbono orgánico del suelo a 15-30cm usando la libreria terra que proporciona métodos para el análisis de datos espaciales con datos raster y vectoriales
archivo <- ("SOC_15-30.TIF")
(soc <- rast(archivo))## class : SpatRaster
## dimensions : 1599, 897, 1 (nrow, ncol, nlyr)
## resolution : 0.002261584, 0.002390448 (x, y)
## extent : -75.79764, -73.769, 6.979202, 10.80153 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : SOC_15-30.tif
## name : SOC_15-30
## min value : 0
## max value : 1277Ahora se convierten los datos SOC, es decir de carbono orgánico a porcentaje
soc.perc<-soc/10Ahora necesitamos reproyectar dichos datos CRS a WGS84 CRS, donde CRS es un sistema de números que definen ubicaciones sobre la superficie de la tierra, mientras que el WGS84 es un sistema de coordenadas geográficas mundial que permite localizar cualquier punto de la Tierra
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))## class : SpatRaster
## dimensions : 1619, 860, 1 (nrow, ncol, nlyr)
## resolution : 0.002360243, 0.002360243 (x, y)
## extent : -75.79764, -73.76783, 6.980295, 10.80153 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : SOC_15-30
## min value : 0.0000
## max value : 139.7799Vamos a convertir la capa ** SpatRaster ** en un objeto ** stars ** de manera que ahora nos va a permitir visualizarlo
stars.soc = st_as_stars(geog.soc)
m <- leaflet() %>%
addTiles() %>%
leafem:::addGeoRaster(
stars.soc,
opacity = 0.8,
colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
domain = 8:130)
) 3. Muestreo del mundo
A continuación se obtiene una muestra de aproximadamente 500 lugares de los datos del mundo real para el departamento de Bolivar, estas muestras tomadas de forma aleatoria
set.seed(123456)
# Muestreo aleatorio de 500 puntos
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))## class : SpatVector
## geometry : points
## dimensions : 500, 1 (geometries, attributes)
## extent : -75.79646, -73.77373, 6.981475, 10.78855 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : SOC_15-30
## type : <num>
## values : 28.06
## 27.69
## 35.99Esta muestra posee 500 geometrias y un atributo, es un SpatVector y arroja tres valores de SOC y posteriormente convertir el ** spatVector ** en un objeto ** simple feature **
(muestras <- sf::st_as_sf(samples))## Simple feature collection with 500 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.79646 ymin: 6.981475 xmax: -73.77373 ymax: 10.78855
## Geodetic CRS: GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
## First 10 features:
## SOC_15-30 geometry
## 1 28.05957 POINT (-74.25286 8.244205)
## 2 27.68818 POINT (-73.96963 8.28669)
## 3 35.98964 POINT (-75.41174 8.088429)
## 4 14.12467 POINT (-74.35435 10.24333)
## 5 16.05506 POINT (-75.67609 7.847685)
## 6 43.62112 POINT (-74.95621 8.211162)
## 7 14.47324 POINT (-74.93261 9.896376)
## 8 36.97206 POINT (-74.19149 7.559735)
## 9 17.00693 POINT (-73.99323 9.237868)
## 10 20.50831 POINT (-73.88702 10.28346)nmuestras<-na.omit(muestras)Ahora visualizamos las muestras que solicitamos anteriormente en el código
longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$`SOC_15-30`
id <-seq(1, 500, 1)(sitios <- data.frame(id, longit, latit, soc))## id longit latit soc
## 1 1 -74.25286 8.244205 28.0595741
## 2 2 -73.96963 8.286690 27.6881828
## 3 3 -75.41174 8.088429 35.9896393
## 4 4 -74.35435 10.243331 14.1246748
## 5 5 -75.67609 7.847685 16.0550575
## 6 6 -74.95621 8.211162 43.6211166
## 7 7 -74.93261 9.896376 14.4732351
## 8 8 -74.19149 7.559735 36.9720573
## 9 9 -73.99323 9.237868 17.0069294
## 10 10 -73.88702 10.283456 20.5083122
## 11 11 -74.05460 9.617867 11.9866705
## 12 12 -74.14901 9.405445 16.2152920
## 13 13 -73.80205 10.045071 14.2562990
## 14 14 -75.10255 9.235508 20.3817558
## 15 15 -75.69261 7.248183 81.5330963
## 16 16 -75.64540 8.881471 24.0316982
## 17 17 -74.53845 9.030166 21.2579098
## 18 18 -75.58404 8.501472 17.0446644
## 19 19 -74.49360 6.988556 35.1194115
## 20 20 -73.86106 9.768923 18.0306587
## 21 21 -73.79497 7.623461 25.5612011
## 22 22 -75.39286 8.425944 15.7945004
## 23 23 -75.47310 8.841347 17.4245758
## 24 24 -75.21584 7.411040 53.0650330
## 25 25 -75.73981 7.389797 55.6309433
## 26 26 -74.67062 9.778364 11.3915243
## 27 27 -74.81932 9.322837 25.9151917
## 28 28 -73.91770 7.408679 35.5900269
## 29 29 -75.46602 8.100231 15.8514414
## 30 30 -75.55335 9.480973 35.6574631
## 31 31 -74.97981 7.285947 50.5240097
## 32 32 -74.34019 10.340101 14.8848000
## 33 33 -73.94131 10.141841 14.5410833
## 34 34 -75.59112 9.728798 16.7827816
## 35 35 -75.23000 10.316499 15.0915537
## 36 36 -75.12143 10.613890 16.1540661
## 37 37 -73.92006 9.563581 15.3579721
## 38 38 -75.35745 9.683954 20.4302769
## 39 39 -75.04354 8.168678 13.9764652
## 40 40 -75.36453 7.564455 27.6942558
## 41 41 -75.61236 9.967183 0.0000000
## 42 42 -73.78081 8.281969 35.6615906
## 43 43 -75.37397 10.071034 58.5289421
## 44 44 -74.99634 8.359857 45.1824837
## 45 45 -75.72329 9.884574 0.0000000
## 46 46 -75.58640 9.532898 0.0000000
## 47 47 -73.99323 7.991659 21.2579823
## 48 48 -75.01522 7.566816 34.5338593
## 49 49 -75.58168 10.307058 0.0000000
## 50 50 -74.51249 10.481716 14.2073469
## 51 51 -73.83982 8.215882 31.5441551
## 52 52 -75.03882 10.741343 13.9842443
## 53 53 -74.88068 7.191537 44.7966690
## 54 54 -74.92789 8.383460 30.7113304
## 55 55 -74.37087 9.948301 15.0799723
## 56 56 -74.43224 9.570662 14.4590788
## 57 57 -75.77050 9.537619 0.0000000
## 58 58 -73.84926 10.585567 27.2293243
## 59 59 -74.39683 8.659608 37.8332825
## 60 60 -75.20404 9.698115 20.4407196
## 61 61 -74.56913 10.432151 13.0187483
## 62 62 -75.62180 7.960976 21.4078789
## 63 63 -74.64466 7.467685 56.4309044
## 64 64 -74.10652 8.496751 25.9415684
## 65 65 -75.28192 9.955382 16.2529278
## 66 66 -75.30317 7.179736 50.8672981
## 67 67 -75.62652 10.394387 0.0000000
## 68 68 -75.51087 8.317373 29.2433605
## 69 69 -74.09000 7.989299 21.9199753
## 70 70 -75.71149 8.992403 50.3206100
## 71 71 -74.76031 9.615507 13.6129875
## 72 72 -74.80280 7.205698 45.7384682
## 73 73 -74.46292 7.540853 54.3553696
## 74 74 -74.53609 7.460605 64.1545792
## 75 75 -75.21584 7.580977 18.4924908
## 76 76 -74.28354 9.433768 27.5855064
## 77 77 -75.22292 8.728055 1.8122197
## 78 78 -75.78230 7.968057 12.2852936
## 79 79 -74.68242 9.018365 0.0000000
## 80 80 -75.04590 9.164700 16.2944107
## 81 81 -74.38503 7.614020 58.6872520
## 82 82 -75.36925 7.451164 21.8577271
## 83 83 -74.65410 9.790165 12.9742393
## 84 84 -74.28826 8.886192 40.1886444
## 85 85 -75.66900 7.800480 19.0562515
## 86 86 -75.33385 9.672153 16.6132736
## 87 87 -75.46602 10.191406 34.3041763
## 88 88 -74.40391 8.803583 24.9862309
## 89 89 -75.24652 8.057746 16.7184219
## 90 90 -75.65012 8.772900 21.1964626
## 91 91 -75.19932 7.521971 37.3589439
## 92 92 -75.35509 7.715511 20.9047966
## 93 93 -75.42590 9.263830 14.2746305
## 94 94 -74.54317 7.665946 33.2875443
## 95 95 -73.91062 7.694269 23.3592606
## 96 96 -74.71783 10.658734 37.9015579
## 97 97 -74.06404 8.171038 34.4506950
## 98 98 -75.42354 10.441592 12.1390247
## 99 99 -73.99087 8.574639 19.0432796
## 100 100 -74.89721 8.529795 42.0425758
## 101 101 -74.85472 7.828803 22.0890369
## 102 102 -75.26540 9.938860 16.4884624
## 103 103 -73.78081 10.092276 16.0631237
## 104 104 -75.60764 8.959359 41.7959061
## 105 105 -74.26702 6.983836 25.5765305
## 106 106 -75.34329 7.722592 20.4725189
## 107 107 -75.48727 10.488797 11.3981791
## 108 108 -74.23870 10.731902 51.3881645
## 109 109 -74.72963 9.877494 15.1712742
## 110 110 -74.21509 10.550163 36.5672302
## 111 111 -74.15137 10.148922 13.9663429
## 112 112 -75.70205 8.827185 38.3547668
## 113 113 -74.40391 9.018365 26.3779716
## 114 114 -75.27248 7.609300 22.0922565
## 115 115 -74.21746 9.488054 15.7312441
## 116 116 -74.70839 8.706813 45.9137840
## 117 117 -74.33783 7.503089 48.9706230
## 118 118 -75.05062 10.311778 22.7992439
## 119 119 -74.02392 9.582463 11.8525381
## 120 120 -75.71621 9.469172 0.0000000
## 121 121 -74.03808 8.097870 32.2169838
## 122 122 -74.18205 8.005821 30.8005371
## 123 123 -74.25758 9.898736 15.7478533
## 124 124 -75.62888 7.694269 24.8907852
## 125 125 -74.72963 9.450290 29.1749477
## 126 126 -74.32367 9.030166 25.0922546
## 127 127 -74.86888 9.443209 13.8199930
## 128 128 -75.56987 7.590418 36.8620491
## 129 129 -74.40863 10.351903 14.5869856
## 130 130 -75.69261 8.175758 22.0624847
## 131 131 -75.39050 8.083709 0.0000000
## 132 132 -74.11832 9.247309 21.7568245
## 133 133 -74.73671 10.399107 27.7742882
## 134 134 -74.24814 6.981475 27.8280029
## 135 135 -75.12615 9.525817 16.4906502
## 136 136 -74.21982 7.840604 28.5663204
## 137 137 -75.74217 9.839730 0.0000000
## 138 138 -75.61472 9.544699 0.0000000
## 139 139 -74.28590 7.691909 23.2809753
## 140 140 -75.40466 8.810664 21.0996265
## 141 141 -75.48255 9.646190 30.0268497
## 142 142 -75.74453 7.977498 16.5385876
## 143 143 -74.03336 8.248926 26.1125393
## 144 144 -75.73981 9.901096 0.0000000
## 145 145 -73.85398 9.627308 12.3227663
## 146 146 -74.75795 7.368555 47.3747635
## 147 147 -75.76577 9.450290 0.0000000
## 148 148 -74.95621 9.270911 11.8435402
## 149 149 -73.96019 7.654145 14.9471560
## 150 150 -74.62814 9.745320 12.8269272
## 151 151 -74.22926 8.466068 31.9356022
## 152 152 -75.28192 9.974264 16.2202759
## 153 153 -74.80044 8.114392 27.7133732
## 154 154 -74.42280 7.344953 32.1170959
## 155 155 -74.45348 7.257624 34.6883888
## 156 156 -75.37633 9.407805 25.4416828
## 157 157 -75.20876 7.090046 48.3874321
## 158 158 -73.81149 8.520354 30.6227131
## 159 159 -74.89012 7.876007 20.0861053
## 160 160 -74.05460 7.892529 26.3588333
## 161 161 -74.94441 9.353520 21.6923180
## 162 162 -75.26776 10.078114 24.3850307
## 163 163 -73.84454 8.348056 17.5060062
## 164 164 -74.53373 7.047562 30.1888561
## 165 165 -74.78391 9.960102 5.1063652
## 166 166 -74.61162 10.179605 12.8269796
## 167 167 -73.92478 7.519611 0.0000000
## 168 168 -75.56279 8.966440 39.1245155
## 169 169 -75.15919 8.940477 27.0464516
## 170 170 -74.27174 9.804326 14.2341280
## 171 171 -73.92951 10.717740 27.7103500
## 172 172 -75.14267 10.474635 9.7211094
## 173 173 -74.08764 8.364578 27.1710606
## 174 174 -74.97273 8.305572 2.9756329
## 175 175 -74.49360 9.945941 15.0044699
## 176 176 -74.89957 9.804326 20.4382572
## 177 177 -74.99161 10.720101 10.4670267
## 178 178 -74.92553 10.786187 11.2140322
## 179 179 -75.07186 7.205698 32.9538422
## 180 180 -75.25360 8.935757 16.5039997
## 181 181 -74.69658 8.923955 0.0000000
## 182 182 -74.51013 10.616250 0.1959799
## 183 183 -74.55497 9.811407 13.1749306
## 184 184 -74.66826 8.879111 39.0596962
## 185 185 -75.20876 9.624948 18.3471413
## 186 186 -74.84528 9.891655 18.3610134
## 187 187 -74.79335 10.675256 12.4976206
## 188 188 -74.35907 10.382586 11.9426918
## 189 189 -73.99795 10.401468 17.8226166
## 190 190 -74.28826 7.335512 50.2104836
## 191 191 -74.34019 7.347313 46.7970085
## 192 192 -75.70913 8.015262 17.7530746
## 193 193 -75.17571 7.578617 23.3168869
## 194 194 -74.18205 7.755635 29.4742527
## 195 195 -75.53447 9.627308 27.2134552
## 196 196 -75.24416 8.447186 32.8687019
## 197 197 -74.92317 10.094636 21.6410675
## 198 198 -74.55969 8.392901 34.4433022
## 199 199 -75.34801 7.106568 56.9591599
## 200 200 -74.40627 8.789422 0.0000000
## 201 201 -74.70367 9.608426 0.0000000
## 202 202 -73.78317 10.663455 66.7288513
## 203 203 -75.45422 10.441592 42.4318657
## 204 204 -74.88776 9.273271 24.1906681
## 205 205 -74.97745 9.516377 14.5041237
## 206 206 -73.86814 7.774517 23.5647144
## 207 207 -75.62416 7.075885 78.3653336
## 208 208 -75.69261 8.879111 31.8892422
## 209 209 -74.50068 9.023086 54.9295578
## 210 210 -74.68714 8.036504 33.1823540
## 211 211 -74.30950 7.264704 36.8051262
## 212 212 -74.54789 8.928676 58.5832901
## 213 213 -74.97037 8.598242 38.5954781
## 214 214 -74.16789 7.262344 57.2308807
## 215 215 -75.14503 9.698115 14.1376953
## 216 216 -75.42354 7.458244 31.6339283
## 217 217 -74.17497 8.199361 27.8660202
## 218 218 -75.11907 8.093150 18.9753475
## 219 219 -74.54789 7.177375 42.5224304
## 220 220 -73.88938 10.337741 22.3210144
## 221 221 -75.74689 10.481716 0.0000000
## 222 222 -74.30950 8.357497 32.6727905
## 223 223 -74.25994 10.127679 15.5015678
## 224 224 -75.15447 7.573896 25.7704239
## 225 225 -75.46130 8.999483 28.7146702
## 226 226 -74.59037 10.328300 15.3485432
## 227 227 -73.86814 7.939734 31.7563953
## 228 228 -75.06714 10.160723 13.2438822
## 229 229 -74.55261 7.668306 57.7816505
## 230 230 -74.30950 10.172524 19.1536789
## 231 231 -75.77050 9.438488 0.0000000
## 232 232 -75.05770 8.782341 46.3950424
## 233 233 -73.79497 9.438488 23.5303459
## 234 234 -73.98615 10.455753 18.1149254
## 235 235 -75.06714 7.283586 37.9295769
## 236 236 -74.77447 7.167935 44.1526260
## 237 237 -75.19932 9.681594 18.9786015
## 238 238 -74.32367 9.419607 16.1867275
## 239 239 -75.39994 10.616250 20.5357113
## 240 240 -75.01522 7.821722 13.8569450
## 241 241 -73.96491 7.363835 21.2167988
## 242 242 -74.17025 10.687057 24.5564690
## 243 243 -75.70441 8.527435 17.2506218
## 244 244 -74.55025 7.762716 54.5596886
## 245 245 -73.87286 9.270911 12.0743866
## 246 246 -73.85634 7.845324 20.6918468
## 247 247 -75.34329 7.052283 38.7425995
## 248 248 -74.94205 8.232404 27.5439510
## 249 249 -74.26466 9.157620 21.8309402
## 250 250 -74.27410 10.644573 26.2515049
## 251 251 -75.78702 9.832649 0.0000000
## 252 252 -75.50143 8.822465 16.7415981
## 253 253 -74.09944 9.912897 15.7468557
## 254 254 -75.27012 8.289050 26.3729992
## 255 255 -75.10019 7.956256 14.5330791
## 256 256 -75.13087 8.525074 35.5236626
## 257 257 -73.82329 8.300851 18.9715214
## 258 258 -74.52429 7.106568 32.3540993
## 259 259 -74.86180 9.830289 28.7962627
## 260 260 -75.69497 9.273271 29.2306957
## 261 261 -75.69025 10.392027 0.0000000
## 262 262 -73.92478 8.763459 37.0896378
## 263 263 -75.63596 7.580977 52.3955765
## 264 264 -73.92242 7.821722 14.6112347
## 265 265 -74.07112 10.243331 12.7419481
## 266 266 -73.86578 8.921595 42.1305695
## 267 267 -75.19932 7.340232 43.4576187
## 268 268 -74.69658 8.067187 30.4639816
## 269 269 -74.04988 10.774386 28.2198029
## 270 270 -75.45658 10.757864 0.0000000
## 271 271 -75.60764 8.843707 19.7743473
## 272 272 -74.13013 8.326814 29.8043385
## 273 273 -75.46602 10.448673 23.1045895
## 274 274 -75.14031 8.742217 17.5192661
## 275 275 -74.82640 10.620970 8.8118439
## 276 276 -73.94839 7.996380 49.1219864
## 277 277 -74.26938 10.455753 18.1505604
## 278 278 -74.19385 10.788548 42.6184616
## 279 279 -74.96093 10.271654 21.7528419
## 280 280 -74.41100 8.475509 33.1210480
## 281 281 -73.80205 9.018365 0.0000000
## 282 282 -75.71385 9.506936 0.0000000
## 283 283 -74.35907 9.308675 24.7012882
## 284 284 -74.23634 8.324454 29.6211014
## 285 285 -74.34019 9.752401 17.8028507
## 286 286 -75.45186 7.153773 59.4960556
## 287 287 -75.75397 10.146561 0.0000000
## 288 288 -73.98143 9.561221 13.9181585
## 289 289 -74.46528 7.106568 53.0538254
## 290 290 -74.51013 9.577743 14.2614975
## 291 291 -74.51721 9.169421 62.3546028
## 292 292 -74.81932 9.403085 16.3266296
## 293 293 -73.88230 10.295257 21.8809929
## 294 294 -74.93733 7.111289 45.7378273
## 295 295 -75.10491 7.668306 21.8612080
## 296 296 -75.75397 9.211905 35.3634262
## 297 297 -74.30478 7.057003 35.2203140
## 298 298 -75.30317 9.023086 27.5038509
## 299 299 -75.30081 10.118239 0.0000000
## 300 300 -73.89410 8.168678 32.5220985
## 301 301 -75.31969 10.670535 16.5048618
## 302 302 -74.66590 8.248926 64.1319962
## 303 303 -75.39758 9.462091 13.1375198
## 304 304 -74.17733 8.048305 29.0594521
## 305 305 -73.89646 8.461348 25.4013367
## 306 306 -75.76341 9.159980 55.1475868
## 307 307 -75.58404 7.694269 24.1824341
## 308 308 -73.90354 10.514759 25.1021252
## 309 309 -73.79969 6.993277 21.8346920
## 310 310 -75.01758 8.940477 19.5638237
## 311 311 -74.86180 8.433025 31.5268192
## 312 312 -74.65646 9.221346 0.0000000
## 313 313 -75.70205 8.687931 26.3330727
## 314 314 -74.85708 10.042711 17.9209137
## 315 315 -75.31969 9.988425 46.4407578
## 316 316 -74.70367 8.463708 33.7927513
## 317 317 -75.68317 9.200104 91.2178116
## 318 318 -74.58565 8.199361 29.0698757
## 319 319 -75.46366 8.451907 16.3342648
## 320 320 -75.64304 10.644573 0.0000000
## 321 321 -75.37633 7.453524 21.7673130
## 322 322 -74.86652 7.642343 17.1898956
## 323 323 -75.67845 8.237125 21.7445793
## 324 324 -75.49435 10.184325 28.1117001
## 325 325 -74.41572 9.412526 15.8163729
## 326 326 -75.71385 9.327557 17.0138741
## 327 327 -74.54317 8.291410 37.6947975
## 328 328 -75.04354 7.507810 28.6770802
## 329 329 -74.64938 7.321350 65.1273193
## 330 330 -75.17807 7.698989 16.6413002
## 331 331 -74.16081 10.370785 14.1314392
## 332 332 -74.44640 7.337872 29.0166893
## 333 333 -74.88304 9.535258 16.0938606
## 334 334 -74.79335 9.473892 0.0000000
## 335 335 -73.96019 9.830289 14.9160128
## 336 336 -74.11360 7.333152 34.0655174
## 337 337 -73.78081 9.502215 22.5179176
## 338 338 -74.23162 7.595139 49.5554733
## 339 339 -75.06478 10.557244 15.3317900
## 340 340 -75.73745 9.768923 0.0000000
## 341 341 -74.88304 10.774386 9.3639059
## 342 342 -73.86106 7.524331 22.6585674
## 343 343 -75.53211 10.776746 0.0000000
## 344 344 -74.35435 9.348799 19.3950233
## 345 345 -73.92715 9.657991 12.4283953
## 346 346 -74.74615 9.799606 13.9366722
## 347 347 -75.44714 9.254390 13.8457527
## 348 348 -74.27174 9.200104 23.7888298
## 349 349 -75.23708 8.310292 29.3501072
## 350 350 -73.94131 9.190663 15.1764717
## 351 351 -74.59745 9.610786 20.2457848
## 352 352 -74.12069 7.193897 39.4994392
## 353 353 -74.61398 7.991659 40.3494148
## 354 354 -74.65174 8.078988 24.5667019
## 355 355 -74.12069 10.458114 14.9580345
## 356 356 -75.18515 7.793399 19.6461868
## 357 357 -75.79646 7.781598 31.5877934
## 358 358 -74.34491 9.988425 16.7148781
## 359 359 -75.50379 10.108798 50.2695160
## 360 360 -74.97037 8.796502 52.0738029
## 361 361 -74.23870 10.229170 14.5109568
## 362 362 -73.92951 7.245823 21.3454742
## 363 363 -75.37161 9.006564 26.9441109
## 364 364 -75.04826 7.772157 14.6193638
## 365 365 -74.19857 10.524200 0.0000000
## 366 366 -75.23472 9.252029 26.6333828
## 367 367 -73.79969 9.334638 11.9939899
## 368 368 -74.84528 10.198487 17.7445107
## 369 369 -75.33621 10.052152 86.1205444
## 370 370 -74.19857 7.175015 77.3114471
## 371 371 -74.38975 8.369298 21.5347557
## 372 372 -74.70130 9.308675 12.9126720
## 373 373 -73.81621 9.551780 13.5404100
## 374 374 -75.18043 6.993277 49.5176315
## 375 375 -74.76031 9.223706 7.8218961
## 376 376 -75.77050 9.259110 3.5032337
## 377 377 -74.13721 8.338615 36.1331635
## 378 378 -75.69497 9.433768 0.0000000
## 379 379 -74.71783 9.429048 25.5360966
## 380 380 -75.08602 10.264574 19.7083492
## 381 381 -73.78553 8.895633 12.0929241
## 382 382 -75.40230 8.754018 25.3357887
## 383 383 -75.65720 10.446312 0.0000000
## 384 384 -74.99161 8.409422 30.2557545
## 385 385 -74.25758 7.493648 26.7725277
## 386 386 -74.69186 7.970417 21.5228310
## 387 387 -75.45422 10.654014 0.0000000
## 388 388 -75.12615 9.997866 15.7536669
## 389 389 -75.68317 9.776003 0.0000000
## 390 390 -74.20093 8.999483 25.6987801
## 391 391 -75.03174 7.696629 26.2509804
## 392 392 -74.81932 7.555014 14.1908169
## 393 393 -73.85162 9.063210 0.0000000
## 394 394 -73.77845 8.775260 36.7131462
## 395 395 -75.63124 8.001100 22.1033897
## 396 396 -74.01447 10.595008 22.9290886
## 397 397 -74.00739 9.322837 13.3228683
## 398 398 -74.44404 7.583337 37.7722855
## 399 399 -75.56043 10.684697 0.0000000
## 400 400 -73.92006 8.274889 22.1767673
## 401 401 -75.23944 7.786318 22.9160423
## 402 402 -75.59584 7.566816 29.1019268
## 403 403 -74.56913 10.717740 0.0000000
## 404 404 -74.39683 8.024703 46.7152901
## 405 405 -74.70839 7.101848 42.7139435
## 406 406 -75.03882 10.089916 12.1551085
## 407 407 -74.55733 7.111289 47.8646240
## 408 408 -75.25596 10.375505 12.5350046
## 409 409 -74.28118 9.233147 0.4303940
## 410 410 -74.10652 7.847685 20.3173199
## 411 411 -75.01286 9.023086 19.1815758
## 412 412 -73.77373 10.037990 16.1455536
## 413 413 -74.08056 10.297617 14.2290277
## 414 414 -74.89012 10.455753 25.2362232
## 415 415 -73.84690 10.158363 20.4094696
## 416 416 -75.15683 7.819362 26.0349236
## 417 417 -75.71621 7.977498 14.9660892
## 418 418 -74.17261 10.406188 12.6749897
## 419 419 -74.28826 7.724952 38.8337288
## 420 420 -74.57857 10.262213 12.8816414
## 421 421 -75.04354 10.153642 10.2246799
## 422 422 -73.99795 8.772900 17.8991261
## 423 423 -75.06242 7.760356 18.0177612
## 424 424 -74.44404 8.897993 7.7517271
## 425 425 -73.78789 10.219729 17.7356777
## 426 426 -74.95621 10.019108 12.2328329
## 427 427 -74.10652 8.376379 31.7338486
## 428 428 -75.73037 9.502215 0.0000000
## 429 429 -74.80280 8.810664 42.8685036
## 430 430 -75.58404 10.498238 0.0000000
## 431 431 -75.44478 9.818488 18.7194519
## 432 432 -75.55571 9.990785 61.0913696
## 433 433 -74.19385 9.978984 13.4338846
## 434 434 -74.32839 10.538362 11.3845530
## 435 435 -75.23472 10.533641 17.8251820
## 436 436 -74.70839 10.231530 23.5194016
## 437 437 -75.69261 9.889295 0.0000000
## 438 438 -75.48019 8.404702 16.7089748
## 439 439 -75.53211 7.647064 21.3341084
## 440 440 -74.81932 10.592647 10.1963873
## 441 441 -74.39211 8.626565 46.3168373
## 442 442 -74.27410 7.614020 53.7703972
## 443 443 -75.50851 10.170164 39.0115700
## 444 444 -73.94603 9.207185 18.7911339
## 445 445 -74.41808 10.691778 65.3523331
## 446 446 -73.93895 9.504575 14.3253098
## 447 447 -74.05460 7.035761 21.4302788
## 448 448 -75.70677 9.414886 27.9034481
## 449 449 -74.13721 9.077371 15.0035095
## 450 450 -73.78317 10.307058 31.5900517
## 451 451 -75.27012 8.050665 12.9159327
## 452 452 -75.10963 7.953896 17.1706924
## 453 453 -74.52429 9.110415 0.0000000
## 454 454 -75.12615 9.660351 15.7549944
## 455 455 -75.74689 9.384203 23.5570164
## 456 456 -74.92789 7.302468 37.3928261
## 457 457 -74.38267 8.270168 40.7100105
## 458 458 -75.59348 7.800480 19.6484795
## 459 459 -74.22454 10.134760 14.3765240
## 460 460 -75.73273 8.536876 15.8556213
## 461 461 -73.83510 10.189046 19.1649475
## 462 462 -74.79571 8.775260 32.1847649
## 463 463 -75.17099 8.905074 30.4974594
## 464 464 -74.54317 10.616250 49.2132263
## 465 465 -74.63286 10.715380 35.4047623
## 466 466 -75.39522 9.908177 18.7799129
## 467 467 -74.96329 7.911411 19.2602749
## 468 468 -75.39522 7.793399 20.4142609
## 469 469 -74.79099 10.167804 14.0400410
## 470 470 -74.06640 9.037247 20.8697414
## 471 471 -75.35273 7.559735 36.4474106
## 472 472 -75.42590 7.524331 24.8207378
## 473 473 -74.94441 7.474766 35.1183205
## 474 474 -75.23472 8.359857 27.8221912
## 475 475 -74.12777 8.015262 29.6537952
## 476 476 -75.52267 7.769797 20.9746494
## 477 477 -74.72727 7.708430 28.4315701
## 478 478 -74.65646 8.643087 30.2252007
## 479 479 -74.45820 7.061724 24.5610142
## 480 480 -75.55807 8.520354 17.5375290
## 481 481 -74.33075 10.590287 16.8964329
## 482 482 -74.08056 7.878368 28.3671131
## 483 483 -75.44950 8.631285 18.5950184
## 484 484 -74.68714 8.227684 51.3112030
## 485 485 -74.34255 9.662712 19.9625359
## 486 486 -75.60764 8.163957 28.6952820
## 487 487 -75.74925 8.463708 16.6750546
## 488 488 -73.90826 10.448673 21.0236073
## 489 489 -75.54863 8.055386 12.3588257
## 490 490 -73.77373 9.455010 24.0866623
## 491 491 -74.61398 10.550163 49.5359840
## 492 492 -74.67534 7.715511 42.3515358
## 493 493 -74.41572 7.106568 40.9320984
## 494 494 -74.25758 9.657991 15.8988886
## 495 495 -75.27012 7.208059 71.0203476
## 496 496 -73.95311 7.269425 24.1373806
## 497 497 -75.25124 7.267065 59.4436722
## 498 498 -75.02938 8.128553 13.0947580
## 499 499 -75.10019 7.750915 19.0530128
## 500 500 -74.38503 7.946815 39.3919868Con la finalidad de limpiar la tabla se corre un código para eliminar datos nulos
sitios<-na.omit(sitios)
head (sitios)## id longit latit soc
## 1 1 -74.25286 8.244205 28.05957
## 2 2 -73.96963 8.286690 27.68818
## 3 3 -75.41174 8.088429 35.98964
## 4 4 -74.35435 10.243331 14.12467
## 5 5 -75.67609 7.847685 16.05506
## 6 6 -74.95621 8.211162 43.62112Ahora nuevamente se vuelve a dar la istrucción para visualizar las muestras
m <- leaflet() %>%
addTiles() %>%
leafem:::addGeoRaster(
stars.soc,
opacity = 0.7,
colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
domain = 8:130)
) %>%
addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())Esta información es la que nos va a ayudar a hacer el ejercico de interpolación
4. Interpolación
Hay que cumplir con ciertos parametros para poder interpolar, inicialmente se debe crear un objeto gstat ¿Por qué? contiene toda la información necesaria para realizar la interpolación que es La definición del modelo y los datos de calibración ya sea para la interpolación OK o IDW ### 4.1 Interpolación IDW Para poder interpolar nuestro SOC por IDW, inicialmente se debe crear el objeto gstat
g1 = gstat(formula = soc ~ 1, data = nmuestras)Teniendo lo anterior es posible para hacer una estimación de los valores, dicha función opta por un raster o un modelo gstat Vamos a crear un objeto ráster con valores de celda iguales a 1:
rrr = aggregate(geog.soc, 4)rrr## class : SpatRaster
## dimensions : 405, 215, 1 (nrow, ncol, nlyr)
## resolution : 0.009440973, 0.009440973 (x, y)
## extent : -75.79764, -73.76783, 6.977935, 10.80153 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : SOC_15-30
## min value : 0.0000
## max value : 119.2192Se definen nuevos valores
values(rrr) <-1Definir nuevos nombres:
names(rrr) <- "valor"Reproyectamos rrr
rrr## class : SpatRaster
## dimensions : 405, 215, 1 (nrow, ncol, nlyr)
## resolution : 0.009440973, 0.009440973 (x, y)
## extent : -75.79764, -73.76783, 6.977935, 10.80153 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : valor
## min value : 1
## max value : 1stars.rrr = st_as_stars(rrr)la siguiente expresión interpola los valores SOC según el modelo definido y stars.rrr
## [inverse distance weighted interpolation]
z1 = predict(g1, stars.rrr)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/Hp/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/Hp/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.## [inverse distance weighted interpolation]z1## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## var1.pred 0.01840842 18.09741 23.2624 25.05296 30.37548 90.07681 0
## var1.var NA NA NA NaN NA NA 87075
## dimension(s):
## from to offset delta refsys x/y
## x 1 215 -75.7976 0.00944097 +proj=longlat +datum=WGS8... [x]
## y 1 405 10.8015 -0.00944097 +proj=longlat +datum=WGS8... [y]A continuación se creará un subconjunto solo del primer atributo y se le cambiar el nombre a “soc”
z1 = z1["var1.pred",,]
names(z1) = "soc"Se configura la paleta de colores
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")Finalmente el raster SOC interpolado por el metodo IDW será posible de visualizar a continuación teniendo en cuenta los ajustes que se hicieron tanto técnicos como físicos para facilitar su interpretación
m <- leaflet() %>%
addTiles() %>%
leafem:::addGeoRaster(
z1,
opacity = 0.7,
colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
domain = 11:55)
) %>%
addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
addLegend("bottomright", pal=paleta, values= z1$soc,
title = "IDW SOC interpolation [%]"
)## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NAInterpolación OK
El modelo kriging ordinary hace uso de la herramienta del r variograma el cual establece una autorrlación en los datos y asignar pesos en consecuencia al hacer predicciones correspondientes a la interpolación, este modelo sigue unos pasos, iniciando por calcular y examinar el variograma empírico, para esto se necesita una formula que debe cumplir con lo siguiente fórmula: especifica la variable dependiente y las covariables, al igual que en gstat datos: la capa de puntos con la variable dependiente y las covariables como atributos de puntos
v_emp_ok = variogram(soc ~ 1, data=nmuestras)Con el código anterior ya es posible dibujar el variograma
plot(v_emp_ok)
Como nuestro objetivo es ajustar un modelo de variograma a un variograma
empírico, hay dos maneras de hacerlo, pero para el presente caso se usa
el metodo mas simple que es ajuste automático usando la función **
autofitVariogram ** del paquete ** automap **
v_mod_ok = autofitVariogram(soc ~ 1, as(nmuestras, "Spatial"))## Warning in autofitVariogram(soc ~ 1, as(nmuestras, "Spatial")): Some models where removed for being either NULL or having a negative sill/range/nugget,
## set verbose == TRUE for more informationComo se mencionó anteriormente, la función elige el modelo que mejor se ajusta. Cabe resaltar que tuvimos que convertir las muestras a ** SpatialPintsDataFrame ** puesto que la función usada no es compatible con objetos sf
Ahora si, revisamos el modelo obtenido:
plot(v_mod_ok)
El objeto resultante de hecho es una lista con muchos componentes,
incluyendo el variograma empírico y ajustado. El componente $var_model
del objeto resultante, contiene al modelo actual, lo podemos ver con el
siguiente código
v_mod_ok$var_model## model psill range kappa
## 1 Nug 0.741551 0.00000 0.0
## 2 Ste 203.684570 29.80368 0.2Ahora, el modelo de variograma puede ser utilizado con la función gstat, y podemos continuar con la interpolación de Kriging ordinario
## [using ordinary kriging]
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)## [using ordinary kriging]Nuevamente, se subagrupan los valores predichos y se les renombra:
z2 = z2["var1.pred",,]
names(z2) = "soc"Ahora observamos el siguiente mapa con los resultados:
m <- leaflet() %>%
addTiles() %>%
leafem:::addGeoRaster(
z2,
opacity = 0.7,
colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
domain = 11:55)
) %>%
addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
addLegend("bottomright", pal = paleta, values= z2$soc,
title = "OK SOC interpolation [%]"
)## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA5. Evaluación de los datos obtenidos
Evaluación cualitativa
Otra vista de los resultados de la interpolación.
Primero procuramos la paleta de colores
colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")m <- leaflet() %>%
addTiles() %>%
addGeoRaster(stars.soc, opacity = 0.8, colorOptions = colores, group="RealWorld") %>%
addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW") %>%
addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK") %>%
# Add layers controls
addLayersControl(
overlayGroups = c("RealWorld", "IDW", "OK"),
options = layersControlOptions(collapsed = FALSE)
) %>%
addLegend("bottomright", pal = paleta, values= z1$soc,
title = "Soil organic carbon [%]"
)## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA