1. introducción

Este cuaderno ilustra dos técnicas de interpolación espacial: distancia ponderada inversa (IDW) y Kriging ordinario (OK). IDW es una técnica determinista. OK es probabilístico. Ambas técnicas se utilizan aquí para obtener una superficie continua de SOC a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m

2. Configuración

Primero es necesario limpiar la memoria

rm(list=ls())

Luego, asegúrese de haber instalado previamente las bibliotecas necesarias.

#install.packages("sp")
library(sp)
#install.packages("terra")
library(terra)
## terra 1.7.55
#install.packages("sf")
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
#install.packages("stars")
library(stars)
## Loading required package: abind
#install.packages("leaflet")
library(leaflet)
#install.packages("dplyr")
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, union
#install.packages("gstat")
library(gstat)
#install.packages("automap")
library(automap)
#install.packages("leafem")
library(leafem)
#install.packages("curl")
library(curl)
## Using libcurl 8.3.0 with Schannel
h <- new_handle()
handle_setopt(h, http_version = 2)

#3. Leer los datos de entrada

Necesitamos leer un conjunto de datos para imitar datos del mundo real. Por lo tanto, leamos la capa SOC que descargamos de ISRIC) usando la biblioteca terra:

archivo <- ("soc_igh_15_30.tif")
(soc <- rast(archivo))
## class       : SpatRaster 
## dimensions  : 1068, 879, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8313500, -8093750, 637000, 904000  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

Ahora, convierta los datos del SOC en porcentaje. Revise el factor de escala de SOC en el sitio web de SoilGrids y anótelo aquí.

soc.perc <-  soc/10

¿Cuál es el CRS de los datos del mundo real?

El Sistema de Referencia de Coordenadas (CRS) es un sistema de coordenadas que se utiliza para describir la ubicación de un objeto en la Tierra. El sistema de referencia de coordenadas que se muestra en la capa SOC que descargamos de ISRIC es el Interrupted Goode Homolosine. Este sistema de referencia de coordenadas se utiliza para representar la superficie terrestre en un mapa y se basa en una proyección interrumpida que minimiza las distorsiones en la forma y el tamaño de los continentes.

El sistema de referencia de coordenadas más utilizado en todo el mundo es el Sistema de Coordenadas Geográficas (WGS84). Este sistema de referencia de coordenadas se utiliza para la navegación, la cartografía y la geodesia, y es el sistema de referencia de coordenadas utilizado por el Sistema de Posicionamiento Global (GPS)

Vamos a una transformación de CRS Interrupted Goode Homolosineal conocido CRS WGS84:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## dimensions  : 1094, 969, 1  (nrow, ncol, nlyr)
## resolution  : 0.002191783, 0.002191783  (x, y)
## extent      : -74.55464, -72.43081, 5.72296, 8.12077  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      9.687282 
## max value   :    128.557999

Convirtamos la capa SpatRaster en un objeto de estrellas:

stars.soc = st_as_stars(geog.soc)
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("yellow", "orange", "green", "blue"), 
                                  domain = 8:130)
    ) 
#
m  #Imprime el mapa

#4. Muestreo del mundo

Consigamos una muestra de aproximadamente 500 sitios a partir de datos del mundo real utilizando una muestra ubicada aleatoriamente:

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      : -74.55355, -72.44067, 5.730631, 8.115291  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         42.45
##                        45.16
##                        19.95

Las caracteristicas del objeto de muestra es de clase SpatVector para representar datos geoespaciales con geometría de puntos. El objeto tiene 500 geometrías y 1 atributo. La extensión geografica del objeto es de -74.55355 a -72.44067 en el eje x y de 5.730631 a 8.115291 en el eje y. La referencia de coordenadas del objeto es +proj=longlat +datum=WGS84 +no_defs.

Vamos a convertir el objeto spatVector en un objeto de característica simple:

(muestras <- sf::st_as_sf(samples))
## Simple feature collection with 500 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.55355 ymin: 5.730631 xmax: -72.44067 ymax: 8.115291
## 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_igh_15_30                   geometry
## 1       42.45292 POINT (-73.85437 5.745973)
## 2       45.15730 POINT (-73.48834 5.785425)
## 3       19.95003  POINT (-73.9793 6.563508)
## 4       41.47083 POINT (-73.03464 7.722962)
## 5       25.04777 POINT (-73.29985 6.495563)
## 6       12.87853 POINT (-73.70533 8.115291)
## 7       17.27739 POINT (-73.62204 7.280221)
## 8       33.53328  POINT (-72.84834 7.59803)
## 9       53.53739 POINT (-72.99081 7.639674)
## 10      54.60721 POINT (-72.68177 7.021591)
nmuestras <- na.omit(muestras)

Visualicemos las muestras:

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$soc_igh_15_30
id <- seq(1,500,1) 
(sitios <- data.frame(id, longit, latit, soc))
##      id    longit    latit      soc
## 1     1 -73.85437 5.745973 42.45292
## 2     2 -73.48834 5.785425 45.15730
## 3     3 -73.97930 6.563508 19.95003
## 4     4 -73.03464 7.722962 41.47083
## 5     5 -73.29985 6.495563 25.04777
## 6     6 -73.70533 8.115291 12.87853
## 7     7 -73.62204 7.280221 17.27739
## 8     8 -72.84834 7.598030 33.53328
## 9     9 -72.99081 7.639674 53.53739
## 10   10 -72.68177 7.021591 54.60721
## 11   11 -73.10040 6.824331 27.56308
## 12   12 -73.87410 7.418304 23.19725
## 13   13 -74.49875 6.666522 38.91019
## 14   14 -73.70533 6.337755 23.58551
## 15   15 -72.75629 6.475837 35.27520
## 16   16 -74.30807 5.984878 32.49445
## 17   17 -73.21437 6.447344 29.95158
## 18   18 -72.81766 6.223782 44.66993
## 19   19 -73.57820 6.300494 24.64726
## 20   20 -74.22916 7.170632 56.47837
## 21   21 -73.13766 6.747618 24.51562
## 22   22 -72.50423 6.894468       NA
## 23   23 -74.03410 7.692277 22.91378
## 24   24 -72.55903 7.508167 38.13704
## 25   25 -72.83081 7.124605 45.39619
## 26   26 -72.56341 7.670359 26.22174
## 27   27 -73.44012 7.946523 22.98853
## 28   28 -72.99738 6.971180 42.55526
## 29   29 -73.45547 7.082961 19.14106
## 30   30 -73.49053 7.933373 14.40394
## 31   31 -72.99738 6.973372 46.98222
## 32   32 -73.98588 5.781042 39.35436
## 33   33 -74.55355 7.442413       NA
## 34   34 -73.20560 5.853371 41.60892
## 35   35 -74.28177 7.269263 34.57869
## 36   36 -73.78642 6.942687 22.82866
## 37   37 -73.10697 7.072002       NA
## 38   38 -72.59190 7.823784 18.14644
## 39   39 -72.49985 8.064880 17.13908
## 40   40 -72.81985 5.875289 35.03060
## 41   41 -74.15903 7.328441 32.23650
## 42   42 -73.85875 6.977755 34.46631
## 43   43 -73.06971 6.947070 26.34098
## 44   44 -74.03190 7.920222 25.99850
## 45   45 -72.57218 6.131727 58.60030
## 46   46 -72.69930 7.096112 51.28063
## 47   47 -73.67026 7.777756 26.45388
## 48   48 -74.12615 6.793646 32.53820
## 49   49 -73.06971 7.911455 36.03264
## 50   50 -73.24725 6.462686 28.34337
## 51   51 -73.53656 7.335016 21.07949
## 52   52 -72.71464 7.505975 50.01265
## 53   53 -73.98149 6.228166 24.30491
## 54   54 -73.59136 6.357481 31.34481
## 55   55 -74.52944 6.440769       NA
## 56   56 -73.01272 7.019399 49.16971
## 57   57 -74.45930 6.589810 45.42380
## 58   58 -73.55190 7.896112 36.29366
## 59   59 -73.89382 6.850632 37.21909
## 60   60 -72.68615 6.195289 55.30422
## 61   61 -74.18752 6.464878 28.02065
## 62   62 -72.62697 6.600769 37.94643
## 63   63 -74.42423 6.521865 20.18045
## 64   64 -73.80615 6.405700 31.53595
## 65   65 -73.12451 7.181591 27.14919
## 66   66 -72.45601 6.342138       NA
## 67   67 -73.52779 7.072002 27.63693
## 68   68 -72.60944 7.554194 37.00108
## 69   69 -74.23574 6.265426 31.74939
## 70   70 -73.93108 6.236933 21.74888
## 71   71 -74.44834 6.692824 48.74871
## 72   72 -73.07190 6.309262 38.17021
## 73   73 -73.84341 6.791454 30.17904
## 74   74 -73.64615 7.786523 21.35337
## 75   75 -72.74752 6.052823 37.15347
## 76   76 -72.67957 6.011179 32.11831
## 77   77 -74.40012 7.319674 39.05693
## 78   78 -72.89875 7.462139 59.01918
## 79   79 -72.89218 6.410084 44.95979
## 80   80 -73.29766 7.961866 46.90643
## 81   81 -72.64670 5.785425 71.51294
## 82   82 -74.33437 8.010085 38.54310
## 83   83 -72.81108 7.887345 29.88874
## 84   84 -74.03629 7.514742 24.04203
## 85   85 -74.36725 6.287344 28.80410
## 86   86 -73.41382 6.464878 45.45023
## 87   87 -73.25382 6.901043 23.78895
## 88   88 -73.54314 6.175563 52.12928
## 89   89 -73.84999 7.665975 35.75359
## 90   90 -72.60067 6.988714 45.05716
## 91   91 -73.40286 6.883509 26.25704
## 92   92 -74.05163 7.405153 22.85619
## 93   93 -73.14642 7.468715 33.51587
## 94   94 -72.98423 6.865974 30.21781
## 95   95 -73.82588 6.475837 45.35667
## 96   96 -72.67519 6.859399 35.13875
## 97   97 -72.96012 7.703236 72.88007
## 98   98 -73.71190 6.677481 33.39132
## 99   99 -72.74314 7.747071 21.51763
## 100 100 -73.95958 6.486796 28.31949
## 101 101 -74.47903 6.936112 45.85521
## 102 102 -73.04560 7.227619 27.05100
## 103 103 -73.15300 6.953646 22.44392
## 104 104 -73.42259 6.272001 31.16130
## 105 105 -73.85437 7.047892 30.17341
## 106 106 -72.88341 5.750357 40.39511
## 107 107 -73.87190 7.284605 20.66852
## 108 108 -72.58971 7.030358 71.94563
## 109 109 -73.29766 6.865974 40.64669
## 110 110 -72.68834 6.699399 49.83183
## 111 111 -73.77327 7.139947 20.75211
## 112 112 -72.98423 5.952001 44.18787
## 113 113 -74.49218 7.352550       NA
## 114 114 -73.63081 6.410084 23.48857
## 115 115 -74.35629 7.915838 36.77057
## 116 116 -74.22697 6.826522 27.92758
## 117 117 -74.32999 6.002412 32.09584
## 118 118 -73.06533 7.172824 31.93351
## 119 119 -73.19026 7.448989 21.75216
## 120 120 -74.03190 5.842412 27.12596
## 121 121 -73.54314 7.339400 23.38866
## 122 122 -74.45711 7.543235 51.40715
## 123 123 -74.25985 7.205701 92.86704
## 124 124 -74.30807 6.392549 25.98410
## 125 125 -73.64615 7.194742 24.11949
## 126 126 -72.66204 8.042962 20.93973
## 127 127 -73.24944 7.817208 39.64945
## 128 128 -73.13985 5.857754 62.71442
## 129 129 -74.28615 5.802960 37.03570
## 130 130 -72.83519 7.326249 61.17689
## 131 131 -73.27355 7.194742 18.38232
## 132 132 -72.82862 8.045154 20.78201
## 133 133 -72.93601 8.106524 50.43324
## 134 134 -74.47026 7.150906 47.61388
## 135 135 -73.01272 8.014469 38.82798
## 136 136 -72.68396 7.508167 53.38914
## 137 137 -74.06478 7.225427 32.01427
## 138 138 -72.46478 7.201317       NA
## 139 139 -73.91574 6.335563 31.34085
## 140 140 -73.37218 7.028166 20.22351
## 141 141 -73.93108 7.275838 31.40185
## 142 142 -74.50752 8.003510       NA
## 143 143 -73.74697 7.731729 34.80710
## 144 144 -73.24067 7.749263 28.10692
## 145 145 -74.36506 6.846248 30.51407
## 146 146 -73.09382 6.675289 28.22809
## 147 147 -74.49437 5.934467 35.40536
## 148 148 -72.73218 7.464331 38.00966
## 149 149 -72.49108 5.884056       NA
## 150 150 -72.44286 6.252275       NA
## 151 151 -73.62423 7.012824 20.11571
## 152 152 -72.54807 7.992551 18.44113
## 153 153 -73.97273 7.786523 17.23624
## 154 154 -73.38971 6.701591 74.97483
## 155 155 -73.47300 6.927344 21.92573
## 156 156 -74.18095 6.756385 32.18096
## 157 157 -73.77985 6.469262 22.71301
## 158 158 -74.18095 8.113099 29.14473
## 159 159 -72.58533 6.381590 39.63258
## 160 160 -72.58094 6.074741 51.78520
## 161 161 -73.61985 7.096112 23.41284
## 162 162 -72.94697 7.402961 52.69487
## 163 163 -72.44067 7.221043       NA
## 164 164 -74.14149 6.909810 30.33812
## 165 165 -74.20944 5.851179 29.68117
## 166 166 -73.76231 7.495016 25.11459
## 167 167 -74.03190 6.447344 34.20715
## 168 168 -73.00396 7.681318 47.18099
## 169 169 -74.29930 7.525701 58.04691
## 170 170 -73.56067 7.536660 35.02535
## 171 171 -73.49492 6.855016 21.91438
## 172 172 -73.31081 6.245700 38.86548
## 173 173 -74.33218 6.855016 20.16607
## 174 174 -73.19245 7.799674 32.85225
## 175 175 -72.96012 6.576659 31.64868
## 176 176 -74.23136 7.080769 53.57404
## 177 177 -73.52341 6.837481 20.29269
## 178 178 -73.20341 7.948715 24.80415
## 179 179 -74.37821 7.878578 54.77054
## 180 180 -74.50971 6.535015 45.95883
## 181 181 -72.80012 6.008987 38.18554
## 182 182 -74.32341 6.699399 30.22901
## 183 183 -72.67519 7.223235 48.45990
## 184 184 -73.45547 6.594193 38.47881
## 185 185 -73.51245 7.975017 21.18150
## 186 186 -73.27136 7.221043 21.27345
## 187 187 -73.12451 6.282960 35.77753
## 188 188 -74.06478 7.295564 25.73138
## 189 189 -72.84177 5.787617 59.86628
## 190 190 -74.04725 7.433646 24.15568
## 191 191 -73.84122 6.006796 40.95007
## 192 192 -73.90697 5.798576 52.10886
## 193 193 -73.24286 6.265426 34.01301
## 194 194 -73.14642 7.218852 23.62639
## 195 195 -74.31903 6.701591 30.73307
## 196 196 -73.11136 7.740496 42.27460
## 197 197 -73.73382 6.228166 21.29712
## 198 198 -73.38094 7.576112 29.87044
## 199 199 -73.23410 6.004604 39.08007
## 200 200 -73.94204 6.375015 32.58807
## 201 201 -73.15300 8.036387 25.22655
## 202 202 -72.48012 6.293919       NA
## 203 203 -73.37656 6.302686 29.41260
## 204 204 -72.51519 5.822686       NA
## 205 205 -72.87903 7.793099 39.94248
## 206 206 -72.99081 6.208440 75.05328
## 207 207 -72.48889 7.953099 18.66853
## 208 208 -74.55355 7.278030       NA
## 209 209 -72.76286 8.108715 15.08638
## 210 210 -73.76670 7.628715 40.50224
## 211 211 -74.44177 5.960768 19.71934
## 212 212 -72.49546 6.464878       NA
## 213 213 -72.46040 6.918577       NA
## 214 214 -73.65492 6.734467 25.26451
## 215 215 -72.78040 5.820494 49.83630
## 216 216 -72.97546 7.146523 52.84117
## 217 217 -73.31081 7.512550 17.93532
## 218 218 -72.83300 6.968988 71.14563
## 219 219 -73.62642 6.688440 33.92425
## 220 220 -72.70588 6.050631 38.64185
## 221 221 -74.03410 7.589263 22.50572
## 222 222 -74.03848 6.822139       NA
## 223 223 -72.49327 7.650633 25.58686
## 224 224 -73.18588 7.602414 29.49251
## 225 225 -74.44177 6.644604 34.76703
## 226 226 -74.11958 6.469262 32.33146
## 227 227 -73.75355 7.486249 28.09793
## 228 228 -73.60889 7.999126 12.85587
## 229 229 -73.10916 5.750357 36.14986
## 230 230 -72.60725 6.876933 54.25474
## 231 231 -73.65273 6.101042 54.50092
## 232 232 -73.46862 5.947617 35.34071
## 233 233 -74.04944 6.596385 32.75922
## 234 234 -73.89382 7.854469 21.92833
## 235 235 -73.64615 6.392549 26.26751
## 236 236 -73.67684 5.921316 38.00555
## 237 237 -72.44724 6.653372       NA
## 238 238 -73.75355 6.158029 21.99532
## 239 239 -74.16341 7.416112 33.82071
## 240 240 -74.45053 7.365701 40.82557
## 241 241 -72.82862 5.949809 34.91240
## 242 242 -74.17656 6.633645 23.06605
## 243 243 -74.08670 5.938850 32.96878
## 244 244 -74.43300 7.975017 32.79679
## 245 245 -72.86807 7.672551 67.88184
## 246 246 -72.91190 7.505975 61.08981
## 247 247 -73.22094 6.127344 34.79502
## 248 248 -73.25820 6.688440 33.13253
## 249 249 -73.48834 6.171179 27.43503
## 250 250 -74.47464 6.385974 44.98693
## 251 251 -73.13766 7.545427 36.29713
## 252 252 -73.33711 5.993645 36.78954
## 253 253 -74.43300 6.944879 62.17144
## 254 254 -72.55683 6.887892 44.53963
## 255 255 -73.27136 7.218852 20.65205
## 256 256 -72.83081 6.914194 46.92439
## 257 257 -72.57437 7.893921 18.91816
## 258 258 -73.77327 7.161865 19.46141
## 259 259 -73.45766 8.095565 15.27970
## 260 260 -72.79793 7.915838 20.17726
## 261 261 -74.46807 7.861044       NA
## 262 262 -73.19464 7.505975 27.74461
## 263 263 -73.71190 7.190358 19.63759
## 264 264 -74.41766 6.684056 29.72932
## 265 265 -73.28012 6.633645 59.91309
## 266 266 -73.35245 5.807343 56.58055
## 267 267 -72.51738 6.624878       NA
## 268 268 -74.23574 7.015016 28.62822
## 269 269 -73.45547 7.615564 35.72246
## 270 270 -73.75574 7.168440 18.32500
## 271 271 -73.03464 7.554194 52.60732
## 272 272 -73.80177 7.637482 24.96370
## 273 273 -73.69437 5.934467 38.15382
## 274 274 -72.67300 6.857207 35.92673
## 275 275 -72.64012 6.050631 46.75638
## 276 276 -73.64396 7.578304 19.39101
## 277 277 -74.12396 6.085700 26.06566
## 278 278 -73.17711 5.952001 54.96245
## 279 279 -73.85875 7.863236 19.84436
## 280 280 -74.16999 6.681865 26.11835
## 281 281 -73.60451 6.758577 20.87590
## 282 282 -74.00560 7.560770 24.66288
## 283 283 -74.14149 7.424879 25.80301
## 284 284 -72.89218 5.862138 37.84090
## 285 285 -74.46369 6.734467 43.98751
## 286 286 -72.58971 6.960221 76.82034
## 287 287 -73.97930 6.655563 23.52059
## 288 288 -74.33437 6.688440 24.56556
## 289 289 -73.17930 5.833645 45.09869
## 290 290 -73.75793 6.850632 36.84744
## 291 291 -73.03245 6.846248 32.84698
## 292 292 -73.29327 7.622140 27.67715
## 293 293 -72.84396 6.350905 56.88662
## 294 294 -72.92725 6.219398 53.24379
## 295 295 -73.86094 7.790907 25.28194
## 296 296 -72.73656 5.899398 55.06007
## 297 297 -73.23629 5.993645 43.28880
## 298 298 -73.37437 7.983784 21.20854
## 299 299 -73.54094 7.374468 19.57732
## 300 300 -73.01711 7.168440 31.42556
## 301 301 -73.53875 6.447344 26.99893
## 302 302 -74.12396 7.054468 27.89057
## 303 303 -73.27355 7.690085 30.71265
## 304 304 -72.69272 6.688440 46.10737
## 305 305 -74.43958 6.158029 30.13487
## 306 306 -72.52615 6.094467       NA
## 307 307 -73.45547 7.889537 18.47433
## 308 308 -73.01711 7.376660 46.58370
## 309 309 -72.59190 8.056113 16.05432
## 310 310 -73.12451 5.774467 43.23463
## 311 311 -73.46423 8.007893 13.38861
## 312 312 -73.78862 7.175016 17.89427
## 313 313 -73.99026 6.272001 27.29986
## 314 314 -73.09163 5.910357 48.51468
## 315 315 -74.37163 6.521865 23.65972
## 316 316 -72.99300 6.438577 30.54139
## 317 317 -74.22697 6.114193 20.85888
## 318 318 -72.85492 6.995290 60.71124
## 319 319 -73.26916 6.206248 40.52767
## 320 320 -72.72122 6.469262 39.04755
## 321 321 -73.97711 7.411728 21.14436
## 322 322 -73.63957 7.652825 20.20422
## 323 323 -73.21437 7.799674 34.57818
## 324 324 -72.94916 7.523509 44.16721
## 325 325 -72.44505 6.342138       NA
## 326 326 -73.25163 5.737206 37.29203
## 327 327 -72.44505 7.619948       NA
## 328 328 -73.57163 7.519126 36.47328
## 329 329 -73.87190 6.236933 27.83380
## 330 330 -73.60451 7.238578 19.27549
## 331 331 -73.71190 7.580496 30.59500
## 332 332 -73.53656 7.394194 15.30825
## 333 333 -74.04944 5.868713 23.98436
## 334 334 -73.59793 6.914194 26.30350
## 335 335 -72.64012 6.272001 28.43031
## 336 336 -72.55245 7.839126 20.43381
## 337 337 -73.49273 7.207893 23.26634
## 338 338 -73.26478 7.367893 20.50091
## 339 339 -72.92505 7.356934 51.56438
## 340 340 -72.82642 7.876386 27.55520
## 341 341 -72.73875 7.872003 20.18484
## 342 342 -73.39629 7.591455 20.71361
## 343 343 -73.59793 7.273646 27.05190
## 344 344 -73.50807 5.895015 51.15445
## 345 345 -72.59409 7.926797 11.93713
## 346 346 -73.02588 6.101042 49.23828
## 347 347 -72.54149 7.118029 76.35188
## 348 348 -73.98807 7.534468 22.35123
## 349 349 -73.92231 6.640221 29.83994
## 350 350 -73.71848 8.018852       NA
## 351 351 -73.48396 6.916385 25.71694
## 352 352 -72.74314 6.024330 37.32923
## 353 353 -73.51026 6.519673 33.12961
## 354 354 -73.36999 7.661592 17.96783
## 355 355 -73.34368 6.550358 33.40836
## 356 356 -72.88341 7.061043 54.24199
## 357 357 -73.15300 6.804605 22.26552
## 358 358 -73.76012 5.770083 82.36882
## 359 359 -72.69492 7.878578 22.68740
## 360 360 -73.31957 6.017754 42.00903
## 361 361 -73.20560 7.552003 22.20754
## 362 362 -73.92451 6.239125 24.37524
## 363 363 -73.91793 6.359673 32.61518
## 364 364 -73.59355 7.948715 14.39087
## 365 365 -74.25766 8.040770 27.03939
## 366 366 -74.15464 7.291180 68.87611
## 367 367 -72.44944 7.532276       NA
## 368 368 -74.10642 6.482412 36.94728
## 369 369 -74.09108 7.633098 23.86676
## 370 370 -74.28396 8.073647 48.78477
## 371 371 -73.79300 6.971180 16.96653
## 372 372 -72.83519 6.456111 42.14939
## 373 373 -73.90916 6.116385 20.78928
## 374 374 -73.56725 7.324057 20.82701
## 375 375 -72.50204 7.924606       NA
## 376 376 -73.37875 6.105426 42.71332
## 377 377 -74.23355 5.730631 48.94191
## 378 378 -72.69492 7.063235 62.50469
## 379 379 -72.78259 5.949809 39.11405
## 380 380 -74.52725 7.793099       NA
## 381 381 -73.79738 6.870358 27.10026
## 382 382 -72.64012 7.887345 19.39812
## 383 383 -73.18368 7.058851 28.30287
## 384 384 -73.10259 7.348167 29.14145
## 385 385 -73.27574 7.262687 21.25363
## 386 386 -73.20560 7.299947 19.63104
## 387 387 -74.04286 5.912549 27.12214
## 388 388 -73.92231 6.552549 34.86199
## 389 389 -73.52122 7.212276 27.89954
## 390 390 -74.28396 7.968441 33.52697
## 391 391 -73.81492 6.140494 44.98314
## 392 392 -74.22697 7.815016 32.45933
## 393 393 -72.50861 6.002412       NA
## 394 394 -74.41327 6.298303 36.96074
## 395 395 -72.52396 7.552003 42.03910
## 396 396 -72.79574 6.276385 37.19548
## 397 397 -73.02149 6.828714 32.12154
## 398 398 -73.93108 6.776111 42.30256
## 399 399 -72.96231 6.712550 23.07893
## 400 400 -73.58478 5.844604 40.54174
## 401 401 -74.09327 8.014469 30.66399
## 402 402 -74.24232 6.212823 18.44606
## 403 403 -72.79793 7.348167 55.63739
## 404 404 -73.57163 7.972825 15.79856
## 405 405 -73.65930 6.173371 27.23950
## 406 406 -73.91574 6.458303 21.51831
## 407 407 -72.62478 6.138303 57.84774
## 408 408 -74.22040 7.203509 40.07147
## 409 409 -73.49273 6.570084 40.97471
## 410 410 -73.49492 6.410084 63.64931
## 411 411 -73.57601 7.028166 42.76866
## 412 412 -73.58478 8.071455 11.99131
## 413 413 -74.11081 7.788715 25.35094
## 414 414 -72.57656 6.839673 55.43125
## 415 415 -73.03683 5.873097 45.66635
## 416 416 -72.54807 5.829261       NA
## 417 417 -74.33218 6.620495 25.32955
## 418 418 -74.49875 5.984878 25.22976
## 419 419 -72.51081 6.890084       NA
## 420 420 -74.01218 6.725700 28.35497
## 421 421 -74.06697 5.846795 25.35579
## 422 422 -73.28012 6.758577 39.72438
## 423 423 -73.54752 7.768989 27.71651
## 424 424 -74.11519 5.938850 45.71989
## 425 425 -74.04725 6.916385 42.33695
## 426 426 -74.19190 6.160220 30.91815
## 427 427 -73.84779 6.613919 32.53971
## 428 428 -72.90094 7.604605 56.13659
## 429 429 -73.58916 8.034195 17.68919
## 430 430 -72.49766 6.879125       NA
## 431 431 -74.43738 6.936112 42.18103
## 432 432 -73.63519 6.574467 38.24611
## 433 433 -73.81273 6.269810 36.23250
## 434 434 -73.84779 7.900496 24.95307
## 435 435 -73.11574 5.914741 51.27871
## 436 436 -73.19026 6.101042 30.76420
## 437 437 -74.15464 6.767344 31.31685
## 438 438 -73.74697 6.673098 25.25042
## 439 439 -73.31738 5.993645 37.28141
## 440 440 -74.04725 5.848987 27.84288
## 441 441 -72.81546 6.962413 63.68837
## 442 442 -74.31684 6.204056 17.37730
## 443 443 -74.28177 5.824878 36.36651
## 444 444 -73.72505 7.251728       NA
## 445 445 -73.45985 5.982686 45.99365
## 446 446 -74.36506 5.892823 36.63497
## 447 447 -73.72944 8.047345       NA
## 448 448 -72.61382 8.051729 13.79646
## 449 449 -74.42862 6.307070 27.47002
## 450 450 -73.04341 7.957482 39.27979
## 451 451 -74.03629 7.505975 27.86212
## 452 452 -73.56067 6.361864 35.24495
## 453 453 -74.42204 6.068165 27.48841
## 454 454 -73.60451 6.188714 42.94498
## 455 455 -73.51245 6.421043 35.89798
## 456 456 -73.17053 7.501591 32.06793
## 457 457 -72.69272 7.271454 41.28005
## 458 458 -73.10259 7.223235 30.93579
## 459 459 -73.47081 6.602961 35.47165
## 460 460 -73.66588 6.515289 38.14496
## 461 461 -74.13930 7.238578 34.02181
## 462 462 -72.57437 6.272001 46.27254
## 463 463 -72.85711 6.289536 54.77951
## 464 464 -73.97930 5.820494 37.48809
## 465 465 -73.58040 6.153645 34.75377
## 466 466 -72.91190 7.424879 63.65310
## 467 467 -74.23355 8.069263 35.70585
## 468 468 -74.46588 6.212823 28.45226
## 469 469 -72.74314 7.091728 57.44058
## 470 470 -74.09766 6.609536 27.26381
## 471 471 -73.62642 7.383235 32.42111
## 472 472 -74.52944 7.262687       NA
## 473 473 -73.08725 6.068165 46.72138
## 474 474 -73.41382 6.811180 32.96577
## 475 475 -74.28396 6.254467 23.06467
## 476 476 -73.56067 6.828714 31.47260
## 477 477 -72.88560 6.938303 48.61845
## 478 478 -74.00999 6.495563 21.80342
## 479 479 -73.10916 8.091181 26.49287
## 480 480 -73.63519 6.892276 25.83645
## 481 481 -73.46862 5.927891 40.45108
## 482 482 -74.14368 6.276385 20.08639
## 483 483 -72.46259 7.442413       NA
## 484 484 -73.26697 7.957482 38.73718
## 485 485 -74.47026 6.513097 38.56111
## 486 486 -74.09547 7.470907 19.70577
## 487 487 -74.21821 5.967343 23.74341
## 488 488 -72.74314 7.324057 43.71346
## 489 489 -74.31245 7.696660 38.13279
## 490 490 -72.55245 6.423234 50.88673
## 491 491 -72.96670 6.044056 56.00127
## 492 492 -74.39574 8.007893 47.85434
## 493 493 -73.26259 7.317482 20.80192
## 494 494 -72.58752 7.475290 48.84172
## 495 495 -73.97711 5.735014 26.44505
## 496 496 -73.87410 6.909810 29.24789
## 497 497 -73.81053 7.819400 35.32628
## 498 498 -73.15957 8.067071 25.93486
## 499 499 -73.30423 7.212276 25.11077
## 500 500 -72.57218 7.361317 74.64268

Eliminemos los valores de NA:

sitios <- na.omit(sitios)
head(sitios)
##   id    longit    latit      soc
## 1  1 -73.85437 5.745973 42.45292
## 2  2 -73.48834 5.785425 45.15730
## 3  3 -73.97930 6.563508 19.95003
## 4  4 -73.03464 7.722962 41.47083
## 5  5 -73.29985 6.495563 25.04777
## 6  6 -73.70533 8.115291 12.87853

Visualicemos las muestras:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("yellow", "orange", "green", "blue"), 
                                  domain = 8:130)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())

m  # Print the map

5. Interpolación

##5.1 Creación del objeto gstat Para interpolar, primero necesitamos crear un objeto de clase gstat, usando una función del mismo nombre: gstat.

Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial:

-La definición del modelo. -Los datos de calibración

Según sus argumentos, la función gstat “entiende” qué tipo de modelo de interpolación queremos usar:

-Sin modelo de variograma → IDW - Modelo de variograma, sin covariables → Kriging ordinario

Vamos a utilizar tres parámetros de la función gstat:

-fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes (también conocidas como covariables) -datos: los datos de calibración (también conocidos como los datos del tren) -modelo: el modelo de variograma

5.2 interpolación IDW

Para interpolar SOC usando el método IDW creamos el siguiente objeto gstat, especificando solo la fórmula y los datos:

g1 = gstat(formula = soc_igh_15_30 ~ 1, data = nmuestras)

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función de predicción para interpolar, es decir, estimar los valores de precipitación.

La función de predicción acepta:

-Un objeto ráster: star, como dem -Un modelo: objeto gstat, como g1

El ráster sirve para dos propósitos:

-Especificar las ubicaciones donde queremos hacer predicciones (en todos los métodos) -Especificación de valores de covariables (solo en Universal Kriging)

Creemos un objeto ráster con valores de celda iguales a 1:

# una copia sencilla
rrr = aggregate(geog.soc, 4)

¿Qué es rrr?

rrr
## class       : SpatRaster 
## dimensions  : 274, 243, 1  (nrow, ncol, nlyr)
## resolution  : 0.008767132, 0.008767132  (x, y)
## extent      : -74.55464, -72.42423, 5.718576, 8.12077  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      11.34278 
## max value   :      95.73760

Definir nuevos valores:

values(rrr) <-1

Definir nuevos nombres:

names(rrr) <- "valor"

¿Qué es rrr ahora?

rrr
## class       : SpatRaster 
## dimensions  : 274, 243, 1  (nrow, ncol, nlyr)
## resolution  : 0.008767132, 0.008767132  (x, y)
## extent      : -74.55464, -72.42423, 5.718576, 8.12077  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1
stars.rrr = st_as_stars(rrr)

Por ejemplo, la siguiente expresión interpola los valores SOC según el modelo definido en g1 y la plantilla ráster definida en stars.sntdr:

z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

¿Qué es z1?

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  12.03249 29.26843 34.02578 35.22126 40.78655 92.57333     0
## var1.var         NA       NA       NA      NaN       NA       NA 66582
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 243 -74.55  0.008767 +proj=longlat +datum=WGS8... [x]
## y    1 274  8.121 -0.008767 +proj=longlat +datum=WGS8... [y]

Tome nota de los nombres de los dos atributos incluidos dentro del objeto z1.

  • var1.pred
  • var1.var

Podemos crear un subconjunto solo del primer atributo y cambiarle el nombre a “soc”:

z1 = z1["var1.pred"]
names(z1) = "soc"

Necesitamos una paleta de colores:

paleta <- colorNumeric(palette = c("yellow", "orange", "green", "blue"), domain = 10:100, na.color = "transparent")

El ráster SOC interpolado, utilizando IDW, se muestra en la siguiente figura:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("yellow", "orange", "green", "blue"), 
                                  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 [%]"
    )
m # Imprimir el mapa

5.3 Interpolación correcta

Los métodos de Kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos y asignar ponderaciones en consecuencia al realizar predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función de variograma.

La función requiere dos argumentos:

-fórmula: especifica la variable dependiente y las covariables, como en gstat -datos: la capa de puntos con la variable dependiente y covariables como atributos de puntos

Por ejemplo, la siguiente expresión calcula el variograma empírico de muestras, sin covariables:

v_emp_ok = variogram(soc_igh_15_30 ~ 1, data=nmuestras)

Tracemos el variograma:

plot(v_emp_ok)

Hay varias formas de ajustar un modelo de variograma a un variograma empírico. Usaremos el más simple: ajuste automático usando la función autofitVariogram del paquete automap:

v_mod_ok = autofitVariogram(soc_igh_15_30 ~ 1, as(nmuestras, "Spatial"))

La función elige el tipo de modelo que mejor se adapta y también ajusta sus parámetros. Puede utilizar show.vgms() para mostrar los tipos de modelos de variograma.

show.vgms()

Tenga en cuenta que la función autofitVariogram no funciona en objetos sf, por lo que convertimos el objeto en un SpatialPointsDataFrame (paquete sp).

El modelo ajustado se puede trazar con plot:

plot(v_mod_ok)

El objeto resultante es en realidad una lista con varios componentes, incluido el variograma empírico y el modelo de variograma ajustado. El componente $var_model del objeto resultante contiene el modelo real:

v_mod_ok$var_model
##   model     psill    range kappa
## 1   Nug  12.72716  0.00000   0.0
## 2   Ste 200.47247 80.59171   0.3

En tu cuaderno, explica el significado de cada elemento del modelo anterior.

Modelo: El modelo de variograma es una función que describe la relación entre la varianza de los datos y la distancia entre dichos objetos, esto quiere decir que es una forma de modelar la autocorrelación espacial de los datos. Para este caso el modelo usado fue Ste

Pepita: Es la varianza de los datos a una distancia de cero, en este caso es 13

Sill: Es la varianza que se puede explicar por la autocorrelación espacial, para este cado es de 213.

Rango: Es una medida de la escala de la autocorrelación espacial. Para este caso es 81. Esto significa que la autocorrelación espacial se reduce a cero a una distancia de 81.

Kappa: Este controla la suavidad de la función de variograma. Un valor más alto de kappa indica una función de variograma más suave. Para este caso es 0.3 lo que significa que es suave. Esto quiere decir que la autocorrelación espacial de los datos disminuye más lentamente con la distancia.

Ahora, el modelo de variograma se puede pasar a la función gstat y podemos continuar con la interpolación Kriging ordinaria:

g2 = gstat(formula = soc_igh_15_30 ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Nuevamente, crearemos un subconjunto del atributo de valores predichos y le cambiaremos el nombre:

z2 = z2["var1.pred",,]
names(z2) = "soc"

Las predicciones de Kriging ordinario se muestran en la siguiente figura:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("yellow", "orange", "green", "blue"), 
                                  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 [%]"
    )
m # Imprimir el mapa

#6. Evaluación de resultados

6.1 Evaluación cualitativa

Otra vista de las tres salidas de interpolación:

colores <- colorOptions(palette = c("yellow", "orange", "green", "blue"), 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 [%]"
)
m # Imprimir el mapa

6.2 Validación cruzada

Hemos estimado las superficies climáticas utilizando dos métodos diferentes: IDW y Kriging Ordinario. Aunque es útil examinar y comparar los resultados gráficamente, también necesitamos una forma objetiva de evaluar la precisión de la interpolación. De esa manera, podemos elegir objetivamente el método más preciso entre los métodos de interpolación disponibles.

En la validación cruzada Leave-One-Out nosotros:

-Saque un punto de los datos de calibración. -Haz una predicción para ese punto. -Repita para todos los puntos. Al final, lo que obtenemos es una tabla con un valor observado y un valor predicho para todos los puntos.

Podemos ejecutar la validación cruzada Leave-One-Out usando la función gstat.cv, que acepta un objeto gstat.

Al escribir el siguiente fragmento, oculte el mensaje y los resultados.

### ejecuta el siguiente código desde la consola, línea por línea
cv1 = gstat.cv(g1)
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
cv2 = gstat.cv(g2)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]

El resultado de gstat.cv tiene los siguientes atributos:

-var1.pred: valor previsto -var1.var: varianza (solo para Kriging) -observado: valor observado -residual: observado-predicho -zscore: puntuación Z (solo para Kriging) -plegar: ID de validación cruzada

cv1 = na.omit(cv1)
cv1
## class       : SpatialPointsDataFrame 
## features    : 465 
## extent      : -74.50971, -72.48889, 5.730631, 8.115291  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 17.5724595411129,       NA, 11.9371318817139, -17.5158108636536,     NA,    1 
## max values  : 62.5729664727844,       NA, 92.8670425415039,  52.3412983970521,     NA,  465

Convirtamos el objeto cv1:

cv1 = st_as_sf(cv1)

Ahora, grafiquemos los residuos:

sp::bubble(as(cv1[, "residual"], "Spatial"))

Ahora, calculamos índices de precisión de predicción, como el error cuadrático medio (RMSE):

# Este es el valor RMSE para la interpolación IDW
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 9.467968

Ahora, repita el proceso con los resultados correctos:

Tiempo de conversión:

cv2 = st_as_sf(cv2)

Calcule RSME para obtener resultados correctos:

# Este es el valor RMSE para la interpolación OK
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 8.895876

Escribe un párrafo que resuma tus resultados. Según ellos, qué técnica de interpolación parece más precisa. Explicar.

interpolación IDW: La interpolación mediante distancia inversa ponderada (IDW) es una técnica que se utiliza para interpolar valores desconocidos de una variable espacial a partir de valores conocidos. IDW supone explícitamente que los elementos cercanos entre sí son más parecidos que los que están más alejados. Para predecir un valor para cualquier ubicación sin medición, IDW usa los valores medidos que rodean a la ubicación de predicción. Dio como resultado en el valor RMSE 9.467968

Kriging ordinario: En este método se utiliza la media de un subconjunto de puntos vecinos para producir un punto de interpolación específico. Es flexible y puede adaptarse a cambios en el valor medio de la superficie, siempre que la distancia entre vecinos de búsqueda no sea demasiado grande. Dio como resultado en el valor RMSE 8.895876

Segun los resultados de error Cuadrático Medio de la Raíz (RMSE) podemos decir que la tecnica de interpolación mas precisa es Kriging ordinario ya que, entre menor sea el resultado es mas bajo y esto sugiere una mejor capacidad de predicción.

7. Referencia

Citar este trabajo como: Lizarazo, I., 2023. Interpolación espacial del carbono orgánico del suelo. Disponible en https://rpubs.com/ials2un/soc_interp.

sessionInfo()
## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22621)
## 
## 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] curl_5.1.0    leafem_0.2.3  automap_1.1-9 gstat_2.1-1   dplyr_1.1.4  
##  [6] ggplot2_3.4.4 leaflet_2.2.1 stars_0.6-4   abind_1.4-5   sf_1.0-14    
## [11] terra_1.7-55  sp_2.1-1     
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       xfun_0.41          bslib_0.6.0        raster_3.6-26     
##  [5] htmlwidgets_1.6.3  lattice_0.21-9     vctrs_0.6.4        tools_4.3.2       
##  [9] crosstalk_1.2.1    generics_0.1.3     parallel_4.3.2     tibble_3.2.1      
## [13] proxy_0.4-27       spacetime_1.3-0    fansi_1.0.5        highr_0.10        
## [17] xts_0.13.1         pkgconfig_2.0.3    KernSmooth_2.23-22 lifecycle_1.0.4   
## [21] farver_2.1.1       compiler_4.3.2     FNN_1.1.3.2        munsell_0.5.0     
## [25] codetools_0.2-19   htmltools_0.5.7    class_7.3-22       sass_0.4.7        
## [29] yaml_2.3.7         pillar_1.9.0       jquerylib_0.1.4    ellipsis_0.3.2    
## [33] classInt_0.4-10    cachem_1.0.8       tidyselect_1.2.0   digest_0.6.33     
## [37] fastmap_1.1.1      grid_4.3.2         colorspace_2.1-0   cli_3.6.1         
## [41] magrittr_2.0.3     base64enc_0.1-3    utf8_1.2.4         e1071_1.7-13      
## [45] withr_2.5.2        scales_1.2.1       rmarkdown_2.25     zoo_1.8-12        
## [49] png_0.1-8          evaluate_0.23      knitr_1.45         rlang_1.1.2       
## [53] Rcpp_1.0.11        glue_1.6.2         DBI_1.1.3          reshape_0.8.9     
## [57] jsonlite_1.8.7     R6_2.5.1           plyr_1.8.9         intervals_0.15.4  
## [61] units_0.8-4