TALLER INTERPOLACION

Como primer paso, se deben analizar y descargar las librerías necesarias para realizar la interpolación..

install.packages( c( "sf", "stars","leaflet", "gstat","automap","raster", "RColorBrewer"),repos = "https://cran.r-project.org")

## Installing packages into 'C:/Users/natal/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)

## package 'sf' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'sf'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problema al copiar
## C:\Users\natal\AppData\Local\R\win-library\4.3\00LOCK\sf\libs\x64\sf.dll a
## C:\Users\natal\AppData\Local\R\win-library\4.3\sf\libs\x64\sf.dll: Permission
## denied

## Warning: restored 'sf'

## package 'stars' successfully unpacked and MD5 sums checked
## package 'leaflet' successfully unpacked and MD5 sums checked
## package 'gstat' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'gstat'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problema al copiar
## C:\Users\natal\AppData\Local\R\win-library\4.3\00LOCK\gstat\libs\x64\gstat.dll
## a C:\Users\natal\AppData\Local\R\win-library\4.3\gstat\libs\x64\gstat.dll:
## Permission denied

## Warning: restored 'gstat'

## package 'automap' successfully unpacked and MD5 sums checked
## package 'raster' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'raster'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problema al copiar
## C:\Users\natal\AppData\Local\R\win-library\4.3\00LOCK\raster\libs\x64\raster.dll
## a C:\Users\natal\AppData\Local\R\win-library\4.3\raster\libs\x64\raster.dll:
## Permission denied

## Warning: restored 'raster'

## package 'RColorBrewer' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\natal\AppData\Local\Temp\Rtmp299Law\downloaded_packages

paquetes = c('terra', 'sp', 'leafem')
install.packages(paquetes, repos = "https://cran.r-project.org")

## Installing packages into 'C:/Users/natal/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)

## package 'terra' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'terra'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problema al copiar
## C:\Users\natal\AppData\Local\R\win-library\4.3\00LOCK\terra\libs\x64\terra.dll
## a C:\Users\natal\AppData\Local\R\win-library\4.3\terra\libs\x64\terra.dll:
## Permission denied

## Warning: restored 'terra'

## package 'sp' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'sp'

## Warning in file.copy(savedcopy, lib, recursive = TRUE): problema al copiar
## C:\Users\natal\AppData\Local\R\win-library\4.3\00LOCK\sp\libs\x64\sp.dll a
## C:\Users\natal\AppData\Local\R\win-library\4.3\sp\libs\x64\sp.dll: Permission
## denied

## Warning: restored 'sp'

## package 'leafem' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\natal\AppData\Local\Temp\Rtmp299Law\downloaded_packages

Procedemos a verificar las librerías..

library(terra)

## Warning: package 'terra' was built under R version 4.3.3

## terra 1.8.21

library(sp)

## Warning: package 'sp' was built under R version 4.3.3

library(leafem)

## Warning: package 'leafem' was built under R version 4.3.3

library(sf) 

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(stars) 

## Warning: package 'stars' was built under R version 4.3.3

## Loading required package: abind

library(leaflet) 

## Warning: package 'leaflet' was built under R version 4.3.3

library(gstat)

## Warning: package 'gstat' was built under R version 4.3.3

library(automap) 

## Warning: package 'automap' was built under R version 4.3.3

library(raster) 

## Warning: package 'raster' was built under R version 4.3.3

library(RColorBrewer)

Confirmamos que el directorio de trabajo es el correcto utilizando el siguiente comando.

getwd()

## [1] "C:/Users/natal/OneDrive/Documentos/interpolacion/Datos"

Definimos las variables que se utilizarán.

ruta_aoi="./municipios_muestreo.shp" 
ruta_puntos="./puntos_muestreo.shp" 
ruta_raster="./dem_srtm_9377.tif"

Procedemos a leer y visualizar los archivos de datos, específicamente aquellos en formato SHP y TIF, para asegurarnos de que se han cargado correctamente.

raster_dem=read_stars(ruta_raster, RasterIO = list(bands = 1))
    raster_dem

## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##                    Min. 1st Qu. Median    Mean 3rd Qu. Max.
## dem_srtm_9377.tif  -247       0      0 131.274       0 1673
## dimension(s):
##   from   to  offset  delta                       refsys point x/y
## x    1 1754 4966240  29.78 MAGNA-SIRGAS 2018 / Orige... FALSE [x]
## y    1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... FALSE [y]

plot(raster_dem)

crs_raster=st_crs(raster_dem)
crs_raster

## Coordinate Reference System:
##   User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
##   wkt:
## PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
##     BASEGEOGCRS["MAGNA-SIRGAS 2018",
##         DATUM["Marco Geocentrico Nacional de Referencia 2018",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",20046]],
##     CONVERSION["Colombia Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-73,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9992,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",2000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Cadastre, topographic mapping."],
##         AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
##         BBOX[-4.23,-84.77,15.51,-66.87]],
##     ID["EPSG",9377]]

puntos = st_read(ruta_puntos)

## Reading layer `puntos_muestreo' from data source 
##   `C:\Users\natal\OneDrive\Documentos\interpolacion\Datos\puntos_muestreo.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 67 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
## Projected CRS: MAGNA-SIRGAS_Origen-Nacional

crs_puntos=st_crs(puntos)
crs_puntos

## Coordinate Reference System:
##   User input: MAGNA-SIRGAS_Origen-Nacional 
##   wkt:
## PROJCRS["MAGNA-SIRGAS_Origen-Nacional",
##     BASEGEOGCRS["GCS_MAGNA-SIRGAS_2018",
##         DATUM["Marco Geocentrico Nacional de Referencia 2018",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",1329]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-73,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9992,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",2000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

Realizamos una comparación detallada de los sistemas de referencia utilizados, asegurándonos de que sean compatibles para el análisis.

if (identical(crs_raster, crs_puntos)) {
  cat("El sistema de referencia de coordenadas del raster y de los puntos es el mismo:\n")
  print(crs_raster)
} else {
  cat("Los sistemas de referencia de coordenadas son diferentes.\n")
  cat("Raster:\n")
  print(crs_raster)
  cat("Puntos:\n")
  print(crs_puntos)
}

## Los sistemas de referencia de coordenadas son diferentes.
## Raster:
## Coordinate Reference System:
##   User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
##   wkt:
## PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
##     BASEGEOGCRS["MAGNA-SIRGAS 2018",
##         DATUM["Marco Geocentrico Nacional de Referencia 2018",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",20046]],
##     CONVERSION["Colombia Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-73,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9992,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",2000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Cadastre, topographic mapping."],
##         AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
##         BBOX[-4.23,-84.77,15.51,-66.87]],
##     ID["EPSG",9377]]
## Puntos:
## Coordinate Reference System:
##   User input: MAGNA-SIRGAS_Origen-Nacional 
##   wkt:
## PROJCRS["MAGNA-SIRGAS_Origen-Nacional",
##     BASEGEOGCRS["GCS_MAGNA-SIRGAS_2018",
##         DATUM["Marco Geocentrico Nacional de Referencia 2018",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",1329]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-73,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9992,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",2000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

puntos_reproj <- st_transform(puntos, crs_raster)

(puntos_reproj)

## Simple feature collection with 67 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
##    Department Municipali Latitude.D Longitude. Altitude.m  pH EC.dSm NH4.ppm
## 1   Santander   Confines   6.331580  -73.25827       1533 4.8   0.08    10.3
## 2   Santander   Confines   6.345463  -73.28759       1489 4.9   0.07    10.3
## 3   Santander   Confines   6.357339  -73.28138       1482 4.8   0.08    18.2
## 4   Santander   Confines   6.360368  -73.24903       1581 4.7   0.07     6.7
## 5   Santander   Confines   6.340922  -73.20868       1878 4.8   0.15    24.0
## 6   Santander   Confines   6.357226  -73.23005       1742 4.8   0.09    11.7
## 7   Santander   Confines   6.339062  -73.23301       1752 4.9   0.09     9.4
## 8   Santander   Confines   6.402768  -73.22524       1715 4.8   0.08     9.4
## 9   Santander   Confines   6.365608  -73.21098       1878 4.8   0.10    19.3
## 10  Santander     Paramo   6.384360  -73.19552       1912 4.9   0.12    19.1
##    NO3.ppm K2O.ppm P2O5.ppm SOC.pct Sand.pct Silt.pct Clay.pct Slope.pct
## 1      8.4   173.5     10.1    2.46     40.7     32.6     26.6      10.5
## 2      5.4   122.9      8.9    1.62     51.2     26.1     22.6       9.2
## 3      0.9   119.3      6.2    2.95     51.9     35.4     12.7       7.4
## 4      2.3   181.9     12.1    1.02     37.3     36.1     26.5      12.6
## 5     12.3   130.1      7.3    2.48     65.6     30.3      4.1       8.5
## 6      1.2   131.3      4.8    1.60     40.8     46.5     12.6       7.5
## 7      2.0    73.5     11.5    3.65     55.0     42.5      2.5       5.6
## 8      5.9   194.0      6.6    2.12     33.7     45.4     20.9      13.7
## 9      0.3   273.5      3.4    3.81     64.4     24.9     10.6       1.5
## 10     0.1   128.9     18.8    3.14     39.1     50.7     10.3      15.8
##                   geometry
## 1  POINT (4971447 2257635)
## 2  POINT (4968206 2259171)
## 3  POINT (4968893 2260483)
## 4  POINT (4972470 2260816)
## 5  POINT (4976930 2258665)
## 6  POINT (4974568 2260467)
## 7  POINT (4974240 2258460)
## 8  POINT (4975101 2265500)
## 9  POINT (4976677 2261393)
## 10 POINT (4978386 2263464)

#Procedemos a interpolar los valores de pH utilizando el método de ponderación por distancia inversa (IDW)

var="pH"
puntos[var]

## Simple feature collection with 67 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 4968206 ymin: 2257635 xmax: 5005157 ymax: 2292361
## Projected CRS: MAGNA-SIRGAS_Origen-Nacional
## First 10 features:
##     pH                geometry
## 1  4.8 POINT (4971447 2257635)
## 2  4.9 POINT (4968206 2259171)
## 3  4.8 POINT (4968893 2260483)
## 4  4.7 POINT (4972470 2260816)
## 5  4.8 POINT (4976930 2258665)
## 6  4.8 POINT (4974568 2260467)
## 7  4.9 POINT (4974240 2258460)
## 8  4.8 POINT (4975101 2265500)
## 9  4.8 POINT (4976677 2261393)
## 10 4.9 POINT (4978386 2263464)

Establecemos y definimos el modelo de interpolación que se utilizará para analizar y estimar los datos de manera precisa.

st_crs(puntos_reproj) 

## Coordinate Reference System:
##   User input: MAGNA-SIRGAS 2018 / Origen-Nacional 
##   wkt:
## PROJCRS["MAGNA-SIRGAS 2018 / Origen-Nacional",
##     BASEGEOGCRS["MAGNA-SIRGAS 2018",
##         DATUM["Marco Geocentrico Nacional de Referencia 2018",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",20046]],
##     CONVERSION["Colombia Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-73,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9992,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",2000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Cadastre, topographic mapping."],
##         AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
##         BBOX[-4.23,-84.77,15.51,-66.87]],
##     ID["EPSG",9377]]

crs(raster_dem)

## [1] NA

g = gstat(formula = pH ~ 1, data = puntos_reproj,set=list(idp=2))
z = predict(g, raster_dem)

## [inverse distance weighted interpolation]

print(z)

## stars object with 2 dimensions and 2 attributes
## attribute(s), summary of first 1e+05 cells:
##               Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  5.04681 5.130046 5.162891 5.174822 5.238671 5.283381 0e+00
## var1.var        NA       NA       NA      NaN       NA       NA 1e+05
## dimension(s):
##   from   to  offset  delta                       refsys x/y
## x    1 1754 4966240  29.78 MAGNA-SIRGAS 2018 / Orige... [x]
## y    1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... [y]

#Dado que el método IDW no genera un atributo de varianza, simplemente extraemos el primer atributo de ‘z’ y lo asignamos a la variable correspondiente.

z_idw = z["var1.pred",,]

#Procedemos a realizar la visualización de los resultados obtenidos.

names(z) = "pH"
b = seq(4, 7.5, 0.1)

Generamos un gráfico para representar visualmente la predicción realizada.

plot(z_idw, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE) 
plot(st_geometry(puntos_reproj), pch = 3, add = TRUE) 
contour(z_idw, breaks = b, add = TRUE)

#Procedemos a añadir los puntos de muestreo al gráfico para complementar la visualización.

#Guardamos el ráster generado en un archivo con formato GeoTIFF para su posterior uso y análisis.

interpolacion_idw="./ph_idw.tif" 

write_stars(z_idw, dsn = interpolacion_idw)

Realizamos la interpolación de los valores de pH utilizando el método de Kriging, que permite obtener estimaciones más precisas considerando la estructura espacial de los datos.

v_emp_ok = variogram(pH ~ 1, puntos_reproj)
plot(v_emp_ok)

#Ajustamos el variograma empírico a una función teórica para modelar la dependencia espacial de los datos.

v_mod_ok = autofitVariogram(pH ~ 1, as(puntos_reproj, "Spatial"))
 plot(v_mod_ok)

Definimos el modelo teórico que se utilizará para realizar la interpolación de los datos.

g2 = gstat(formula = pH ~ 1, model = v_mod_ok$var_model, data = puntos_reproj)

Generamos la predicción de los valores mediante el modelo de interpolación definido.

z_kr = predict(g2, raster_dem)

## [using ordinary kriging]

print(z_kr)

## stars object with 2 dimensions and 2 attributes
## attribute(s), summary of first 1e+05 cells:
##                Min.   1st Qu.    Median      Mean 3rd Qu.     Max.
## var1.pred  4.973535 5.0923192 5.1475910 5.1348409 5.18763 5.238132
## var1.var   0.604284 0.7349258 0.8883815 0.9006311 1.04760 1.302663
## dimension(s):
##   from   to  offset  delta                       refsys x/y
## x    1 1754 4966240  29.78 MAGNA-SIRGAS 2018 / Orige... [x]
## y    1 1428 2297460 -29.78 MAGNA-SIRGAS 2018 / Orige... [y]

Visualizamos la superficie resultante de la predicción para analizar la distribución espacial de los valores estimados.

  prediccion = z_kr["var1.pred",,]
  names(prediccion) = "pH"
  b_predict = seq(4, 7, 0.1)
  plot(prediccion, breaks = b_predict, col = hcl.colors(length(b_predict)-1, "Spectral"), reset = FALSE)

#Procedemos a visualizar la superficie de la varianza para evaluar la incertidumbre asociada a las predicciones.

varianza = z_kr["var1.var",,]
  names(varianza) = "varianza pH"

b_var = seq(0.1, 1.4, 0.1)
 plot(varianza, breaks = b_var, col = hcl.colors(length(b_var)-1, "Spectral"), reset = FALSE) 
 plot(st_geometry(puntos_reproj), pch = 3, add = TRUE)

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]

cv2 = st_as_sf(cv2)

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

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

## [1] 0.7066087

cv3 = gstat.cv(g)

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

cv3 = st_as_sf(cv3)

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

sqrt(sum((cv3$var1.pred - cv3$observed)^2) / nrow(cv3))

## [1] 0.7196336