Natalia Martínez Z Agosto de 2018

Nota: Para la elaboración de este cuaderno se utilizó el tutorial Advanced Techniques With Raster DAta: Part3 - Excercises, realizado por Joao Goncalves publicado el 28 de Febrero de 2018.

Introducción

Las imagenes digitales contienen gran cantidad de datos que gracias al procesamiento digital de imagenes se obtiene informacion sobre los territorios. La extracción de información temática, hace parte del procesamiento de imagenes, donde se busca mapear o estimar las categorias presentes en la imagen. Para esto se utilizan métodos cualitativos como los son la clasificacion supervisada y no supervisada que permite establecer categorias biofisicas; y métodos cualitativos como son las regresiones las cuales estiman propiedades biofísicas.

La regresión Kriging se considerada una técnica de interpolación hiblida, ya que combina interpolacion basada en observaciones puntales de la variable objetivo e interpolacion basada en la regresion de la variable objetivo sobre información auxiliar.

Este cuaderno presenta una aplicación de la regresión kriging, utilizando una imagen Landsat 8 perteneciente a los llanos orientales, con la cual se busca obtener resultados de una clasificación temática donde se etimen los contenidos de carbono organico en el suelo.

Regresión Kriging

La regresion kriging es una técnica de interpolación espacial que combina una regresion lineal o multiple de la variable dependiente en predictores y de los kriging de los residuos de predicción (Hengl et al, 2007).

La ecuacion de regresion predice el valor de la variable dependiente utilizando una funcion lineal en las variables independientes (FAO, 2018).

“Kriging se basa en modelos estadisticos que incluyen una autocorrelación, es decir la relacion estadística entre los puntos medidos. La distancia o direccion entre los puntos de muestra, en kriging, reflejan una correlación espacial que se puede utilizar para explicar la variacion en la superficie” (https://pro.arcgis.com/es/pro-app/tool-reference/3d-analyst/how-kriging-works.htm).

Librerias

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-7
library(ggplot2)
library(raster)
## Warning: package 'raster' was built under R version 3.4.4
library(sp)
library(landsat8)
## Warning: package 'landsat8' was built under R version 3.4.4
library(corrplot)
## Warning: package 'corrplot' was built under R version 3.4.4
## corrplot 0.84 loaded
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.4.4
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(mgcv)
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## This is mgcv 1.8-22. For overview type 'help("mgcv-package")'.
library(gstat)
## Warning: package 'gstat' was built under R version 3.4.4
Datos

Los datos utilizados para la elaboración de este cuaderno fueron datos de carbono organico en suelo tomados en el centro de investigación la libertad tomados en el año 2008.

setwd("C:/Maestria en Geomatica/Percepcion remota avanzada/Taller 1/Landsat_8")
valores <- readOGR("valores5.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "valores5.shp", layer: "valores5"
## with 140 features
## It has 11 fields
## Integer64 fields read as strings:  OBJECTID_1 OBJECTID X Y COBERTURA
lalibertad <- readOGR("C:/Maestria en Geomatica/Percepcion remota avanzada/Taller 1/Landsat_8/Libertadgeo.shp")
## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv
## = use_iconv, : ogrInfo: C:/Maestria en Geomatica/Percepcion remota
## avanzada/Taller 1/Landsat_8/Libertadgeo.dbf not found
## OGR data source with driver: ESRI Shapefile 
## Source: "C:/Maestria en Geomatica/Percepcion remota avanzada/Taller 1/Landsat_8/Libertadgeo.shp", layer: "Libertadgeo"
## with 267 features
## It has 0 fields
COS <- read.csv("C:/Maestria en Geomatica/Percepcion remota avanzada/Taller 1/Landsat_8/COS.csv")
banda1 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B1.TIF")
banda2 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B2.TIF")
banda3 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B3.TIF")
banda4 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B4.TIF")
banda5 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B5.TIF")
banda6 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B6.TIF")
banda7 <- raster("LC08_L1TP_007057_20170814_20170825_01_T1_B7.TIF")

CORRECCIÓN ATMOSFÉRICA

Radiancia

banda3.dn <- as(banda3, 'SpatialGridDataFrame')
banda3.rad <- radconv(banda3.dn, 1.1548E-02, -57.74164)
banda4.dn <- as(banda4, 'SpatialGridDataFrame')
banda4.rad <- radconv(banda4.dn, 9.7382E-03, -48.69100)
banda5.dn <- as(banda5, 'SpatialGridDataFrame')
banda5.rad <- radconv(banda5.dn, 5.9593E-03, -29.79646)
banda6.dn <- as(banda6, 'SpatialGridDataFrame')
banda6.rad <- radconv(banda6.dn, 1.4820E-03, -7.41011)

Reflectancia

refb3 <- reflconvS(banda3.dn, 2.0000E-05, -0.100000, 61.24856156)
refb4 <- reflconvS(banda4.dn, 2.0000E-05, -0.100000, 61.24856156)
refb5 <- reflconvS(banda5.dn, 2.0000E-05, -0.100000, 61.24856156)
refb6 <- reflconvS(banda6.dn, 2.0000E-05, -0.100000, 61.24856156)
B3 <- as(refb3, 'RasterLayer')
B4 <- as(refb4, 'RasterLayer')
B5 <- as(refb5, 'RasterLayer')
B6 <- as(refb6, 'RasterLayer')

INDICES

Indice de vegetación ajustado al suelo SAVI

Dado por la relacion entre las bandas NIR y red de la siguiente manera: [(NIR - Red)/(NIR + Red + L)] * (1 + L)

Indice de vegetación normalizada NDVI

Dado por la relacion [(NIR - Red)/(NIR + Red)]

Indice diferencial de agua normalizado NDWI_1

En el método Gao 1996 la relacion multiespectral es entre las bandas NIR (infroarrojo cercano) y SWIR asi: [(NIR - SWIR)/(NIR + SWIR)]

Indice diferencial de agua normalizado NDWI

El método McFeeters, 1996 relaciona las bandas NIR y Green de la siguiente manera: [(Green - NIR)/(Green + NIR)]

savi <- overlay(B5, B4, fun=function(x,y){((x-y)/(x+y+1))*(2)})
ndvi <- overlay(B5, B4, fun=function(x,y){((x-y)/(x+y))})
ndwi_1 <- overlay(B5, B6, fun=function(x,y){((x-y)/(x+y))})
ndwi <- overlay(B3, B5, fun=function(x,y){((x-y)/(x+y))})

RECORTE POR ZONA DE ESTUDIO LA LIBERTAD y REPROYECCION

crs(banda5)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
crs(lalibertad)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
lalibertad <- spTransform(lalibertad, CRS = crs(banda5))
cortar savi, ndvi, ndwi con la zona de estudio de la libertad
savi <- crop(savi, extent(lalibertad))
savi <- mask(savi, lalibertad)
ndvi <- crop(ndvi, extent(lalibertad))
ndvi <- mask(ndvi, lalibertad)
ndwi_1<- crop(ndwi_1, extent(lalibertad))
ndwi_1 <- mask(ndwi_1, lalibertad)
ndwi <- crop(ndwi, extent(lalibertad))
ndwi <- mask(ndwi, lalibertad)
plot(savi, xlab="INDICE SAVI")

plot(ndvi, xlab="INDICE NDVI")

plot(ndwi_1, xlab="INDICE NDWI_1")

plot(ndwi, xlab="INDICE NDWI")

knitr::kable(head(COS, n=146))
FID X Y COS10 Altitud SAVI NDVI NDWI_1 NDWI COBERTURA
0 669030 448722 20.925 334 0.353243 0.587312 0.137751 -0.553431 231
1 669330 448734 24.894 330 0.393609 0.657186 0.261281 -0.604291 231
2 670292 448697 21.712 335 0.409268 0.689646 0.206797 -0.606751 231
3 669330 449102 14.399 327 0.418529 0.741943 0.380933 -0.642464 231
4 669547 448912 12.127 325 0.452179 0.800162 0.413086 -0.684220 231
5 669709 448952 13.835 324 0.356344 0.657743 0.211107 -0.586219 231
6 671321 448700 13.343 330 0.353791 0.704545 0.209911 -0.579770 231
7 670881 448699 18.370 334 0.436663 0.702336 0.276042 -0.611923 231
8 670548 448700 22.340 337 0.394027 0.736339 0.280111 -0.617732 231
9 671340 448969 17.785 331 0.375891 0.696917 0.238469 -0.601935 231
10 671082 448966 21.150 330 0.363304 0.675810 0.173271 -0.585470 231
11 670862 448968 18.005 332 0.413202 0.762726 0.333302 -0.638345 231
12 672432 449002 20.304 325 0.254349 0.520124 0.163089 -0.492877 231
13 672162 448976 15.640 329 0.317503 0.518371 -0.010307 -0.514196 231
14 671899 448965 15.904 327 0.445024 0.709578 0.244262 -0.591312 231
15 671611 448967 18.353 328 0.403133 0.685907 0.248822 -0.587568 231
16 670561 448982 15.088 332 0.355428 0.586971 0.110262 -0.545278 231
17 669213 449523 24.532 322 0.407725 0.691820 0.292292 -0.599574 231
18 669838 449510 32.055 322 0.368960 0.657006 0.203906 -0.582404 231
19 670021 449522 18.819 321 0.476896 0.767343 0.409524 -0.663990 231
20 670292 449525 23.808 319 0.449450 0.719974 0.299453 -0.625634 231
21 670613 449523 21.071 318 0.354067 0.581747 0.153978 -0.547255 231
22 670831 449522 18.597 317 0.393127 0.609508 0.224889 -0.557771 231
23 670831 449792 21.335 318 0.376794 0.644467 0.247351 -0.572206 231
24 670827 449248 18.478 319 0.369568 0.662246 0.225181 -0.589265 231
25 671101 449252 19.187 323 0.376040 0.616577 0.171704 -0.562406 231
26 671371 449252 20.536 325 0.450477 0.730380 0.307737 -0.647119 231
27 671641 449252 21.876 323 0.395321 0.632739 0.178372 -0.572923 231
28 671911 449252 18.312 326 0.432446 0.791467 0.391374 -0.670543 314
29 672170 449298 16.752 327 0.460050 0.799672 0.305734 -0.675221 231
30 672721 451142 24.838 319 0.389464 0.686564 0.320402 -0.600341 231
31 672451 450872 31.979 311 0.402370 0.682287 0.252101 -0.601329 231
32 672181 450872 23.565 313 0.424125 0.722491 0.430200 -0.634541 231
33 672181 450602 42.775 312 0.354797 0.629910 0.237792 -0.566882 231
34 672181 450332 29.642 312 0.498781 0.772914 0.380233 -0.668732 231
35 670561 449252 18.844 317 0.425549 0.792063 0.402859 -0.665526 314
36 671891 449798 22.361 314 0.403277 0.668197 0.253331 -0.607087 231
37 671911 450062 28.073 312 0.439456 0.667830 0.258292 -0.599750 231
38 671614 449794 28.278 317 0.405457 0.748296 0.358618 -0.640229 231
39 671371 449792 25.113 318 0.459357 0.729783 0.273924 -0.630775 231
40 671641 450062 21.875 313 0.465677 0.739623 0.323462 -0.630254 231
41 671364 450079 28.411 314 0.346487 0.629794 0.251038 -0.547580 231
42 671095 450076 16.020 321 0.351663 0.629198 0.241896 -0.571675 231
43 670838 450012 20.806 321 0.410554 0.737982 0.312735 -0.626602 231
44 670561 449792 24.728 318 0.378887 0.656744 0.218727 -0.574914 231
45 670021 449252 17.646 322 0.453824 0.748153 0.362745 -0.641805 231
46 669737 449210 14.532 325 0.377855 0.672112 0.280779 -0.581167 231
47 669482 449250 16.451 324 0.348069 0.629493 0.178243 -0.553295 231
48 669223 449236 24.577 326 0.336629 0.613292 0.125226 -0.571735 231
49 670559 450872 21.581 318 0.429160 0.723798 0.352389 -0.626935 231
50 670291 450602 19.579 320 0.388234 0.673350 0.309533 -0.597047 231
51 670014 450329 23.828 325 0.405008 0.674429 0.342941 -0.602222 231
52 669786 448532 23.685 334 0.424088 0.700284 0.285391 -0.615764 231
53 669502 448517 15.915 335 0.373145 0.633463 0.185423 -0.563206 231
54 669373 448526 19.882 335 0.454919 0.704512 0.250211 -0.610788 231
55 669298 448687 15.231 330 0.453642 0.704109 0.256726 -0.617909 231
56 668910 448520 17.516 334 0.347707 0.595135 0.175014 -0.538084 231
57 670682 448659 14.292 337 0.342129 0.598905 0.042715 -0.546104 231
58 670694 448961 16.942 333 0.465684 0.746340 0.270153 -0.632660 231
59 672136 449657 16.414 316 0.435489 0.798941 0.302071 -0.673480 231
60 672646 449669 16.371 317 0.461163 0.786653 0.378555 -0.648817 231
61 672595 448835 24.906 324 0.378764 0.696822 0.288885 -0.604171 231
62 671745 453293 18.353 319 0.508426 0.805631 0.363755 -0.688193 2232
63 672183 453305 14.800 316 0.456695 0.763116 0.294044 -0.656337 2232
64 672463 453303 21.588 312 0.397630 0.705189 0.226872 -0.628297 2232
65 672713 453296 17.419 314 0.449008 0.791792 0.411049 -0.676374 2232
66 672437 453025 14.827 320 0.544047 0.803458 0.400935 -0.698267 2232
67 672179 453034 26.020 315 0.476489 0.767143 0.315811 -0.654530 2232
68 671903 453028 13.010 322 0.472119 0.766255 0.327276 -0.672276 2232
69 671637 453033 13.418 321 0.569678 0.816865 0.414919 -0.702497 2232
70 671649 452756 21.804 322 0.533246 0.800548 0.374412 -0.685372 2232
71 671923 452767 22.021 319 0.584918 0.826721 0.387240 -0.717515 2232
72 672184 452752 20.159 319 0.565424 0.826236 0.417725 -0.708466 2232
73 672451 452762 16.837 320 0.572650 0.821765 0.426653 -0.713274 2232
74 672451 452492 19.138 321 0.594096 0.831541 0.444153 -0.720519 2232
75 672181 452492 19.233 322 0.547533 0.806140 0.416928 -0.698862 2232
76 671911 452492 18.332 321 0.543908 0.804512 0.397474 -0.695136 2232
77 671641 452492 20.162 322 0.570154 0.826781 0.415491 -0.711106 2232
78 670561 450602 25.950 317 0.352344 0.638246 0.228023 -0.562202 231
79 670609 450449 15.446 319 0.388156 0.702204 0.283610 -0.603032 314
80 670291 450332 13.805 320 0.352854 0.627457 0.244304 -0.561271 231
81 670291 450062 14.042 319 0.408659 0.697363 0.337265 -0.605301 231
82 670021 450062 26.048 321 0.397053 0.674414 0.233931 -0.591426 231
83 669758 450059 48.195 322 0.360725 0.660873 0.278862 -0.566421 231
84 669481 450062 19.262 326 0.379822 0.652174 0.185277 -0.581681 231
85 672451 451682 26.981 313 0.376993 0.629207 0.129827 -0.562175 231
86 672394 451403 26.110 313 0.445624 0.774751 0.378000 -0.652182 231
87 672724 451375 28.764 311 0.447395 0.737286 0.294692 -0.645654 231
88 672753 451514 17.842 309 0.447611 0.787767 0.337462 -0.676659 231
89 672893 451459 21.083 308 0.463231 0.788342 0.389545 -0.649797 231
90 672721 451952 23.972 314 0.407886 0.648962 0.171979 -0.575677 231
91 672713 452223 21.993 311 0.520969 0.806590 0.366099 -0.705918 231
92 672275 452237 57.849 313 0.425405 0.770257 0.349993 -0.647554 231
93 671911 452222 17.492 314 0.485422 0.798934 0.402944 -0.694154 231
94 671641 452222 18.569 315 0.474679 0.805240 0.422635 -0.693707 231
95 671371 452222 18.948 316 0.477962 0.802323 0.418186 -0.696342 314
96 672241 451267 17.623 320 0.441143 0.794939 0.403698 -0.666606 314
97 672429 451159 25.914 317 0.475987 0.806374 0.432348 -0.671142 314
98 672451 450600 21.268 310 0.445691 0.742093 0.397314 -0.646479 231
99 672757 450285 32.839 311 0.430347 0.775718 0.341715 -0.628525 314
100 672478 450266 35.191 309 0.413162 0.750496 0.310078 -0.643489 314
101 672695 449780 8.466 312 0.473870 0.808989 0.410946 -0.688640 231
102 672398 449762 13.075 317 0.477265 0.797904 0.379610 -0.676465 231
103 672452 450068 19.682 313 0.406116 0.659622 0.187899 -0.582510 231
104 672761 450050 29.078 312 0.497715 0.794082 0.357743 -0.673982 231
105 672162 450020 26.505 313 0.454715 0.792056 0.408067 -0.666412 314
106 669752 449794 26.877 322 0.360423 0.635578 0.195349 -0.570019 231
107 670018 449792 19.238 321 0.418271 0.716166 0.335140 -0.621342 231
108 670292 449795 21.179 320 0.402942 0.701790 0.309161 -0.598456 231
109 670565 450066 19.432 321 0.389195 0.661449 0.292765 -0.587694 231
110 670582 450379 18.482 319 0.403971 0.761420 0.365364 -0.626681 231
111 670833 450330 21.318 324 0.362620 0.632176 0.184423 -0.568000 231
112 671068 450295 31.290 324 0.419297 0.769746 0.393712 -0.653078 231
113 671371 450602 14.098 318 0.457289 0.771150 0.385117 -0.672791 231
114 671640 450601 21.080 321 0.415368 0.744064 0.370209 -0.641725 231
115 671375 451968 11.702 323 0.446105 0.794078 0.382849 -0.671517 314
116 671572 452064 21.341 322 0.486782 0.811024 0.418333 -0.694289 314
117 671604 452032 15.778 321 0.480304 0.790856 0.366463 -0.681825 314
118 671738 451851 19.407 318 0.479059 0.781910 0.353156 -0.678568 231
119 672122 451430 14.326 312 0.447555 0.793564 0.395931 -0.680730 314
120 672061 451300 31.427 315 0.433772 0.788139 0.381808 -0.664860 314
121 671992 451371 17.485 311 0.471863 0.776871 0.328952 -0.692034 231
122 671757 451436 18.598 314 0.428231 0.714270 0.251946 -0.639311 231
123 671604 451444 16.005 315 0.447868 0.763677 0.324633 -0.673329 231
124 671632 451295 16.653 317 0.450181 0.730547 0.305484 -0.645854 231
125 671400 451164 12.972 316 0.491871 0.807638 0.425098 -0.693968 231
126 671002 451307 17.604 318 0.464433 0.792417 0.390217 -0.672653 231
127 671912 450331 20.608 315 0.374940 0.676222 0.294091 -0.589333 231
128 672013 450625 45.979 312 0.394533 0.723324 0.297129 -0.637026 231
129 671639 450332 12.672 313 0.404392 0.677115 0.327227 -0.594113 231
130 671369 450332 23.660 317 0.383869 0.676672 0.338331 -0.590800 231
131 671066 450591 36.170 319 0.466748 0.762636 0.358941 -0.665118 231
132 670869 450823 16.645 321 0.336832 0.549770 0.061304 -0.517071 231
133 670648 451050 20.737 316 0.463722 0.790006 0.384942 -0.676422 231
134 669752 450335 24.577 328 0.403004 0.715316 0.309239 -0.605015 231
135 669479 450333 15.214 327 0.451672 0.792795 0.389066 -0.676172 231
136 669208 450061 18.039 329 0.385911 0.648637 0.228996 -0.581590 231
137 669210 449790 20.365 324 0.384374 0.639783 0.193134 -0.577914 231
138 669479 449789 23.731 325 0.425058 0.702916 0.265678 -0.619234 231
139 669956 449030 42.055 323 0.469126 0.782861 0.330706 -0.662040 231

MATRIZ DE CORRELACION

library(corrplot)
corMat <- cor(COS)
corMat <- cor(COS[,3:ncol(COS)])
corMat
##                     Y       COS10     Altitud       SAVI        NDVI
## Y          1.00000000  0.02584042 -0.55945329  0.6121277  0.54000547
## COS10      0.02584042  1.00000000 -0.26837236 -0.1149885 -0.05558253
## Altitud   -0.55945329 -0.26837236  1.00000000 -0.2614483 -0.32497181
## SAVI       0.61212770 -0.11498851 -0.26144826  1.0000000  0.86730024
## NDVI       0.54000547 -0.05558253 -0.32497181  0.8673002  1.00000000
## NDWI_1     0.50263155 -0.02253699 -0.39346608  0.7634227  0.90099341
## NDWI      -0.60241581  0.09112130  0.33720921 -0.9074901 -0.97346023
## COBERTURA  0.68818444 -0.13615460 -0.07917889  0.5967216  0.39922200
##                NDWI_1       NDWI   COBERTURA
## Y          0.50263155 -0.6024158  0.68818444
## COS10     -0.02253699  0.0911213 -0.13615460
## Altitud   -0.39346608  0.3372092 -0.07917889
## SAVI       0.76342267 -0.9074901  0.59672161
## NDVI       0.90099341 -0.9734602  0.39922200
## NDWI_1     1.00000000 -0.8758713  0.31083025
## NDWI      -0.87587134  1.0000000 -0.45866726
## COBERTURA  0.31083025 -0.4586673  1.00000000
corrplot.mixed(corMat, number.cex=0.8, tl.cex = 0.9, tl.col = "black", 
               outline=FALSE, mar=c(0,0,2,2), upper="square", bg=NA)

Esta matriz de correlacion muestra que las variables seleccionadas tienen una muy baja correlación con el objeto de estudio. Sinembargo se realizara la regresion con fines practicos de aprendizaje.

RegresiÓn

Validacion cruzada K-fold

K-fold, proporciona indices de tren-prueba para dividir los datos. Divide el conjunto de datos en K pliegues consecutivos y cada pliegue se usa como una validación, y los k-1 pliegues forman el conjunto de entrenamiento. El resultado de esta funcion es la lista con indices de prueba. x son las filas de indices.

kfoldSplit <- function(x, k=10, train=TRUE){
  x <- sample(x, size = length(x), replace = FALSE)
  out <- suppressWarnings(split(x, factor(1:k)))
  if(train) out <- lapply(out, FUN = function(x, len) (1:len)[-x], len=length(unlist(out)))
  return(out)
}

## Residuos de la regresion del objeto RF
resid.RF <- function(x) return(x$y - x$predicted)

Es necesario definir unos parametros adicionales para obtener los splits de entrenamiento con la funcion kfoldSplit e inicializar la matriz que almacenara los valores RMSE

set.seed(1982)

k <- 10

kfolds <- kfoldSplit(1:nrow(COS), k = 10, train = TRUE)

evalData <- matrix(NA, nrow=k, ncol=7, dimnames = list(1:k, c("OK","RF","GLM","GAM","RF_OK","GLM_OK","GAM_OK")))

Ok = Interpolación Ordinary Kriging que consiste en considerar la estimacion de los puntos que contienen la informacion, para este caso de COS, como una combinacion lineal de las observaciones, bajo el cririo de que dicha estimación es optima

Regresiones: RF = Random Forest GLM = Modelo lineal generalizdo GAM = Modelo aditivo generalizado

Regresion kriging: GLM + residuos de OK GAM + residuos de OK RF + residuos de OK

Regression residuals from RF object

resid.RF <- function(x) return(x$y - x$predicted)

MODELO

for(i in 1:k){
  
  cat("K-fold...",i,"of",k,"....\n")
  
  # Indices TRAIN como entero
  idx <- kfolds[[i]]
  
  # Indices TRAIN como vector booleano
  idxBool <- (1:nrow(COS)) %in% idx
  
  # Datos de prueba observados para la variable objetivo (Carbono organico en suelo).
  obs.test <- COS[!idxBool, "COS10"]
  
  
  ## ----------------------------------------------------------------------------- ##
  ## Kriging Ordinario ----
  ## ----------------------------------------------------------------------------- ##
  
  # Variograma
  formMod <- COS10 ~ 1
  mod <- vgm(model  = "Exp", psill  = 3, range  = 100, nugget = 0.5)
  variog <- variogram(formMod, valores[idxBool, ])
  
  # Ajuste de variograma por mínimo cuadrado
  variogFitOLS<-fit.variogram(variog, model = mod,  fit.method = 6)
  plot(variog, variogFitOLS, main="OLS Model")
  
  # Predicciones kriging 
  OK <- krige(formula = formMod ,
              locations = valores[idxBool, ], 
              model = variogFitOLS,
              newdata = valores[!idxBool, ],
              debug.level = 0)
  
  ok.pred.test <- OK@data$var1.pred
  evalData[i,"OK"] <- sqrt(mean((ok.pred.test - obs.test)^2))
  
  
  
  ## ----------------------------------------------------------------------------- ##
  ## Calibración RF----
  ## ----------------------------------------------------------------------------- ##
  
  RF_COS <- randomForest(y = COS[idx, "COS10"], 
                     x = COS[idx, c("Altitud","SAVI", "NDVI", "NDWI", "COBERTURA")],
                     ntree = 500,
                     mtry = 2)
  
  rf.pred.test <- predict(RF_COS, newdata = COS[-idx,], type="response")
  evalData[i,"RF"] <- sqrt(mean((rf.pred.test - obs.test)^2))
  
  # Kriging ordinario de los residuos de Random Forest
  
  statPointsTMP <- valores[idxBool, ]
  statPointsTMP@data <- cbind(statPointsTMP@data, residRF = resid.RF(RF_COS))
  
  formMod <- residRF ~ 1
  mod <- vgm(model  = "Exp", psill  = 0.6, range  = 10, nugget = 0.01)
  variog <- variogram(formMod, statPointsTMP)
  
  # Ajuste de variograma por mínimo cuadrado
  variogFitOLS<-fit.variogram(variog, model = mod,  fit.method = 6)
  #plot(variog, variogFitOLS, main="OLS Model")
  
  # Predicciones kriging
  RF.OK <- krige(formula = formMod ,
                 locations = statPointsTMP, 
                 model = variogFitOLS,
                 newdata = valores[!idxBool, ],
                 debug.level = 0)
  
  rf.ok.pred.test <- rf.pred.test + RF.OK@data$var1.pred
  evalData[i,"RF_OK"] <- sqrt(mean((rf.ok.pred.test - obs.test)^2))
  
  
  
  ## ----------------------------------------------------------------------------- ##
  ## Calibration GLM----
  ## ----------------------------------------------------------------------------- ##
  
  GLM <- glm(formula = COS10 ~ Altitud + SAVI + NDVI + NDWI + COBERTURA, data = COS[idx, ])
  
  glm.pred.test <- predict(GLM, newdata = COS[-idx,], type="response")
  evalData[i,"GLM"] <- sqrt(mean((glm.pred.test - obs.test)^2))
  
  # Kriging ordinario de los residuos de GLM
  #
  statPointsTMP <- valores[idxBool, ]
  statPointsTMP@data <- cbind(statPointsTMP@data, residGLM = resid(GLM))
  
  formMod <- residGLM ~ 1
  mod <- vgm(model  = "Exp", psill  = 0.4, range  = 10, nugget = 0.01)
  variog <- variogram(formMod, statPointsTMP)
  
  # Ajuste de variograma por minimo cuadrado
  variogFitOLS<-fit.variogram(variog, model = mod,  fit.method = 6)
  #plot(variog, variogFitOLS, main="OLS Model")
  
  #  Prediccion kriging
  GLM.OK <- krige(formula = formMod ,
                  locations = statPointsTMP, 
                  model = variogFitOLS,
                  newdata = valores[!idxBool, ],
                  debug.level = 0)
  
  glm.ok.pred.test <- glm.pred.test + GLM.OK@data$var1.pred
  evalData[i,"GLM_OK"] <- sqrt(mean((glm.ok.pred.test - obs.test)^2))
  
  
  
  ## ----------------------------------------------------------------------------- ##
  ## Calibracion GAM ----
  ## ----------------------------------------------------------------------------- ##
  
  GAM <- gam(formula = COS10 ~ Altitud + SAVI + NDVI + NDWI + COBERTURA, data = COS[idx, ])
  
  gam.pred.test <- predict(GAM, newdata = COS[-idx,], type="response")
  evalData[i,"GAM"] <- sqrt(mean((gam.pred.test - obs.test)^2))
  
  # Kriging ordinario de los residuos de GAM
  #
  statPointsTMP <- valores[idxBool, ]
  statPointsTMP@data <- cbind(statPointsTMP@data, residGAM = resid(GAM))
  
  formMod <- residGAM ~ 1
  mod <- vgm(model  = "Exp", psill  = 0.3, range  = 10, nugget = 0.01)
  variog <- variogram(formMod, statPointsTMP)
  
  # Ajuste de variograma por minimo cuadrado
  variogFitOLS<-fit.variogram(variog, model = mod,  fit.method = 6)
  #plot(variog, variogFitOLS, main="OLS Model")
  
  # Predicciones kriging
  GAM.OK <- krige(formula = formMod ,
                  locations = statPointsTMP, 
                  model = variogFitOLS,
                  newdata = valores[!idxBool, ],
                  debug.level = 0)
  
  gam.ok.pred.test <- gam.pred.test + GAM.OK@data$var1.pred
  evalData[i,"GAM_OK"] <- sqrt(mean((gam.ok.pred.test - obs.test)^2))
  
}
## K-fold... 1 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 2 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 3 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 4 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 5 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 6 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 7 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 8 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 9 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## K-fold... 10 of 10 ....
## Warning in fit.variogram(variog, model = mod, fit.method = 6): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit
## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix

## Warning in fit.variogram(variog, model = mod, fit.method = 6): linear model
## has singular covariance matrix
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

Union de bandas

Altitud <- raster("DEM_zona.tif")
cobertura <- raster("Cober.tif")
REPROYECTAR
crs(Altitud)
## CRS arguments:
##  +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666
## +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs
crs(savi)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
Altitud <- projectRaster(Altitud, crs = crs(savi))
crs(Altitud)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
Altitud <- resample(Altitud, savi)

Altitud <- crop(Altitud, extent(lalibertad))
Altitud <- mask(Altitud, lalibertad)
plot(Altitud)

cobertura <- crop(cobertura, extent(lalibertad))
cobertura <- mask(cobertura, lalibertad)
plot(cobertura)

cobertura <- resample(cobertura, savi)
library(mgcv)
library(gstat)
library(randomForest)


union <- stack(Altitud, savi, ndvi, ndwi)
union
## class       : RasterStack 
## dimensions  : 214, 141, 30174, 4  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 668895, 673125, 448485, 454905  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :     DEM_zona,      layer.1,      layer.2,      layer.3 
## min values  : 298.56344625,  -0.01979234,  -0.05802729,  -0.75564268 
## max values  : 340.06785651,   0.65928459,   0.86359717,   0.09409625
names(union) <- c("Altitud", "SAVI", "NDVI", "NDWI")

round(apply(evalData,2,FUN = function(x,...) c(mean(x,...),sd(x,...))),3)
##         OK    RF   GLM   GAM RF_OK GLM_OK GAM_OK
## [1,] 6.593 6.387 6.566 6.566 6.389  6.566  6.566
## [2,] 2.071 1.705 2.056 2.056 1.711  2.056  2.056
GAM <- gam(formula = COS10 ~ Altitud + SAVI + NDVI + NDWI, data = COS)

rstPredGAM <- predict(union, GAM, type="response")

rstPixDF <- as(union[[1]], "SpatialPixelsDataFrame")

Crear un objeto Spatial points para almacenar los residuos de GAM

statPointsTMP <- valores
crs(statPointsTMP) <- crs(rstPixDF)
statPointsTMP@data <- cbind(statPointsTMP@data, residGAM = resid(GAM))

Definicion de los parametros de kriging y ajuste del variograma utilizando OLS

formMod <- residGAM ~ 1
mod <- vgm(model  = "Exp", psill  = 0.15, range  = 10, nugget = 0.01)
variog <- variogram(formMod, statPointsTMP)
variogFitOLS <- fit.variogram(variog, model = mod,  fit.method = 6)
## Warning in fit.variogram(variog, model = mod, fit.method = 6): singular
## model in variogram fit

Ploteo de resultados. Las variables predictoras no se ajustan para el calculo de COS

plot(variog, variogFitOLS, main="Semi-variogram of GAM residuals")

residKrigMap <- krige(formula = formMod ,
                      locations = statPointsTMP, 
                      model = variogFitOLS,
                      newdata = rstPixDF)
## [using ordinary kriging]
residKrigRstLayer <- as(residKrigMap, "RasterLayer")

gamKrigMap <- rstPredGAM + residKrigRstLayer

plot(gamKrigMap, main="Carbono Orgánico en el suelo (GAM regression-kriging)",
     xlab="Longitude", ylab="Latitude", cex.main=0.8, cex.axis=0.7, cex=0.8)

RESULTADOS y CONCLUSIONES

Como se puede observar en el semivariograma los resultados obtenidos en esta practica de kriging no son los esperados, ya que las variables predictoras seleccionadas no se ajustan para el calculo de COS. Es importante tener datos q

DISCUSIÓN

Este ejercicio como práctica para entender los procesos y la información necesaria para llevar a cabo una regresión kriging es pertinente. Sin embargo los resultados obtenidos no son optimos como aplicación real ya que varios factores influyen en la calidada de los resultados dentro de estos se pueden mencionar: En primer lugar que la toma de los datos de muestra para la regresión no coinciden con la fecha de la imagen utilizada, ya que las imagenes disponibles para la zona en la epoca requerida no tienen caracteristicas adecuadas para el analisis (cobertura de nubes sobre la zona de estudio) Otro factor importante es que las variables predictoras utilizadas no son las adecuadas para calcular el carbono organico en el suelo

REFERENCIAS

  • FAO. 2018. Regresssion Kriging. Global Soil Partnership training on Digital Soil Organic Carbon Mapping
  • Gonzaga, C. 2014. Aplicación de Índices de Vegetación Derivados de Imágenes Satelitales Landsat 7 ETM+ y ASTER para la Caracterización de la Cobertura Vegetal en la Zona Centro de la Provincia De Loja, Ecuador.
  • Hengl, T., Heuvelink, G y Rossiter, D. 2007. About regresssion-kriging: From equations to case studies. Computers & Geosciences.