Introducción

Este archivo rMarkdown tiene como objetivo presentar el taller-parcial propuesto en clase sobre el tema Interpolación espacial en Rstudio con los métodos IDW (Inverse Distance Weighted) y Kriging Ordinario (OK) basado en el cuaderno en Rpubs del profesor Ivan Lizarazo con datos obtenidos de SoilGrids, el cual es un sistema de mapeo digital global de suelos que utiliza información de perfiles de suelos globales y datos de covariables para modelar la distribución espacial de las propiedades del suelo en todo el mundo. SoilGrids es una colección de mapas de propiedades del suelo de todo el mundo producidos mediante aprendizaje automático con una resolución de 250 m. Los datos obtenidos están relacionados con un límite que le hayamos suministrado a SoildGrids, como un shapefile de un municipio. En estas notas se empleó un shapefile del municipio de Boyacá, del departamento con el mismo nombre (Colombia), el cual proviene del filtro y exportación como archivo ESRI shapefile con el sistema de referencia de coordenadas EPSG: 3116 anteriormente usado en Colombia y reemplazado por el EPSG: 9377. Con los datos obtenidos se empleó una función para generar un muestreo aleatorio de 1000 puntos y aplicar los métodos de interpolación espacial por IDW y OK para generar superficies continuas con valores predecidos a partir de valores observados o muestreados como resultado del empleo de estos métodos. Después se calcula la validación cruzada para evaluar el rendimiento de los métodos de interpolación, asegurando que las predicciones sean confiables. Por último, se calculó el Error Cuadrático Medio (RMSE) para cada método, es una métrica utilizada para evaluar la precisión de un modelo predictivo, comparando los valores observados con los reales, nos va a permitir conocer cual es el método más confiable para interpolar datos espaciales.

Actividades

1. Instalar las librerías/paquetes que contienen las funciones que vamos a utilizar

Primero vamos a instalar aquellos paquetes o librerías que no hemos instalado previamente. Si ya las posee instaladas no es necesario volver a instalarlas.

#install.packages("sf", "stars", "leaflet", "gstat", "automap", "raster", "RcolorBrewer", "ggplot2", "leafmen", "tmap", "gdalUtilities", "terra")

2. Traer las librerías/paquetes con contienen las funciones que vamos a utilizar

Ya instaladas vamos a llamar a las librerias correspondientes

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
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.3
library(RColorBrewer)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.3.3
library(tmap)
## Warning: package 'tmap' was built under R version 4.3.3
library(gdalUtilities)
## Warning: package 'gdalUtilities' was built under R version 4.3.3
## 
## Attaching package: 'gdalUtilities'
## The following object is masked from 'package:sf':
## 
##     gdal_rasterize
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21

Vamos a cargar los datos desde R basado en el Notebook del profe desde la fuente de datos de suelos SoilGrids. Primero establecemos el link de donde vamos a tomar los datos

url = "https://files.isric.org/soilgrids/latest/data/" # Path to the webDAV data. #Ruta a la fuente de datos

#indicamos que necesitamos de la fuente

voi= "soc" 
depth= "15-30cm"
quantile= "mean"
(variable= paste(url, voi, sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc"
#Nos generará la dirección de donde obtendremos los datos
(layer = paste(variable,depth,quantile, sep="_")) # Capa de interés
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"

Ahora, la especificación del VRT estará completa

(vrt_layer = paste(layer, '.vrt', sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"

3. Descargar un ráster de SOC para una región de interés

Para esta nota, vamos a utilizar un shapefile del municipio de Boyacá, del departamento del mismo nombre (Colombia) para este caso, proviene del filtro y exportación como archivo ESRI Shapefile con el sistema de coordenadas de referencia EPSG: 3116. El municipio se escogió por un enorme valor familiar y por ser un municipio no muy extenso en área para hacer más breves los cálculos.

Vamos a establecer el mismo CRS de la proyección Homolosine

igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
byaca_igh= st_transform(byaca, igh) #Cambiar el CRS

Ahora tenemos que conseguir el cuadro delimitador de nuestra área de interés

bbox= st_bbox(byaca_igh)

Ahora ya con los límites, vamos a establecer el cuadro de límites tal cual lo establece la librería GDA.

ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
##       xmin       ymax       xmax       ymin 
## -8188601.3   609574.2 -8176948.9   601634.2

Ahora podemos con la función gdal_translate descargar la capa vrt

sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"

Ya con el obtenemos el archivo el cual nos proveerá la capa con los datos de SOC (Soil Organic Carbon) para la extensión del municipio de Boyacá.

# DESCOMENTAR CUANDO LO USE POR PRIMERA VEZ
# LUEGO, COMENTÉLO DE NUEVO
#gdal_translate(paste0(sg_url,datos), file,
#              tr=c(250,250),
#              projwin=bb,
#              projwin_srs =igh)

Vamos a transformar las unidades a las del sistema internacional, de dg/kg a g/kg diviendo por 10

(byaca_soc = terra::rast(file)/10)
## class       : SpatRaster 
## dimensions  : 32, 47, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8188750, -8177000, 601750, 609750  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source(s)   : memory
## varname     : soc_igh_15_30 
## name        : soc_igh_15_30 
## min value   :          29.3 
## max value   :          73.6

4. Exploración de datos

Haremos un histograma de los datos de SOC en el municipio de Boyacá para conocer sus frecuencias

terra::hist(byaca_soc, main= "Soil Organic Carbon", col= "orange")

Con summary() podemos generar un resumen del objeto que contiene los datos

(names(byaca_soc) <-  "soc")
## [1] "soc"
summary(byaca_soc)
##       soc       
##  Min.   :29.30  
##  1st Qu.:40.60  
##  Median :44.10  
##  Mean   :45.09  
##  3rd Qu.:48.90  
##  Max.   :73.60  
##  NA's   :3

Asignemos los valores de SOC en una nueva variable ignorando los NA

valores <- values(byaca_soc, na.rm=TRUE)

Creemos una paleta de colores para un gráfico, lo que nos permitirá reconocer y diferenciar las clases de valores.

orangecyan <- colorNumeric(c("orange","yellow2", "darkseagreen", "cyan"), valores,
  na.color = "transparent")

Y a continuación el gráfico

leaflet::leaflet(byaca) %>% 
  addTiles() %>% 
  setView(-73.35, 5.45, 12) %>%
   addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.5, fillOpacity = 0.10,
    popup = paste("Municipio: ", byaca$MPIO_CNMBR)) %>%
  addRasterImage(byaca_soc, colors ="Spectral", opacity = 0.8)  %>%
  addLegend(pal = orangecyan, values = valores, title = "Soil Organic Carbon (SOC) [g/kg]") 
## Warning: sf layer is not long-lat data
## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'
set.seed(2025)

#Conversión de el sistema de coordenadas de SOC

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(byaca_soc, geog))
## class       : SpatRaster 
## dimensions  : 33, 49, 1  (nrow, ncol, nlyr)
## resolution  : 0.002206672, 0.002206672  (x, y)
## extent      : -73.44269, -73.33456, 5.404657, 5.477477  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :      soc 
## min value   : 31.34987 
## max value   : 72.03896

Muestreo aleatorio de 1000 puntos en el espacio

(samples= terra::spatSample(geog.soc, 1000, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1000, 1  (geometries, attributes)
##  extent      : -73.44158, -73.33566, 5.405761, 5.476374  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      :  47.4
##                42.26
##                40.21
#muestras

(muestras <- sf::st_as_sf(samples))
## Simple feature collection with 1000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.44158 ymin: 5.405761 xmax: -73.33566 ymax: 5.476374
## 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                   geometry
## 1  47.40226 POINT (-73.39304 5.458721)
## 2  42.25752 POINT (-73.35773 5.438861)
## 3  40.21034 POINT (-73.39083 5.441067)
## 4  43.84718 POINT (-73.39966 5.441067)
## 5  42.94074 POINT (-73.42614 5.460927)
## 6  63.28583 POINT (-73.39524 5.469754)
## 7  48.30150 POINT (-73.43276 5.414587)
## 8  47.95037 POINT (-73.33787 5.452101)
## 9  47.21868 POINT (-73.39966 5.445481)
## 10 44.11455 POINT (-73.40628 5.436654)
#Omitamos los NA

nmuestras= na.omit(muestras)
df= data.frame(nmuestras)
View(df)

Ahora visualicemos las muestras eliminadas

longit <- st_coordinates(nmuestras)[,1]
latit <- st_coordinates(nmuestras)[,2]
soc <- nmuestras$soc


#Que obtuvimos???

summary(soc) #Resumen estadístico del objeto
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   31.35   40.84   44.42   45.15   48.62   72.04
length(soc)
## [1] 974

¿Cuántas muestras no contenían datos válidos?

id <- seq(1,974,1)
(sitios <- data.frame(id, longit, latit, soc))
##      id    longit    latit      soc
## 1     1 -73.39304 5.458721 47.40226
## 2     2 -73.35773 5.438861 42.25752
## 3     3 -73.39083 5.441067 40.21034
## 4     4 -73.39966 5.441067 43.84718
## 5     5 -73.42614 5.460927 42.94074
## 6     6 -73.39524 5.469754 63.28583
## 7     7 -73.43276 5.414587 48.30150
## 8     8 -73.33787 5.452101 47.95037
## 9     9 -73.39966 5.445481 47.21868
## 10   10 -73.40628 5.436654 44.11455
## 11   11 -73.42834 5.474167 42.90641
## 12   12 -73.41952 5.452101 39.81597
## 13   13 -73.43717 5.438861 47.80098
## 14   14 -73.36656 5.467547 41.86921
## 15   15 -73.35773 5.432241 43.89245
## 16   16 -73.34890 5.432241 35.50576
## 17   17 -73.43276 5.460927 40.39505
## 18   18 -73.41731 5.412381 49.04708
## 19   19 -73.35111 5.430034 40.08926
## 20   20 -73.39304 5.407967 43.66805
## 21   21 -73.35332 5.456514 39.90910
## 22   22 -73.40848 5.430034 42.01193
## 23   23 -73.39304 5.471961 42.46704
## 24   24 -73.36656 5.460927 36.63972
## 25   25 -73.39304 5.419001 42.55425
## 26   26 -73.41510 5.467547 55.41597
## 27   27 -73.43938 5.423414 59.03062
## 28   28 -73.39083 5.434447 40.27308
## 29   29 -73.40407 5.414587 49.94452
## 30   30 -73.36876 5.438861 32.49722
## 31   31 -73.39524 5.445481 47.78348
## 32   32 -73.34890 5.452101 35.79914
## 33   33 -73.38862 5.412381 44.50109
## 34   34 -73.36656 5.412381 38.93064
## 35   35 -73.38200 5.416794 43.75704
## 36   36 -73.34008 5.476374 43.34033
## 37   37 -73.43717 5.434447 59.36267
## 38   38 -73.37097 5.471961 47.56343
## 39   39 -73.35994 5.460927 38.65915
## 40   40 -73.41510 5.463134 51.09368
## 41   41 -73.39966 5.476374 42.46424
## 42   42 -73.38200 5.463134 55.44702
## 43   43 -73.35552 5.458721 37.49004
## 44   44 -73.37318 5.414587 45.00768
## 45   45 -73.42834 5.436654 43.51689
## 46   46 -73.37097 5.443274 39.24706
## 47   47 -73.43055 5.412381 44.98114
## 48   48 -73.39524 5.458721 47.26334
## 49   49 -73.33787 5.416794 45.79919
## 50   50 -73.38642 5.414587 46.45528
## 51   51 -73.35111 5.443274 38.12983
## 52   52 -73.43055 5.430034 57.62075
## 53   53 -73.43938 5.430034 62.07600
## 54   54 -73.40628 5.447687 49.57507
## 55   55 -73.35111 5.460927 42.46015
## 56   56 -73.40407 5.447687 50.49238
## 57   57 -73.35111 5.416794 47.44403
## 58   58 -73.39966 5.458721 50.26239
## 59   59 -73.42393 5.474167 43.05605
## 60   60 -73.41731 5.421207 47.16410
## 61   61 -73.35552 5.412381 47.01493
## 62   62 -73.41290 5.460927 53.75257
## 63   63 -73.41731 5.456514 48.97870
## 64   64 -73.41069 5.427827 44.33731
## 65   65 -73.41069 5.423414 48.71599
## 66   66 -73.37538 5.419001 49.25046
## 67   67 -73.42614 5.412381 44.66874
## 68   68 -73.41952 5.476374 51.58291
## 69   69 -73.38862 5.449894 43.55551
## 70   70 -73.43276 5.436654 45.05023
## 71   71 -73.39304 5.410174 44.22521
## 72   72 -73.34008 5.441067 34.78304
## 73   73 -73.40848 5.465341 50.71845
## 74   74 -73.35994 5.447687 38.49290
## 75   75 -73.37318 5.465341 41.82494
## 76   76 -73.39966 5.447687 50.21832
## 77   77 -73.40407 5.469754 45.20358
## 78   78 -73.35552 5.421207 38.46742
## 79   79 -73.39083 5.471961 40.53274
## 80   80 -73.42172 5.410174 45.43159
## 81   81 -73.40628 5.423414 50.30973
## 82   82 -73.44158 5.430034 65.44894
## 83   83 -73.40186 5.430034 44.17772
## 84   84 -73.39966 5.419001 45.09396
## 85   85 -73.40628 5.410174 45.92765
## 86   86 -73.42393 5.443274 48.55131
## 87   87 -73.38862 5.467547 61.21534
## 88   88 -73.35332 5.474167 38.67582
## 89   89 -73.36656 5.427827 39.11295
## 90   90 -73.42834 5.469754 40.82159
## 91   91 -73.42172 5.445481 49.34609
## 92   92 -73.39524 5.452101 44.20728
## 93   93 -73.35552 5.463134 38.45116
## 94   94 -73.41290 5.463134 50.04340
## 95   95 -73.42393 5.454307 45.38306
## 96   96 -73.35773 5.423414 36.63765
## 97   97 -73.41952 5.412381 45.29595
## 98   98 -73.38200 5.436654 37.13637
## 99   99 -73.40628 5.463134 47.02560
## 100 100 -73.43055 5.467547 40.53519
## 101 101 -73.37318 5.471961 46.47459
## 102 102 -73.34449 5.467547 37.34631
## 103 103 -73.41731 5.414587 48.99667
## 104 104 -73.33566 5.463134 50.28038
## 105 105 -73.37980 5.419001 46.24824
## 106 106 -73.37980 5.423414 47.98033
## 107 107 -73.38862 5.456514 53.11383
## 108 108 -73.38642 5.438861 36.00539
## 109 109 -73.34670 5.416794 49.78814
## 110 110 -73.33787 5.456514 50.82382
## 111 111 -73.37097 5.452101 40.89188
## 112 112 -73.37538 5.458721 49.42133
## 113 113 -73.41731 5.445481 48.85902
## 114 114 -73.40407 5.454307 50.48294
## 115 115 -73.36214 5.456514 39.72451
## 116 116 -73.34890 5.471961 40.11304
## 117 117 -73.34449 5.469754 39.35167
## 118 118 -73.44158 5.414587 57.29712
## 119 119 -73.43276 5.456514 41.14709
## 120 120 -73.40628 5.405761 46.48852
## 121 121 -73.37980 5.454307 47.23079
## 122 122 -73.41290 5.419001 48.09236
## 123 123 -73.37759 5.436654 39.77884
## 124 124 -73.39304 5.430034 48.75830
## 125 125 -73.37538 5.449894 43.22063
## 126 126 -73.43717 5.474167 45.78143
## 127 127 -73.34449 5.414587 51.53500
## 128 128 -73.33566 5.430034 35.43871
## 129 129 -73.34670 5.456514 36.68510
## 130 130 -73.41510 5.438861 36.95765
## 131 131 -73.43938 5.471961 45.96081
## 132 132 -73.36435 5.421207 42.75901
## 133 133 -73.34890 5.474167 42.70000
## 134 134 -73.33787 5.419001 39.92393
## 135 135 -73.38421 5.454307 56.91447
## 136 136 -73.43717 5.405761 61.91897
## 137 137 -73.44158 5.425621 52.44045
## 138 138 -73.36435 5.465341 49.16417
## 139 139 -73.41069 5.467547 50.84723
## 140 140 -73.39524 5.425621 49.51553
## 141 141 -73.41952 5.467547 51.34970
## 142 142 -73.34670 5.471961 41.32803
## 143 143 -73.36435 5.434447 38.28012
## 144 144 -73.36656 5.421207 42.21025
## 145 145 -73.37538 5.405761 43.39868
## 146 146 -73.34449 5.423414 38.48518
## 147 147 -73.38642 5.465341 60.62045
## 148 148 -73.41069 5.414587 49.27946
## 149 149 -73.41290 5.430034 44.19110
## 150 150 -73.39745 5.471961 41.02042
## 151 151 -73.42834 5.423414 58.35677
## 152 152 -73.42172 5.447687 52.37943
## 153 153 -73.39524 5.414587 44.99604
## 154 154 -73.43276 5.441067 42.77404
## 155 155 -73.36435 5.443274 39.82456
## 156 156 -73.41952 5.445481 53.15687
## 157 157 -73.42172 5.469754 50.33886
## 158 158 -73.37318 5.460927 48.44168
## 159 159 -73.39745 5.416794 45.84288
## 160 160 -73.37759 5.412381 42.10381
## 161 161 -73.37759 5.465341 41.83662
## 162 162 -73.41510 5.434447 42.53495
## 163 163 -73.33566 5.458721 50.30307
## 164 164 -73.33566 5.469754 54.36223
## 165 165 -73.39083 5.419001 41.65295
## 166 166 -73.36656 5.474167 42.88740
## 167 167 -73.42614 5.458721 43.33985
## 168 168 -73.39083 5.443274 42.32495
## 169 169 -73.33787 5.430034 36.46385
## 170 170 -73.37980 5.421207 46.78765
## 171 171 -73.36214 5.414587 43.86512
## 172 172 -73.38642 5.407967 48.13283
## 173 173 -73.41069 5.447687 50.92939
## 174 174 -73.34890 5.410174 50.50776
## 175 175 -73.41510 5.405761 65.80184
## 176 176 -73.37097 5.419001 38.90776
## 177 177 -73.41731 5.458721 51.81182
## 178 178 -73.39966 5.407967 47.64399
## 179 179 -73.43496 5.434447 55.57758
## 180 180 -73.35994 5.467547 40.61130
## 181 181 -73.38200 5.469754 59.64106
## 182 182 -73.34670 5.438861 37.30963
## 183 183 -73.37538 5.412381 42.23064
## 184 184 -73.36656 5.449894 42.46398
## 185 185 -73.35332 5.438861 41.82409
## 186 186 -73.40407 5.471961 37.99984
## 187 187 -73.33787 5.412381 46.40542
## 188 188 -73.36656 5.465341 40.39015
## 189 189 -73.41510 5.471961 43.89641
## 190 190 -73.43496 5.471961 45.93013
## 191 191 -73.35773 5.449894 39.89607
## 192 192 -73.44158 5.419001 55.21532
## 193 193 -73.35773 5.469754 35.92253
## 194 194 -73.36435 5.452101 41.85371
## 195 195 -73.34449 5.412381 53.32421
## 196 196 -73.40407 5.427827 47.43759
## 197 197 -73.37318 5.427827 41.49697
## 198 198 -73.35773 5.443274 44.17723
## 199 199 -73.34449 5.432241 32.52046
## 200 200 -73.41952 5.419001 43.83247
## 201 201 -73.36876 5.405761 44.66120
## 202 202 -73.42614 5.445481 46.51902
## 203 203 -73.41510 5.436654 40.91756
## 204 204 -73.35994 5.405761 48.64407
## 205 205 -73.38421 5.407967 46.96172
## 206 206 -73.42393 5.410174 45.02430
## 207 207 -73.41510 5.476374 38.51451
## 208 208 -73.37097 5.438861 34.07639
## 209 209 -73.37318 5.474167 45.61082
## 210 210 -73.38200 5.434447 38.60172
## 211 211 -73.37759 5.405761 40.60492
## 212 212 -73.36214 5.432241 39.73065
## 213 213 -73.38421 5.425621 46.05936
## 214 214 -73.34890 5.419001 44.23465
## 215 215 -73.35332 5.469754 33.16253
## 216 216 -73.34008 5.425621 40.42813
## 217 217 -73.34228 5.410174 51.62139
## 218 218 -73.39745 5.407967 47.18217
## 219 219 -73.35111 5.447687 36.95216
## 220 220 -73.43938 5.449894 43.58358
## 221 221 -73.37538 5.438861 41.44854
## 222 222 -73.34008 5.474167 41.50867
## 223 223 -73.35773 5.410174 44.72430
## 224 224 -73.43055 5.460927 39.83907
## 225 225 -73.39083 5.438861 37.30058
## 226 226 -73.36435 5.449894 40.38155
## 227 227 -73.36656 5.407967 42.34519
## 228 228 -73.36876 5.414587 40.04138
## 229 229 -73.40848 5.436654 38.77866
## 230 230 -73.36876 5.443274 36.73262
## 231 231 -73.41731 5.460927 57.32101
## 232 232 -73.37097 5.441067 35.90907
## 233 233 -73.37538 5.425621 45.34912
## 234 234 -73.42393 5.421207 48.64487
## 235 235 -73.35332 5.421207 40.55392
## 236 236 -73.41069 5.410174 46.73595
## 237 237 -73.37538 5.454307 40.84047
## 238 238 -73.44158 5.416794 57.17313
## 239 239 -73.41731 5.434447 45.18116
## 240 240 -73.34670 5.474167 43.59105
## 241 241 -73.40628 5.438861 42.17712
## 242 242 -73.35773 5.416794 43.67574
## 243 243 -73.39083 5.425621 45.09306
## 244 244 -73.33787 5.438861 34.33809
## 245 245 -73.41069 5.465341 50.53415
## 246 246 -73.35111 5.414587 51.03080
## 247 247 -73.41290 5.465341 52.11037
## 248 248 -73.39304 5.434447 41.97742
## 249 249 -73.34890 5.423414 46.71787
## 250 250 -73.40628 5.469754 46.60375
## 251 251 -73.44158 5.405761 63.20000
## 252 252 -73.34890 5.436654 37.72066
## 253 253 -73.41952 5.465341 48.99281
## 254 254 -73.40628 5.458721 53.71375
## 255 255 -73.40848 5.458721 55.00873
## 256 256 -73.33787 5.445481 42.39257
## 257 257 -73.40407 5.432241 45.26445
## 258 258 -73.35111 5.423414 42.96600
## 259 259 -73.42393 5.469754 44.92487
## 260 260 -73.44158 5.412381 54.95525
## 261 261 -73.40186 5.441067 45.28738
## 262 262 -73.42834 5.465341 40.80034
## 263 263 -73.40407 5.474167 41.98326
## 264 264 -73.39304 5.454307 46.13153
## 265 265 -73.34228 5.443274 36.68582
## 266 266 -73.40628 5.443274 44.73183
## 267 267 -73.37980 5.456514 53.55440
## 268 268 -73.41510 5.427827 45.61827
## 269 269 -73.40628 5.427827 45.52209
## 270 270 -73.35773 5.436654 40.44595
## 271 271 -73.37097 5.414587 40.09993
## 272 272 -73.41952 5.432241 45.60139
## 273 273 -73.34008 5.449894 50.68363
## 274 274 -73.37318 5.410174 43.60749
## 275 275 -73.37318 5.434447 37.35321
## 276 276 -73.38862 5.407967 48.92815
## 277 277 -73.40186 5.407967 47.06449
## 278 278 -73.38642 5.458721 58.52032
## 279 279 -73.37538 5.407967 42.61964
## 280 280 -73.33787 5.467547 45.61734
## 281 281 -73.43496 5.412381 43.97414
## 282 282 -73.39745 5.441067 44.30735
## 283 283 -73.41069 5.471961 42.22831
## 284 284 -73.40407 5.449894 51.84035
## 285 285 -73.33566 5.471961 45.66602
## 286 286 -73.40848 5.476374 42.29724
## 287 287 -73.39524 5.443274 45.74105
## 288 288 -73.43055 5.438861 42.06468
## 289 289 -73.41290 5.438861 36.56130
## 290 290 -73.37538 5.456514 43.15615
## 291 291 -73.41731 5.407967 55.30029
## 292 292 -73.37318 5.438861 38.24463
## 293 293 -73.42393 5.445481 50.01042
## 294 294 -73.43938 5.427827 54.73120
## 295 295 -73.42393 5.434447 46.09565
## 296 296 -73.40407 5.407967 46.21045
## 297 297 -73.43055 5.407967 44.79149
## 298 298 -73.34890 5.476374 40.69977
## 299 299 -73.42172 5.463134 48.08392
## 300 300 -73.36435 5.463134 43.45834
## 301 301 -73.38421 5.467547 61.92939
## 302 302 -73.33566 5.460927 46.25029
## 303 303 -73.42393 5.465341 47.13096
## 304 304 -73.39524 5.465341 66.40761
## 305 305 -73.38642 5.445481 40.97109
## 306 306 -73.34890 5.405761 47.49922
## 307 307 -73.37318 5.447687 45.40754
## 308 308 -73.37318 5.476374 45.20362
## 309 309 -73.38862 5.438861 37.68571
## 310 310 -73.39966 5.469754 46.62315
## 311 311 -73.42393 5.423414 50.30486
## 312 312 -73.34449 5.463134 33.04419
## 313 313 -73.41731 5.463134 52.04364
## 314 314 -73.35552 5.425621 36.36029
## 315 315 -73.36876 5.407967 42.78098
## 316 316 -73.38200 5.460927 58.17245
## 317 317 -73.34449 5.416794 46.16203
## 318 318 -73.43055 5.405761 42.42271
## 319 319 -73.42834 5.456514 41.08259
## 320 320 -73.39745 5.405761 46.50095
## 321 321 -73.39745 5.463134 58.49112
## 322 322 -73.34449 5.438861 34.20279
## 323 323 -73.43276 5.471961 44.10366
## 324 324 -73.41290 5.458721 52.95012
## 325 325 -73.38862 5.434447 37.96468
## 326 326 -73.33787 5.410174 47.70758
## 327 327 -73.42393 5.414587 43.44710
## 328 328 -73.42834 5.427827 60.02974
## 329 329 -73.38200 5.458721 60.40133
## 330 330 -73.41290 5.474167 44.20406
## 331 331 -73.43717 5.452101 46.63444
## 332 332 -73.37538 5.460927 49.84540
## 333 333 -73.35332 5.463134 35.22337
## 334 334 -73.41290 5.432241 43.35870
## 335 335 -73.40407 5.458721 54.56923
## 336 336 -73.38862 5.452101 46.65857
## 337 337 -73.37538 5.434447 38.82454
## 338 338 -73.36435 5.458721 35.79274
## 339 339 -73.41731 5.427827 45.45077
## 340 340 -73.39745 5.449894 48.68552
## 341 341 -73.39304 5.436654 39.19377
## 342 342 -73.37097 5.449894 41.44303
## 343 343 -73.37759 5.454307 39.16077
## 344 344 -73.43717 5.469754 41.74032
## 345 345 -73.36435 5.469754 44.93378
## 346 346 -73.35552 5.443274 40.09453
## 347 347 -73.43938 5.434447 57.92505
## 348 348 -73.43717 5.443274 46.20481
## 349 349 -73.33787 5.405761 42.55492
## 350 350 -73.35111 5.469754 36.41941
## 351 351 -73.39524 5.421207 45.35098
## 352 352 -73.34670 5.432241 31.34987
## 353 353 -73.34228 5.438861 33.20765
## 354 354 -73.35994 5.436654 39.65306
## 355 355 -73.40407 5.416794 48.15697
## 356 356 -73.38862 5.414587 44.64019
## 357 357 -73.37318 5.407967 44.26708
## 358 358 -73.43938 5.425621 53.80500
## 359 359 -73.36876 5.441067 31.94203
## 360 360 -73.41510 5.412381 47.30489
## 361 361 -73.42614 5.463134 38.94324
## 362 362 -73.41731 5.467547 53.76693
## 363 363 -73.42614 5.465341 40.39139
## 364 364 -73.39966 5.467547 48.04932
## 365 365 -73.39304 5.427827 49.69413
## 366 366 -73.39083 5.465341 63.14013
## 367 367 -73.34890 5.449894 37.86358
## 368 368 -73.37318 5.441067 41.83089
## 369 369 -73.37759 5.458721 52.99266
## 370 370 -73.40848 5.456514 53.04221
## 371 371 -73.39745 5.438861 39.96612
## 372 372 -73.42834 5.430034 59.90351
## 373 373 -73.43938 5.458721 44.68183
## 374 374 -73.42834 5.443274 43.36547
## 375 375 -73.35332 5.414587 48.86217
## 376 376 -73.37097 5.430034 38.54255
## 377 377 -73.41731 5.419001 46.92850
## 378 378 -73.39304 5.456514 50.46093
## 379 379 -73.43717 5.445481 47.18367
## 380 380 -73.39745 5.447687 51.35655
## 381 381 -73.39745 5.430034 48.12897
## 382 382 -73.34228 5.445481 42.90467
## 383 383 -73.43276 5.469754 42.41547
## 384 384 -73.38862 5.419001 43.69330
## 385 385 -73.35332 5.476374 42.84847
## 386 386 -73.39083 5.449894 44.83392
## 387 387 -73.43717 5.449894 44.79073
## 388 388 -73.33787 5.436654 34.31415
## 389 389 -73.35552 5.416794 42.60156
## 390 390 -73.35111 5.432241 40.76215
## 391 391 -73.41069 5.469754 50.97889
## 392 392 -73.41069 5.449894 49.06330
## 393 393 -73.42834 5.405761 49.14398
## 394 394 -73.35332 5.416794 45.91870
## 395 395 -73.34228 5.458721 43.91211
## 396 396 -73.37980 5.445481 41.85647
## 397 397 -73.38421 5.423414 49.58113
## 398 398 -73.42614 5.405761 62.80564
## 399 399 -73.40628 5.414587 49.25579
## 400 400 -73.33566 5.452101 48.12123
## 401 401 -73.43938 5.452101 42.65212
## 402 402 -73.38642 5.412381 47.07903
## 403 403 -73.34228 5.419001 40.75622
## 404 404 -73.37538 5.430034 40.74992
## 405 405 -73.37759 5.449894 41.12054
## 406 406 -73.43496 5.427827 53.89401
## 407 407 -73.41290 5.441067 35.38420
## 408 408 -73.34890 5.443274 39.95263
## 409 409 -73.40186 5.452101 49.95084
## 410 410 -73.39524 5.423414 45.09411
## 411 411 -73.42172 5.434447 44.82993
## 412 412 -73.37980 5.469754 59.01687
## 413 413 -73.41290 5.423414 50.46518
## 414 414 -73.41731 5.454307 45.78771
## 415 415 -73.36214 5.463134 41.77924
## 416 416 -73.42834 5.414587 48.82677
## 417 417 -73.34008 5.454307 48.09661
## 418 418 -73.34890 5.407967 48.71019
## 419 419 -73.40186 5.434447 43.64517
## 420 420 -73.34890 5.456514 40.86304
## 421 421 -73.36876 5.460927 39.02888
## 422 422 -73.43717 5.456514 45.78860
## 423 423 -73.36656 5.471961 41.93270
## 424 424 -73.40628 5.430034 44.97050
## 425 425 -73.42614 5.427827 55.86046
## 426 426 -73.34008 5.416794 46.39150
## 427 427 -73.43496 5.476374 54.86926
## 428 428 -73.36214 5.423414 38.68820
## 429 429 -73.35111 5.454307 43.66240
## 430 430 -73.39745 5.445481 47.74953
## 431 431 -73.41069 5.416794 50.22818
## 432 432 -73.37097 5.474167 48.95340
## 433 433 -73.40848 5.412381 48.50496
## 434 434 -73.43496 5.410174 43.03337
## 435 435 -73.35994 5.458721 39.19220
## 436 436 -73.41731 5.425621 48.87547
## 437 437 -73.36435 5.405761 44.52366
## 438 438 -73.39083 5.467547 63.74224
## 439 439 -73.38200 5.412381 43.45617
## 440 440 -73.41731 5.416794 47.61671
## 441 441 -73.35552 5.476374 40.92491
## 442 442 -73.38862 5.474167 43.83415
## 443 443 -73.37097 5.469754 43.57721
## 444 444 -73.39524 5.416794 45.59480
## 445 445 -73.43496 5.425621 52.13102
## 446 446 -73.37097 5.445481 43.41771
## 447 447 -73.40186 5.447687 50.40054
## 448 448 -73.36435 5.430034 41.26608
## 449 449 -73.39966 5.460927 50.23481
## 450 450 -73.43055 5.445481 45.04114
## 451 451 -73.35773 5.441067 44.14257
## 452 452 -73.38200 5.452101 43.06685
## 453 453 -73.34670 5.463134 31.54183
## 454 454 -73.43055 5.410174 46.54795
## 455 455 -73.35111 5.463134 39.95552
## 456 456 -73.41510 5.430034 45.05549
## 457 457 -73.42172 5.454307 45.33210
## 458 458 -73.42172 5.449894 46.73508
## 459 459 -73.33787 5.463134 48.50465
## 460 460 -73.41069 5.463134 50.90913
## 461 461 -73.37318 5.405761 45.06622
## 462 462 -73.38642 5.467547 59.47989
## 463 463 -73.35994 5.434447 41.68695
## 464 464 -73.36656 5.416794 39.62183
## 465 465 -73.36435 5.410174 37.76670
## 466 466 -73.34670 5.441067 39.02359
## 467 467 -73.40407 5.412381 49.67339
## 468 468 -73.43055 5.423414 57.20042
## 469 469 -73.41510 5.452101 49.57157
## 470 470 -73.40407 5.465341 48.57569
## 471 471 -73.34228 5.474167 42.89139
## 472 472 -73.39083 5.416794 43.50027
## 473 473 -73.39745 5.443274 49.45320
## 474 474 -73.40186 5.425621 50.52450
## 475 475 -73.39304 5.441067 40.84990
## 476 476 -73.39966 5.471961 40.33585
## 477 477 -73.42614 5.469754 40.57530
## 478 478 -73.36656 5.445481 40.04842
## 479 479 -73.34228 5.449894 49.20699
## 480 480 -73.40186 5.465341 49.58255
## 481 481 -73.37759 5.416794 46.46277
## 482 482 -73.39524 5.432241 44.67379
## 483 483 -73.39083 5.430034 46.81084
## 484 484 -73.34670 5.423414 44.79997
## 485 485 -73.34228 5.469754 46.84334
## 486 486 -73.35773 5.447687 38.27044
## 487 487 -73.41069 5.434447 39.94437
## 488 488 -73.37538 5.414587 42.64684
## 489 489 -73.33787 5.460927 48.46933
## 490 490 -73.37097 5.458721 45.80659
## 491 491 -73.38862 5.436654 37.73594
## 492 492 -73.36876 5.458721 36.94270
## 493 493 -73.41731 5.430034 44.34857
## 494 494 -73.40407 5.430034 47.47446
## 495 495 -73.41952 5.471961 42.77400
## 496 496 -73.40407 5.467547 46.83284
## 497 497 -73.34890 5.465341 38.98183
## 498 498 -73.36656 5.432241 40.35704
## 499 499 -73.41731 5.443274 42.14808
## 500 500 -73.33566 5.443274 36.64793
## 501 501 -73.35552 5.414587 45.58502
## 502 502 -73.35332 5.441067 41.88775
## 503 503 -73.39083 5.445481 44.73139
## 504 504 -73.41952 5.436654 46.60120
## 505 505 -73.35332 5.458721 37.34685
## 506 506 -73.39083 5.447687 50.63959
## 507 507 -73.41952 5.405761 63.19221
## 508 508 -73.33566 5.476374 42.81404
## 509 509 -73.34670 5.419001 43.95638
## 510 510 -73.34890 5.438861 40.15751
## 511 511 -73.37538 5.427827 42.43738
## 512 512 -73.34228 5.460927 45.24831
## 513 513 -73.41510 5.414587 51.09039
## 514 514 -73.39304 5.465341 61.15969
## 515 515 -73.42614 5.449894 47.73718
## 516 516 -73.37759 5.421207 46.44115
## 517 517 -73.37097 5.454307 42.00598
## 518 518 -73.39304 5.432241 47.84618
## 519 519 -73.36876 5.467547 40.73384
## 520 520 -73.39966 5.443274 48.16201
## 521 521 -73.34008 5.414587 49.75625
## 522 522 -73.36876 5.421207 40.92021
## 523 523 -73.36214 5.460927 41.01856
## 524 524 -73.43938 5.467547 40.26270
## 525 525 -73.38200 5.414587 42.61871
## 526 526 -73.35773 5.456514 36.89941
## 527 527 -73.40848 5.438861 38.57365
## 528 528 -73.43938 5.476374 51.00000
## 529 529 -73.38421 5.427827 41.66317
## 530 530 -73.38862 5.458721 54.13151
## 531 531 -73.35773 5.414587 44.97250
## 532 532 -73.34890 5.460927 38.23660
## 533 533 -73.43496 5.454307 44.96142
## 534 534 -73.38421 5.469754 59.76800
## 535 535 -73.36214 5.405761 45.94600
## 536 536 -73.41290 5.467547 53.24514
## 537 537 -73.39524 5.405761 45.81394
## 538 538 -73.36656 5.476374 40.39832
## 539 539 -73.43276 5.427827 54.94616
## 540 540 -73.37759 5.474167 48.11385
## 541 541 -73.42614 5.414587 45.87954
## 542 542 -73.35552 5.474167 39.04434
## 543 543 -73.39524 5.430034 50.32883
## 544 544 -73.42834 5.447687 48.94560
## 545 545 -73.37538 5.465341 42.47095
## 546 546 -73.42834 5.445481 46.51774
## 547 547 -73.42172 5.416794 43.27385
## 548 548 -73.37759 5.414587 43.26127
## 549 549 -73.41069 5.452101 50.95281
## 550 550 -73.43717 5.458721 41.59227
## 551 551 -73.38642 5.474167 45.65512
## 552 552 -73.34890 5.427827 37.93399
## 553 553 -73.38200 5.441067 40.50872
## 554 554 -73.42393 5.405761 72.03896
## 555 555 -73.38642 5.449894 40.79681
## 556 556 -73.37538 5.476374 50.35324
## 557 557 -73.39745 5.452101 43.47449
## 558 558 -73.38200 5.419001 45.42834
## 559 559 -73.40848 5.471961 40.92390
## 560 560 -73.43276 5.423414 54.67770
## 561 561 -73.38200 5.443274 39.48820
## 562 562 -73.41510 5.425621 53.12481
## 563 563 -73.41069 5.419001 48.83263
## 564 564 -73.41069 5.441067 37.66321
## 565 565 -73.38421 5.405761 45.32223
## 566 566 -73.43717 5.476374 53.01751
## 567 567 -73.43938 5.460927 41.03379
## 568 568 -73.43717 5.416794 59.58319
## 569 569 -73.40628 5.412381 48.60100
## 570 570 -73.35994 5.476374 41.56091
## 571 571 -73.40186 5.445481 50.20349
## 572 572 -73.34228 5.432241 35.61015
## 573 573 -73.40628 5.445481 47.31370
## 574 574 -73.41069 5.436654 37.39680
## 575 575 -73.37980 5.476374 42.46927
## 576 576 -73.40848 5.449894 51.20391
## 577 577 -73.36214 5.421207 37.49866
## 578 578 -73.42834 5.471961 41.66516
## 579 579 -73.35773 5.463134 38.92830
## 580 580 -73.41510 5.423414 53.22040
## 581 581 -73.38642 5.434447 38.92068
## 582 582 -73.41952 5.474167 42.21921
## 583 583 -73.36435 5.425621 42.29145
## 584 584 -73.41731 5.423414 49.26432
## 585 585 -73.34449 5.427827 35.17061
## 586 586 -73.34890 5.454307 38.39804
## 587 587 -73.41731 5.449894 40.18979
## 588 588 -73.39966 5.452101 45.11057
## 589 589 -73.40628 5.421207 48.53515
## 590 590 -73.42393 5.419001 45.88554
## 591 591 -73.43055 5.454307 44.23693
## 592 592 -73.35994 5.474167 41.81166
## 593 593 -73.41731 5.432241 43.95482
## 594 594 -73.41952 5.438861 44.87574
## 595 595 -73.33566 5.474167 44.14516
## 596 596 -73.35552 5.465341 38.49903
## 597 597 -73.43496 5.467547 43.13591
## 598 598 -73.39966 5.432241 40.63776
## 599 599 -73.43496 5.458721 42.09468
## 600 600 -73.40407 5.463134 47.02023
## 601 601 -73.34670 5.414587 56.42136
## 602 602 -73.43055 5.436654 44.07700
## 603 603 -73.41290 5.405761 64.77224
## 604 604 -73.43276 5.432241 55.39346
## 605 605 -73.36214 5.407967 45.18634
## 606 606 -73.37318 5.443274 42.72950
## 607 607 -73.37759 5.425621 45.42418
## 608 608 -73.33787 5.447687 45.60631
## 609 609 -73.39083 5.469754 58.97887
## 610 610 -73.39966 5.430034 46.28376
## 611 611 -73.43938 5.407967 55.15732
## 612 612 -73.38642 5.423414 48.64305
## 613 613 -73.40407 5.419001 46.31112
## 614 614 -73.37759 5.419001 45.78542
## 615 615 -73.42834 5.410174 45.37014
## 616 616 -73.41510 5.474167 43.37180
## 617 617 -73.43055 5.458721 41.08956
## 618 618 -73.43938 5.443274 47.68595
## 619 619 -73.33566 5.465341 46.67115
## 620 620 -73.41952 5.458721 50.13385
## 621 621 -73.40186 5.469754 45.43066
## 622 622 -73.42393 5.438861 47.00589
## 623 623 -73.39745 5.436654 39.06750
## 624 624 -73.35994 5.463134 39.95602
## 625 625 -73.35994 5.452101 35.59626
## 626 626 -73.35332 5.443274 37.76377
## 627 627 -73.35994 5.430034 39.38533
## 628 628 -73.35773 5.427827 39.85223
## 629 629 -73.41510 5.443274 38.87191
## 630 630 -73.35111 5.419001 43.31575
## 631 631 -73.42393 5.456514 45.13888
## 632 632 -73.33566 5.436654 35.41222
## 633 633 -73.38421 5.419001 44.88336
## 634 634 -73.34670 5.443274 39.01403
## 635 635 -73.42614 5.436654 44.10576
## 636 636 -73.35994 5.449894 37.60000
## 637 637 -73.37759 5.441067 39.45863
## 638 638 -73.34228 5.465341 40.01391
## 639 639 -73.34228 5.441067 35.56310
## 640 640 -73.40407 5.410174 47.62204
## 641 641 -73.43276 5.467547 40.25428
## 642 642 -73.35552 5.445481 35.79804
## 643 643 -73.42172 5.427827 46.46237
## 644 644 -73.37980 5.405761 42.14483
## 645 645 -73.36435 5.476374 39.04353
## 646 646 -73.42172 5.471961 41.96011
## 647 647 -73.42393 5.416794 43.19308
## 648 648 -73.41952 5.456514 47.34280
## 649 649 -73.39524 5.449894 46.74857
## 650 650 -73.43055 5.471961 43.57611
## 651 651 -73.41290 5.447687 45.20364
## 652 652 -73.40628 5.452101 51.06462
## 653 653 -73.39304 5.443274 43.38546
## 654 654 -73.36656 5.436654 33.15445
## 655 655 -73.33787 5.454307 49.95414
## 656 656 -73.39083 5.454307 45.88565
## 657 657 -73.38421 5.430034 42.53691
## 658 658 -73.34449 5.425621 37.49005
## 659 659 -73.43055 5.469754 40.67803
## 660 660 -73.36656 5.463134 39.98840
## 661 661 -73.35552 5.447687 36.97034
## 662 662 -73.43496 5.441067 44.80688
## 663 663 -73.36435 5.414587 37.45103
## 664 664 -73.39083 5.432241 45.12048
## 665 665 -73.42834 5.419001 54.92656
## 666 666 -73.34008 5.421207 40.34056
## 667 667 -73.36876 5.469754 41.51690
## 668 668 -73.37318 5.432241 39.48777
## 669 669 -73.35773 5.452101 39.33462
## 670 670 -73.41510 5.419001 49.31665
## 671 671 -73.40407 5.476374 39.07393
## 672 672 -73.40848 5.467547 49.93373
## 673 673 -73.35994 5.432241 41.34052
## 674 674 -73.38421 5.456514 57.58487
## 675 675 -73.42172 5.423414 49.49359
## 676 676 -73.43055 5.441067 41.10757
## 677 677 -73.35332 5.419001 42.57067
## 678 678 -73.41069 5.432241 42.07642
## 679 679 -73.40186 5.436654 46.61944
## 680 680 -73.38642 5.447687 45.13958
## 681 681 -73.43496 5.460927 41.00821
## 682 682 -73.34449 5.476374 44.21895
## 683 683 -73.36214 5.443274 41.89944
## 684 684 -73.35111 5.407967 50.63838
## 685 685 -73.37318 5.467547 42.96256
## 686 686 -73.33787 5.414587 48.45238
## 687 687 -73.37759 5.460927 49.44168
## 688 688 -73.35773 5.445481 40.06643
## 689 689 -73.35994 5.427827 42.12574
## 690 690 -73.42614 5.441067 42.15796
## 691 691 -73.36214 5.416794 38.73055
## 692 692 -73.37759 5.427827 42.40567
## 693 693 -73.34228 5.436654 33.49640
## 694 694 -73.34008 5.460927 47.20557
## 695 695 -73.34670 5.469754 38.70785
## 696 696 -73.36876 5.476374 41.71130
## 697 697 -73.36214 5.419001 36.18899
## 698 698 -73.40407 5.438861 42.36098
## 699 699 -73.41290 5.416794 50.55104
## 700 700 -73.42834 5.452101 43.93568
## 701 701 -73.42172 5.452101 42.35290
## 702 702 -73.43938 5.447687 43.10969
## 703 703 -73.41510 5.469754 52.05958
## 704 704 -73.33787 5.465341 46.28184
## 705 705 -73.43496 5.405761 52.71515
## 706 706 -73.36214 5.436654 37.80834
## 707 707 -73.42172 5.474167 43.84246
## 708 708 -73.42614 5.425621 55.96672
## 709 709 -73.36656 5.405761 45.73365
## 710 710 -73.39304 5.460927 48.62519
## 711 711 -73.36214 5.476374 39.04472
## 712 712 -73.37097 5.434447 35.76555
## 713 713 -73.43717 5.471961 46.66442
## 714 714 -73.36435 5.407967 40.80494
## 715 715 -73.36656 5.434447 36.43335
## 716 716 -73.37759 5.445481 42.42230
## 717 717 -73.36656 5.419001 41.88402
## 718 718 -73.39304 5.463134 59.95493
## 719 719 -73.33787 5.443274 38.17376
## 720 720 -73.38642 5.425621 44.00524
## 721 721 -73.41069 5.438861 37.33911
## 722 722 -73.34228 5.452101 48.78969
## 723 723 -73.41952 5.430034 46.03648
## 724 724 -73.43276 5.416794 55.09826
## 725 725 -73.43276 5.447687 44.26948
## 726 726 -73.38642 5.441067 39.77383
## 727 727 -73.41290 5.443274 37.53572
## 728 728 -73.42614 5.423414 54.70033
## 729 729 -73.34008 5.407967 50.91446
## 730 730 -73.40186 5.427827 47.07913
## 731 731 -73.43055 5.421207 56.13047
## 732 732 -73.35552 5.430034 42.63944
## 733 733 -73.34890 5.469754 38.47083
## 734 734 -73.41510 5.456514 52.94777
## 735 735 -73.37318 5.421207 47.05555
## 736 736 -73.35773 5.412381 43.64334
## 737 737 -73.41510 5.447687 49.97802
## 738 738 -73.35552 5.423414 35.43890
## 739 739 -73.41069 5.425621 46.83241
## 740 740 -73.36876 5.436654 33.31333
## 741 741 -73.38200 5.427827 41.26371
## 742 742 -73.34670 5.436654 35.86624
## 743 743 -73.44158 5.407967 57.05515
## 744 744 -73.36876 5.434447 32.84466
## 745 745 -73.41290 5.449894 48.75785
## 746 746 -73.43496 5.449894 45.84858
## 747 747 -73.42834 5.425621 58.80851
## 748 748 -73.34890 5.414587 51.82598
## 749 749 -73.36656 5.441067 33.52817
## 750 750 -73.39524 5.474167 42.87407
## 751 751 -73.38200 5.474167 41.88282
## 752 752 -73.38421 5.416794 45.57155
## 753 753 -73.41069 5.458721 54.69836
## 754 754 -73.41069 5.412381 47.27434
## 755 755 -73.39304 5.421207 44.83553
## 756 756 -73.33787 5.441067 34.80117
## 757 757 -73.38200 5.425621 41.90321
## 758 758 -73.39745 5.458721 47.70685
## 759 759 -73.43055 5.427827 58.44082
## 760 760 -73.37759 5.407967 41.32238
## 761 761 -73.38200 5.432241 41.79301
## 762 762 -73.35994 5.445481 43.48042
## 763 763 -73.36876 5.447687 45.21205
## 764 764 -73.36214 5.434447 40.04918
## 765 765 -73.38642 5.405761 47.93385
## 766 766 -73.40848 5.405761 52.57043
## 767 767 -73.34008 5.452101 48.47608
## 768 768 -73.41731 5.476374 45.47673
## 769 769 -73.38421 5.463134 58.41381
## 770 770 -73.41952 5.447687 53.76208
## 771 771 -73.38200 5.407967 44.99913
## 772 772 -73.38642 5.476374 41.00567
## 773 773 -73.38642 5.452101 45.15289
## 774 774 -73.42172 5.458721 49.88632
## 775 775 -73.43276 5.443274 41.85799
## 776 776 -73.36435 5.436654 34.73332
## 777 777 -73.36214 5.410174 45.21509
## 778 778 -73.38421 5.452101 43.91466
## 779 779 -73.36435 5.467547 43.65867
## 780 780 -73.38200 5.449894 42.64629
## 781 781 -73.35552 5.407967 47.57987
## 782 782 -73.41510 5.460927 55.17216
## 783 783 -73.39745 5.414587 48.02464
## 784 784 -73.35111 5.438861 41.12825
## 785 785 -73.34670 5.452101 38.07976
## 786 786 -73.33566 5.441067 34.78856
## 787 787 -73.35994 5.416794 40.34465
## 788 788 -73.35773 5.434447 42.63831
## 789 789 -73.42393 5.425621 50.49945
## 790 790 -73.42614 5.430034 55.08355
## 791 791 -73.34890 5.467547 40.44024
## 792 792 -73.43276 5.410174 43.63637
## 793 793 -73.40186 5.476374 42.21601
## 794 794 -73.35332 5.449894 38.95352
## 795 795 -73.40848 5.414587 49.01724
## 796 796 -73.37318 5.454307 43.77250
## 797 797 -73.39966 5.434447 37.90884
## 798 798 -73.36435 5.456514 40.20930
## 799 799 -73.42834 5.438861 41.28530
## 800 800 -73.39083 5.410174 46.65903
## 801 801 -73.41510 5.421207 48.63790
## 802 802 -73.37097 5.436654 34.56004
## 803 803 -73.35552 5.449894 39.79958
## 804 804 -73.35552 5.436654 40.96066
## 805 805 -73.34008 5.412381 48.96122
## 806 806 -73.36876 5.430034 37.91605
## 807 807 -73.38642 5.410174 48.31697
## 808 808 -73.34008 5.405761 46.51920
## 809 809 -73.36656 5.425621 40.43364
## 810 810 -73.34008 5.471961 39.80046
## 811 811 -73.36656 5.430034 40.67228
## 812 812 -73.42614 5.471961 40.91446
## 813 813 -73.41510 5.407967 54.04846
## 814 814 -73.33787 5.458721 50.06989
## 815 815 -73.37097 5.447687 44.69796
## 816 816 -73.38421 5.414587 44.26747
## 817 817 -73.40848 5.441067 38.20001
## 818 818 -73.34449 5.447687 43.06257
## 819 819 -73.41510 5.445481 51.91584
## 820 820 -73.36435 5.460927 38.86589
## 821 821 -73.37097 5.463134 40.30448
## 822 822 -73.42393 5.447687 53.00161
## 823 823 -73.42393 5.432241 44.55154
## 824 824 -73.39524 5.412381 47.00049
## 825 825 -73.34228 5.471961 41.95512
## 826 826 -73.37759 5.476374 52.77203
## 827 827 -73.41731 5.471961 43.33064
## 828 828 -73.33566 5.445481 39.66276
## 829 829 -73.40848 5.419001 47.97746
## 830 830 -73.35332 5.467547 36.52803
## 831 831 -73.39966 5.474167 42.59058
## 832 832 -73.38200 5.476374 39.14374
## 833 833 -73.39304 5.425621 46.64893
## 834 834 -73.42834 5.449894 45.63355
## 835 835 -73.39304 5.414587 45.22411
## 836 836 -73.36435 5.471961 43.72422
## 837 837 -73.34008 5.465341 39.56614
## 838 838 -73.35994 5.471961 41.58664
## 839 839 -73.39304 5.476374 41.88146
## 840 840 -73.36656 5.456514 39.53939
## 841 841 -73.39524 5.460927 52.10156
## 842 842 -73.36214 5.430034 40.33406
## 843 843 -73.38421 5.476374 41.03047
## 844 844 -73.35111 5.441067 41.77139
## 845 845 -73.35994 5.407967 46.95840
## 846 846 -73.43938 5.441067 47.32513
## 847 847 -73.40628 5.432241 43.39884
## 848 848 -73.41290 5.445481 47.56993
## 849 849 -73.41290 5.425621 49.73243
## 850 850 -73.35994 5.414587 43.90847
## 851 851 -73.43055 5.465341 39.64259
## 852 852 -73.43496 5.474167 45.28473
## 853 853 -73.42614 5.432241 53.15202
## 854 854 -73.36876 5.419001 41.83884
## 855 855 -73.36214 5.474167 39.10812
## 856 856 -73.42834 5.476374 50.47069
## 857 857 -73.34008 5.432241 36.45231
## 858 858 -73.35552 5.456514 39.58107
## 859 859 -73.36876 5.463134 40.26513
## 860 860 -73.41510 5.432241 43.96128
## 861 861 -73.35773 5.425621 37.31580
## 862 862 -73.38862 5.427827 45.23312
## 863 863 -73.41510 5.449894 52.93638
## 864 864 -73.38642 5.460927 59.72892
## 865 865 -73.41952 5.421207 45.94319
## 866 866 -73.43938 5.432241 64.22940
## 867 867 -73.37097 5.432241 36.98088
## 868 868 -73.39966 5.421207 48.39948
## 869 869 -73.37759 5.467547 45.64384
## 870 870 -73.41510 5.458721 54.69753
## 871 871 -73.35773 5.471961 39.69841
## 872 872 -73.35773 5.465341 42.62413
## 873 873 -73.41731 5.405761 63.06640
## 874 874 -73.37318 5.456514 43.29504
## 875 875 -73.41510 5.454307 50.36596
## 876 876 -73.39745 5.454307 43.03474
## 877 877 -73.37538 5.467547 42.13133
## 878 878 -73.40407 5.423414 51.72997
## 879 879 -73.43496 5.423414 53.96479
## 880 880 -73.39083 5.407967 47.09309
## 881 881 -73.36214 5.441067 40.98919
## 882 882 -73.34228 5.416794 46.17708
## 883 883 -73.37318 5.416794 46.73497
## 884 884 -73.42393 5.467547 48.74043
## 885 885 -73.36876 5.432241 35.86665
## 886 886 -73.39524 5.434447 40.60213
## 887 887 -73.39304 5.416794 41.86492
## 888 888 -73.33566 5.456514 51.89601
## 889 889 -73.41952 5.410174 45.52514
## 890 890 -73.34670 5.449894 38.93642
## 891 891 -73.43717 5.432241 63.81195
## 892 892 -73.34008 5.469754 47.81345
## 893 893 -73.36214 5.449894 37.07016
## 894 894 -73.36656 5.414587 38.49582
## 895 895 -73.42172 5.425621 49.51917
## 896 896 -73.41069 5.443274 40.42509
## 897 897 -73.35111 5.436654 39.78691
## 898 898 -73.39304 5.449894 45.72525
## 899 899 -73.42393 5.476374 47.32809
## 900 900 -73.39524 5.471961 42.40784
## 901 901 -73.37980 5.434447 40.23276
## 902 902 -73.41952 5.434447 46.31617
## 903 903 -73.36656 5.458721 36.20028
## 904 904 -73.34228 5.421207 40.18346
## 905 905 -73.34008 5.463134 40.92346
## 906 906 -73.33787 5.476374 42.03540
## 907 907 -73.37318 5.423414 46.58962
## 908 908 -73.38642 5.456514 59.21278
## 909 909 -73.41290 5.414587 49.67503
## 910 910 -73.39745 5.427827 49.77406
## 911 911 -73.37980 5.410174 44.71944
## 912 912 -73.37759 5.432241 42.63743
## 913 913 -73.43717 5.414587 54.50983
## 914 914 -73.37759 5.456514 45.44194
## 915 915 -73.34008 5.467547 41.97190
## 916 916 -73.35332 5.434447 41.70628
## 917 917 -73.41290 5.427827 44.28273
## 918 918 -73.40848 5.407967 48.18471
## 919 919 -73.40186 5.443274 48.03509
## 920 920 -73.39966 5.463134 51.28894
## 921 921 -73.40848 5.410174 46.70950
## 922 922 -73.35111 5.421207 42.61778
## 923 923 -73.35994 5.423414 38.15535
## 924 924 -73.38200 5.430034 43.24614
## 925 925 -73.39524 5.456514 48.19067
## 926 926 -73.40407 5.421207 49.03904
## 927 927 -73.39083 5.436654 38.65682
## 928 928 -73.33566 5.434447 35.32829
## 929 929 -73.41069 5.430034 43.07022
## 930 930 -73.43496 5.419001 57.90462
## 931 931 -73.43717 5.447687 42.90632
## 932 932 -73.42834 5.432241 58.20057
## 933 933 -73.42834 5.421207 57.30176
## 934 934 -73.42172 5.432241 45.13182
## 935 935 -73.37318 5.469754 43.06047
## 936 936 -73.37538 5.421207 47.79208
## 937 937 -73.35111 5.467547 41.87487
## 938 938 -73.40407 5.425621 50.22632
## 939 939 -73.43938 5.469754 41.30481
## 940 940 -73.35332 5.405761 50.69178
## 941 941 -73.37759 5.443274 39.43056
## 942 942 -73.43938 5.445481 46.31640
## 943 943 -73.39745 5.434447 39.38156
## 944 944 -73.37980 5.436654 37.67792
## 945 945 -73.38642 5.432241 43.62373
## 946 946 -73.39524 5.463134 64.33628
## 947 947 -73.43276 5.407967 43.30033
## 948 948 -73.42614 5.467547 41.10035
## 949 949 -73.39966 5.449894 51.43884
## 950 950 -73.37980 5.438861 38.53874
## 951 951 -73.34449 5.430034 33.01776
## 952 952 -73.39966 5.414587 51.76929
## 953 953 -73.42172 5.441067 45.61003
## 954 954 -73.37318 5.449894 43.04182
## 955 955 -73.41290 5.471961 43.13226
## 956 956 -73.35994 5.410174 44.59565
## 957 957 -73.39083 5.456514 49.05302
## 958 958 -73.43055 5.419001 56.26617
## 959 959 -73.42614 5.416794 47.01779
## 960 960 -73.34228 5.414587 50.30021
## 961 961 -73.37980 5.467547 57.76506
## 962 962 -73.36656 5.454307 39.80543
## 963 963 -73.35773 5.405761 49.58494
## 964 964 -73.39745 5.460927 50.47000
## 965 965 -73.39304 5.474167 43.23456
## 966 966 -73.40628 5.425621 47.94738
## 967 967 -73.36876 5.452101 39.72850
## 968 968 -73.35994 5.441067 43.32435
## 969 969 -73.37538 5.441067 41.74892
## 970 970 -73.35552 5.460927 36.87333
## 971 971 -73.40186 5.414587 51.49862
## 972 972 -73.39966 5.427827 47.73553
## 973 973 -73.36435 5.419001 41.03088
## 974 974 -73.42172 5.443274 47.27322

Ahora el mapa de estas muestras

m1 <- leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m1  # Print the map
###Hay que cambiar el sistema de coordenadas de todo, ya que nos están dando en arabia ambos mapas

5. Guardar las muestras

# Verificar y crear directorio si no existe
if (!dir.exists("datos")) {
  dir.create("datos")
}

# Renombrar columnas si es necesario
names(nmuestras) <- substr(names(nmuestras), 1, 10)

# Asegurar que tiene CRS definido
if (is.na(st_crs(nmuestras))) {
  st_crs(nmuestras) <- 3116
}

sf::st_write(nmuestras, "datos/soc_byaca.shp", driver = "ESRI Shapefile", append = FALSE) #Se sobreescribe el archivo ya existente
## Deleting layer `soc_byaca' using driver `ESRI Shapefile'
## Writing layer `soc_byaca' to data source 
##   `datos/soc_byaca.shp' using driver `ESRI Shapefile'
## Writing 974 features with 1 fields and geometry type Point.

6. Interpolación Espacial

La interpolación predice valores para las celdas de un ráster a partir de una cantidad limitada de puntos de datos de muestra. Puede utilizarse para predecir valores desconocidos de cualquier dato de un punto geográfico. Como se mencionó al principio de estas notas, vamos a emplear los métodos de IDW y Kriging Ordinario (OK).

Primero, verifiquemos que los datos están guardados donde los dejamos en la carpeta “datos”

list.files(path="datos", pattern = "*.shp")
## [1] "soc_byaca.shp"

Ahora revisemos que están los valores de las muestras

samples <- sf::st_read("datos/soc_byaca.shp")
## Reading layer `soc_byaca' from data source 
##   `C:\Users\crist\Desktop\Geomática\Taller de Interpolación en Rstudio\datos\soc_byaca.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 974 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.44158 ymin: 5.405761 xmax: -73.33566 ymax: 5.476374
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Recordemmos que nuestra área de estudio que contiene los puntos es el municipio de Boyacá

munic <- sf::st_read("datos/soc_byaca.shp")
## Reading layer `soc_byaca' from data source 
##   `C:\Users\crist\Desktop\Geomática\Taller de Interpolación en Rstudio\datos\soc_byaca.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 974 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.44158 ymin: 5.405761 xmax: -73.33566 ymax: 5.476374
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

6.1 Explorar los datos de entrada (in put)

Obtengamos, antes que todo, un resumen de estadísticas de las muestras (samples)

summary(samples)
##       soc                 geometry  
##  Min.   :31.35   POINT        :974  
##  1st Qu.:40.84   epsg:NA      :  0  
##  Median :44.42   +proj=long...:  0  
##  Mean   :45.15                      
##  3rd Qu.:48.62                      
##  Max.   :72.04

Vamos a generar un histograma para observar la distribución de los datos en términos de frecuencia, cada barra del histograma representa la cantidad de observaciones que caen dentro de un rango específico de valores.

hist(samples$soc, main= "SOC en suelos de Boyacá", col= "skyblue", ylab= "Frecuencia", xlab= "Muestras (g/kg)", breaks = 12)

Vamos a redondear a dos dígitos tales valores muestreados

samples$soc= round(samples$soc, 2) #Redondeo a dos cifras decimales

Conozcamos la distribución espacial de los datos de muestras

pal= colorNumeric(c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$soc,
  )

Una visualización

# This a simple visualization
# Change the code to suit your preferences
leaflet() %>%
  addCircleMarkers(
    data = munic, #Los datos son únicamente de polígonos
    color = "gray",
    # set the opacity of the outline
    opacity = 1,
    # set the stroke width in pixels
    weight = 1,
    # set the fill opacity
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~pal(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "SOC:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

6.2 Proceso de Interpolación

Creación del Objeto tipo gstat

Primero, antes de interpolar, necesitamos un objeto de tipo gstat usando la función gstat(). Este contiene toda la información necesaria para conducir y llevar a cabo una interpolación espacial, es decir, gstat entiende cuál tipo y modelo de interpolación queremos usar: Kriging Ordinario o IDW.

Los tres parámetros que vamos a usar para la función son:

  • Fórmula: la fórmula de predicción especificando la variable dependiente e independiente (aka covariantes).

  • Datos: los datos de calibración (aka datos de tren)

  • Modelo: el modelo del variograma

6.3 Interpolación con IDW (Inverse Distance Weighted)

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

g1= gstat(formula= soc~1, data= samples)

Ya definido, tenemos que especificar la función de predicción para interpolar, o estimar los valores SOC en puntos donde no fueron muestreados.

La función de predicción acepta:

  • Un ráster, objetos stars, tales como DEM
  • Un modelo, objetos gstat, tales como g1

Los ráster tienen dos propósitos principales:

  • Especificar las localizaciones donde queremos hacer las predicciones
  • Especificar los valores covariantes (in Universal Kriging únicamente)

Creemos un ráster con los valores de celda igual a 1. La extensión de este ráster debe cubrir toda nuestra área de interés.

Nota: podemos usar los metadatos de cualquier ráster que cubra nuestra área para definir para definir los parámetros de creación del ráster:

(rrr = terra::rast(xmin=-73.44158, xmax=-73.33566, ymin=5.405761, ymax=5.476374, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326")) ##xmin: -73.44158 ymin: 5.405761 xmax: -73.33566 ymax: 5.476374
## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.0003219453, 0.0001908459  (x, y)
## extent      : -73.44158, -73.33566, 5.405761, 5.476374  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Ahora tenemos que convertir este SpatRaster en un objeto stars

stars.rrr <- stars::st_as_stars(rrr)

La siguiente expresión interpola los valores SOC de acuerdo al modelo g1 y al ráster definido en stars.rrr

Salidainterpolación por IDW

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

Ya con el modelo de interpolación corrido, podemos observar que posee dos atributos, el que nos interesa es el primero ya que el segundo no posee valores, y se rellena con “NA” ya que no está teniendo en cuenta propiedades estadísticas como la varianza. Podemos crear un subconjunto del primer atributo y renombrarlo:

# dealing with NA values
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0

Vamos a renombrar el atributo que nos interesa

#Renombrar
names(z1)= "soc"

Para visualizar la superficie interpolada vamos a crear una paleta de colores

# Paleta de colores, cambielos a su preferencia
paleta <- colorNumeric(palette = c("red", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")

Ahora visualicemos la salida de la interpolación

m2 <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("red", "yellow", "cyan", "green"), 
                                  domain = 10:100)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m2  # Print the map

Se puede observar las muestras obtenidas y los puntos interpolados, obteniendo una superficie de interpolación como si todo el área hubiera sido muestreada, claramente basados en el vecino más cercano, entre más cerca estén dos puntos, mayores similitudes van a poseer.

6.4 Interpolación por Kriging Ordinario

Requiere un modelo del variograma. Es una forma de cuantificar los patrones de autocorrelación en los datos, y asigna ponderaciones en consecuencia de las predicciones. Como primer paso, tenemos que calcular y examinar el variograma empirico usando la funcion variogram().

Variogram() requiere dos argumentos:

  • fórmula: específique la variable dependiente y la covariable, tal cual en gstat.
  • Datos: la capa de puntos con la variable dependiente y covariables como atributos de puntos.
### Expresión que calcula el variograma empíico de **samples**, sin covariables:

v_emp_ok = variogram(soc ~ 1, data=samples)

Ahora grafiquemos el variograma

plot(v_emp_ok) #Grafico del Variograma empírico

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

La función escoge el mejor tipo de ajuste de modelo, y también puede afinar sus parámetros. POdemos usar show.vgms() para visualizar varios tipos de modelos de variograma.

Tenga en cuenta que autofitVariogram() no trabaja en sf_objects, lo cual es porque convertimos samples en un dataframe de puntos espaciales.

El modelo ajustado puede ser graficado con plot()

### Modelo de OK

v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

A continuación el gráfico

plot(v_mod_ok)

Conforme aumenta la distancia aumenta la varianza de los datos, hasta un punto en el que no cambia debido a que ya se alejó demasiado como para compararlos y relacionarlos.

El objeto de salida es realmente una lista con varios componentes, incluyendo el variograma empírico y el modelo ajustado del variograma. El $var_modelcomponent del objeto de salida contiene el modelo actual:

v_mod_ok$var_model
##   model    psill    range kappa
## 1   Nug  0.00000 0.000000   0.0
## 2   Ste 35.10786 1.472366   0.7

Ahora el modelo del variograma puede ser pasado a la función gstat, y podemos continuar con la interpolación Kriging Ordinario.

Salida Interpolación por OK

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

Haremos lo mismo, crearemos un subconjunto con los valores de los atributos predecidos y renombrarlos.

a2 <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"

Los valores interpolados son mostrados a continuación

m3 <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 10:100)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "OK SOC interpolation [%]"
    )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m3 #Imprimir el gráfico

6.5 Evaluación de Resultados

Evaluación cualitativa

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")

A continuación, un gráfico que muestra un ráster con los controles de los modelos

m4 <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Soil organic carbon [%]"
)
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored

## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m4 #Imprimir el Gráfico

Según la salida, los rásters muestran diferencias claras en como interpolan los valores, uno por un modelo de variograma y el otro por el vecino más próximo.

7. Validación Cruzada

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

En Validación Cruzada Leave-One-Out nosotros:

  • Sacamos un punto de los datos de calibración
  • Hacemos una predicción para ese punto
  • Repetimos para todos los puntos en el final, lo que obtenemos es una tabla con un valor observado y un valor predecido para todos los puntos.

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

Cuando escribamos el siguiente chunk, esconda mensajes y resultados:

## Oculte este código
# Puede ser ejecutado desde la consola, no desde aquí
#cv2= gstat.cv(g2) #NO CORRER DESDE LA FUENTE

El resultado posee los siguientes

  • var1.pred - valor predecido
  • VAR1.VAR - Varianza (únicamente para Kriging)
  • Observed - Valor Observado
  • Residual - Observado - Predecido
  • zscore - Z-score (únicamente para Kriging)
  • fold - Validación Cruzada ID

Ahora convirtamos el objeto cv2 en st_as_sf el cual es una superficie espacial compatible con el paquete sf (Simple Features)

## Igualmente correr en la consola, no desde la fuente

# cv2= st_as_sf(cv2)
# Simple Features

Ahora grafiquemos los valores residuales en un gráfico de burbujas

# run from the console
# sp::bubble(as(cv2[, "OK residual"], "Spatial"))

Ahora vamos a calcular el valor del RMSE para la interpolación por OK.

# This is the RMSE value for the OK interpolation
# run from the console
#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2)) #OK interpolation

El valor del RMSE es de 1.914539, no es un valor muy alto, podemos decir que es exacto.

Ahora vamos a correr el código para el IDW

8. RMSE para IDW

#Correr desde la consola
##cv1= gstat.cv(g1)

Lo transformamos en sf

#cv1= st_as_sf(cv1)

Gráfica de Burbujas de los residuales

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

Calculamos y reportamos el RMSE para el método IDW

#sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1)) ##IDW

El calculo nos arrojó 3.653438 para el IDW, lo que nos quiere decir que el método de Kriging Ordinario es mucho más exacto para la interpolación espacial que el IDW, dada su facilidad el resultado era de esperarse, ya que utiliza un método mediado por estadísticas y modelo como el Variograma.

Guardar las Salidas en archivos ráster

write_stars(
  z1, dsn = "datos/IDW_soc_byaca.tif", layer = 1
) 
write_stars(
  z2, dsn = "datos/OK_soc_byaca.tif", layer = 1
)

Conclusiones

Aprendimos teoría y como utilizar en R los métodos de interpolación espacial determinística y estadística determinística por medio de datos de Carbón Orgánico en Suelos provenientes de SoilGrids basados en un shapefile del municipio de Boyacá, en el departamento de Boyacá (Colombia) y tuvimos la oportunidad de compararlos por medio de un Error Cuadrático Medio (RMSE) donde se observó que el método más confiable en nuestro ejemplo patra este tipo de datos es el Kriging Ordinario, con un RMSE fundamentalmente menor que el de IDW (Inverse Distance Weighted).