cobertura fraccional

El presente script muestra la manera en que los indices espectrales de vegetacion y suelos muestran el comportamiento de la vegetacion y se compara con los resultados propios de la cobertura fraccional.

Se tiene en cuenta variables del terreno como la pendiente y el aspecto partiendo de la hipotesis que tiene una reaccion directa con la formacion del dosel arboreo. Se seleccionan cinco indices espectrales que representan la formacion vegetal.

La zona de estudio es la zona sur del departamento del Meta donde el cambio engre vegtacion y suelos expuestos marca una tendencia y un reto en la funcionalidad de los indices espectrales.

Los datos pueden descargarse de https://goo.gl/7JxtDu para poder replicar el ejercicio.

Area de estudio

El area de estudio comprende la escena 6-58 de landsat que hace referencia a los departamentos del Meta, Caqueta y Guaviare. Esta zona muestra la transicion entre la region de la Orinoquia y la Amazonia.

library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
library(tmap)
setwd("E:\\percepcion avanzada\\final")
dir()
##  [1] "cobertura_fraccional_2014.tif"         
##  [2] "codigo original.R"                     
##  [3] "codigo_final.html"                     
##  [4] "codigo_final.Rmd"                      
##  [5] "codigo_final_files"                    
##  [6] "codigofinal.Rmd"                       
##  [7] "delimitacion.dbf"                      
##  [8] "delimitacion.prj"                      
##  [9] "delimitacion.sbn"                      
## [10] "delimitacion.sbx"                      
## [11] "delimitacion.shp"                      
## [12] "delimitacion.shx"                      
## [13] "dem"                                   
## [14] "indices"                               
## [15] "LC80060582014035LGN00_refl.tfw"        
## [16] "LC80060582014035LGN00_refl.tif"        
## [17] "LC80060582014035LGN00_refl.tif.aux.xml"
## [18] "LC80060582014035LGN00_refl.tif.ovr"    
## [19] "msavi2.tfw"                            
## [20] "msavi2.tif"                            
## [21] "msavi2.tif.aux.xml"                    
## [22] "msavi2.tif.ovr"                        
## [23] "muestreo.dbf"                          
## [24] "muestreo.prj"                          
## [25] "muestreo.sbn"                          
## [26] "muestreo.sbx"                          
## [27] "muestreo.shp"                          
## [28] "muestreo.shx"                          
## [29] "muestreo2.dbf"                         
## [30] "muestreo2.lock"                        
## [31] "muestreo2.prj"                         
## [32] "muestreo2.sbn"                         
## [33] "muestreo2.sbx"                         
## [34] "muestreo2.shp"                         
## [35] "muestreo2.shx"                         
## [36] "ndvi.tfw"                              
## [37] "ndvi.tif"                              
## [38] "ndvi.tif.aux.xml"                      
## [39] "ndvi.tif.ovr"                          
## [40] "rdvi.tfw"                              
## [41] "rdvi.tif"                              
## [42] "rdvi.tif.aux.xml"                      
## [43] "rdvi.tif.ovr"                          
## [44] "rf_fraccional_index.tif"               
## [45] "rf_nofraccional_2014.tif"              
## [46] "rsconnect"                             
## [47] "sarvi.tfw"                             
## [48] "sarvi.tif"                             
## [49] "sarvi.tif.aux.xml"                     
## [50] "sarvi.tif.ovr"                         
## [51] "tndvi.tfw"                             
## [52] "tndvi.tif"                             
## [53] "tndvi.tif.aux.xml"                     
## [54] "tndvi.tif.ovr"                         
## [55] "tps20171119t185412_srtm30m.tfw"        
## [56] "tps20171119t185412_srtm30m.tif"        
## [57] "tps20171119t185412_srtm30m.tif.aux.xml"
## [58] "tps20171119t185412_srtm30m.tif.ovr"    
## [59] "tps20171119t185412_srtm30m.tif.vat.dbf"
cob_frac <- shapefile("muestreo.shp")
delimitacion <- shapefile("delimitacion.shp")
plot(delimitacion)

plot(cob_frac)

Metodologia.

El estudio se propone en la realizacion de arboles aleatorios a partir del muestreo de puntos dentro del area de estudio tomando como valores los indices espectrales de vegetacion. El estudio se realiza para el año 2014 con datos tematicos y fisicos.

Datos de cobertura fraccional

Los datos son con una resolucion de 30 metros y resampleados para mantener la misma unidad de pixel.

Ahora veremos un poco de la estadistica basica de los datos insumos:

hist(cob_frac$b4_b)

Digital Elevation

A continuacion la visualizacion del DEM para la zona de estudio:

dem2 <- raster("tps20171119t185412_srtm30m.tif")
dem2
## class       : RasterLayer 
## dimensions  : 7203, 6526, 47006778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 782013, 977793, 205202, 421292  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\tps20171119t185412_srtm30m.tif 
## names       : tps20171119t185412_srtm30m 
## values      : 0, 613  (min, max)
## attributes  :
##         ID OBJECTID Value Count
##  from:   0        1     0     1
##  to  : 496      497   496     1
mapTheme <- rasterTheme(region=rev(brewer.pal(11,"RdYlGn")))
plt <- levelplot(dem2, margin=F, par.settings=mapTheme,
                 main="SRTM Digital Elevation Model (30m)")
plt + layer(sp.lines(delimitacion, col="gray", lwd=0.5))

Como se requiere una reproyeccion de los datos se ejecuta las siguientes instrucciones:

library(sp)
proy_dem <- proj4string(dem2)
crs_dem <- CRS(proy_dem)
frac_rep <- spTransform(cob_frac, crs_dem)
plot(dem2, col=terrain.colors(10), axes = TRUE)
plot(frac_rep, pch=21, main="Localizacion de puntos en el area de estudio",cex=0.5, add=TRUE)

insumos principales. Indices espectrales

Los indices espectrales permiten identificar lasvariables espectrales que permiten diferenciar la cobertura fraccional, separar el dosel arboreo de vegetacion no fotosintetica y de los suelos.

ndvi = raster("ndvi.tif")
ndvi
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\ndvi.tif 
## names       : ndvi 
## values      : -0.1496321, 0.9561768  (min, max)
rdvi = raster("rdvi.tif")
rdvi
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\rdvi.tif 
## names       : rdvi 
## values      : -7.40036, 68.5067  (min, max)
sarvi = raster("sarvi.tif")
sarvi
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\sarvi.tif 
## names       : sarvi 
## values      : -0.5252225, 1.598389  (min, max)
tndvi = raster("tndvi.tif")
tndvi
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\tndvi.tif 
## names       : tndvi 
## values      : 0, 1.206722  (min, max)
msavi2  = raster("msavi2.tif")
msavi2
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : E:\percepcion avanzada\final\msavi2.tif 
## names       : msavi2 
## values      : -1, 0.9775933  (min, max)

Variables topograficas

La cobertura boscosa tiene diferentes grados de pendientes y el terreno permite que un ecosistema boscoso pueda coexistir o no. La libreria dynatopmodel permite modelar estas variables de terreno.

###install.packages("dynatopmodel")
library(dynatopmodel)

Como el DEM no cuenta con las mimas dimensiones de los indices espectrales se realiza una reproyeccion.

dem <- resample(dem2, msavi2, method="ngb")
dem
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : tps20171119t185412_srtm30m 
## values      : 0, 613  (min, max)
####pendiente
##slope = terrain(dem, opt=c('slope'), unit='degrees', neighbors=8)
##asp = terrain(dem, opt=c('aspect'), unit='degrees', neighbors=8)

### the upslope function has been modified
##  based on Robert Hijmans advise

upslope <- function (dem, log = TRUE, atb = FALSE, deg = 0.12, fill.sinks = TRUE) 
{
  if (!all.equal(xres(dem), yres(dem))) {
    stop("Raster has differing x and y cell resolutions. Check that it is in a projected coordinate system (e.g. UTM) and use raster::projectRaster to reproject to one if not. Otherwise consider using raster::resample")
  }
  if (fill.sinks) {
    capture.output(dem <- invisible(raster::setValues(dem, 
                                                      topmodel::sinkfill(raster::as.matrix(dem), res = xres(dem), 
                                                                         degree = deg))))
  }
  topidx <- topmodel::topidx(raster::as.matrix(dem), res = xres(dem))
  a <- raster::setValues(dem, topidx$area)
  if (log) {
    a <- log(a)
  }
  if (atb) {
    atb <- raster::setValues(dem, topidx$atb)
    a <- addLayer(a, atb)
    names(a) <- c("a", "atb")
  }
  return(a)
}

create_layers <- function (dem, fill.sinks = TRUE, deg = 0.1) 
{
  layers <- stack(dem)
  message("Building upslope areas...")
  a.atb <- upslope(dem, atb = TRUE, fill.sinks = fill.sinks, 
                   deg = deg)
  layers <- addLayer(layers, a.atb)
  names(layers) <- c("dem", "a", "atb")
  return(layers)
}

topo <- create_layers(dem, fill.sinks = TRUE, deg = 0.10)
## Building upslope areas...
# topo is a multi-band raster (stack) comprising, in order, 
#the filled elevations, upslope area and topographic wetness
#index values.

### calculo de variables topograficas
slope_n <- terrain(dem, "slope", "tangent", 8)

slope <- slope_n * 100

mapTheme <- rasterTheme(region=rev(brewer.pal(11,"RdYlBu")))
plt <- levelplot(slope, margin=F, par.settings=mapTheme,
                 main="Slope")
plt + layer(sp.lines(delimitacion, col="gray", lwd=0.5))

aspect <- terrain(dem, "aspect", "tangent", 8)
aspect
## class       : RasterLayer 
## dimensions  : 7251, 6582, 47726082  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 781325, 978785, 204675.2, 422205.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : aspect 
## values      : 0, 6.283185  (min, max)
mapTheme <- rasterTheme(region=rev(brewer.pal(11,"Spectral")))
plt <- levelplot(aspect, margin=F, par.settings=mapTheme,
                 main="Aspect")
plt + layer(sp.lines(delimitacion, col="gray", lwd=0.5))

topo$slope <- slope

topo$aspect <- aspect

declaracion de variables tematicas

NVDI (Índice de vegetación de diferencia normalizada) Es un índice usado para estimar la cantidad, calidad y desarrollo de la vegetación con base a la medición, por medio de sensores remotos instalados comúnmente desde una plataforma espacial, de la intensidad de la radiación de ciertas bandas del espectro electromagnético que la vegetación emite o refleja.

Inclusion de las variables topograficas

topo$ndvi <- ndvi

topo$rdvi <- rdvi

topo$sarvi <- sarvi

topo$tndvi <- tndvi

topo$msavi2 <- msavi2

library(rasterVis)

proj <- CRS('+proj=longlat +ellps=WGS84')
levelplot(topo$dem, par.settings = BuRdTheme, xlab.top='DEM')

levelplot(topo$a, par.settings = BTCTheme, xlab.top='Upslope Area')

levelplot(topo$atb, par.settings = viridisTheme, xlab.top='TWI')

hist(ndvi)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(11,"RdYlGn"))
plt <- levelplot(ndvi, margin=F, par.settings=mapTheme, at = do.breaks(c(0.1,1),10),
                 main="NDVI para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=1))

hist(rdvi)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
plt <- levelplot(rdvi, margin=F, par.settings=mapTheme, at = do.breaks(c(-7,68),10),
                 main="RDVI para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

hist(sarvi)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
plt <- levelplot(sarvi, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,2),10),
                 main="SARVI para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

hist(tndvi)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
plt <- levelplot(tndvi, margin=F, par.settings=mapTheme, at = do.breaks(c(0.1,1.5),10),
                 main="TNDVI para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

hist(msavi2)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
plt <- levelplot(msavi2, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),10),
                 main="MSAVI2 para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

#  raw NBRT seems to be useless with no transformation
power2 <- function(x) {x**2}
ndvi2 <- calc(ndvi, power2)
hist(ndvi2, main="Square NDVI para cobertura fraccional")

#
mapTheme <- rasterTheme(region=brewer.pal(11,"YlOrRd"))
## Warning in brewer.pal(11, "YlOrRd"): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors
plt <- levelplot(ndvi2, margin=F, par.settings=mapTheme,
                 at = do.breaks(c(0.1,2),10),
                 main="NDVI  para cobertura fraccional")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

preparacion de los datos

Se realiza la preparacion de los datos originales, proyectados y resampleados para iniciar el conjunto de datos de muestra para las predicciones

frac_rep$x <- frac_rep$longitud
frac_rep$y <- frac_rep$latitud
frac_pts <- frac_rep
df.frac_pts <- data.frame(frac_pts)
coordinates(df.frac_pts) <- ~longitud+latitud

#df.cob_pts@coords

### extraction of explanatory variables data
### at points where soil moisture was measured in the field

proj4string(df.frac_pts) <- proj4string(cob_frac)

df.frac_rep <- spTransform(df.frac_pts, crs_dem)

### extraction of explanatory variables data
### at points where soil moisture was measured in the field

data <- data.frame(coordinates(df.frac_rep),
                   frac_rep, 
                   extract(topo, df.frac_rep))

coordinates(data) <- ~longitud+latitud

proj4string(data) <- proj4string(dem)

generacion del conjunto de muestras

#Sample Indexes
indexes = sample(1:nrow(data), size=0.4*nrow(data))

# Split data
test = data[indexes,]
dim(test)  # 80 35
## [1] 120  19
train = data[-indexes,]
dim(train) # 120 35
## [1] 180  19
###se colocan las columnas del shp de muestreo
drops <- c("longitud","latitud")
drops
## [1] "longitud" "latitud"
ntrain <- train[,!(names(train) %in% drops)] #remove columns "AREA" and "PERIMETER"

ntrain[1,0:16]
## class       : SpatialPointsDataFrame 
## features    : 1 
## extent      : 812505.6, 812505.6, 252968.9, 252968.9  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       : OBJECTID, longitud.1, latitud.1, b4_b,        x,       y,        coords.x1,        coords.x2, optional, dem,                a,              atb,            slope,           aspect,              ndvi, ... 
## value       :        1,   -72.1905,   2.28591,   89, -72.1905, 2.28591, 812501.253008783, 252968.696437463,        1, 223, 7.15305163493748, 7.15305163493748, 3.77307714089062, 4.82304620155858, 0.897284507751465, ...

Prediccion de arboles aleatorios para determinar cobertura fraccional

se realiza la inclusion de las variables procesadas que podrian tener la informacion util para generar prediciones bajo diferentes escenarios de cobertura fraccional boscosa para el area de estudio.

set.seed(42)
#install.packages("randomForest")
library(randomForest)
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
xx <- data.frame(ntrain$b4_b, ntrain$dem, ntrain$a, ntrain$atb, ntrain$slope, ntrain$aspect, ntrain$ndvi, ntrain$msavi2, ntrain$rdvi, ntrain$sarvi, ntrain$tndvi) 

xx <- xx[complete.cases(xx),]

y <- xx[1]

x <- xx[-1]

fit1 <- randomForest(b4_b ~ dem + a + atb + slope + aspect + ndvi + rdvi + sarvi + msavi2, data=train, 
                     ntrees=5000, mtry=5, na.action=na.omit, importance=TRUE)

print(fit1) # view results 
## 
## Call:
##  randomForest(formula = b4_b ~ dem + a + atb + slope + aspect +      ndvi + rdvi + sarvi + msavi2, data = train, ntrees = 5000,      mtry = 5, importance = TRUE, na.action = na.omit) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##           Mean of squared residuals: 863.3144
##                     % Var explained: 65.41
importance(fit1) # importance of each predictor
##           %IncMSE IncNodePurity
## dem     2.5654441      3569.351
## a       0.9782655      3031.033
## atb    -2.9615025      8151.905
## slope  -0.2225078      3371.366
## aspect  3.0118966      5728.003
## ndvi   15.6757849     42887.141
## rdvi    3.3452772     10178.896
## sarvi  14.4600662     26076.685
## msavi2 14.8633204     38459.847
par(mar=c(1,1,1,1))

varImpPlot(fit1,type=2)

rf1 <- predict(fit1, test)

summary(rf1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    67.5   110.3   122.3   126.3   131.6   237.3      75
compare1 <- cbind (actual=test$b4_b, rf1) 

#
mean (apply(compare1, 1, min, na.rm=T)/apply(compare1, 1, max, na.rm=T)) # calculate accuracy
## [1] 0.8996086

construccion del arbol aleatorio caracteristico de las variables.

#install.packages("party")
library(party)
## Loading required package: grid
## Loading required package: mvtnorm
## Loading required package: modeltools
## Loading required package: stats4
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
cf <- cforest(b4_b ~ dem + a + atb + slope + aspect + ndvi + rdvi + sarvi + msavi2 + tndvi
              , data = train) 
pt <- party:::prettytree(cf@ensemble[[1]], names(cf@data@get("input"))) 
pt 
## 1) ndvi <= 0.8636503; criterion = 1, statistic = 13.193
##   2)*  weights = 0 
## 1) ndvi > 0.8636503
##   3) sarvi <= 1.256241; criterion = 0.999, statistic = 12.003
##     4)*  weights = 0 
##   3) sarvi > 1.256241
##     5) msavi2 <= 0.9451706; criterion = 0.983, statistic = 5.722
##       6) sarvi <= 1.301302; criterion = 1, statistic = 13.936
##         7) a <= 7.166483; criterion = 0.876, statistic = 2.372
##           8) a <= 6.802395; criterion = 0.474, statistic = 0.403
##             9)*  weights = 0 
##           8) a > 6.802395
##             10)*  weights = 0 
##         7) a > 7.166483
##           11)*  weights = 0 
##       6) sarvi > 1.301302
##         12)*  weights = 0 
##     5) msavi2 > 0.9451706
##       13) sarvi <= 1.325933; criterion = 0.983, statistic = 5.649
##         14) a <= 6.802395; criterion = 0.892, statistic = 2.588
##           15)*  weights = 0 
##         14) a > 6.802395
##           16)*  weights = 0 
##       13) sarvi > 1.325933
##         17) msavi2 <= 0.9557341; criterion = 0.961, statistic = 4.268
##           18)*  weights = 0 
##         17) msavi2 > 0.9557341
##           19)*  weights = 0
nt <- new("BinaryTree") 
nt@tree <- pt 
nt@data <- cf@data 
nt@responses <- cf@responses 
nt 
## 
##   Conditional inference tree with 0 terminal nodes
## 
## Response:  b4_b 
## Inputs:  dem, a, atb, slope, aspect, ndvi, rdvi, sarvi, msavi2, tndvi 
## Number of observations:  180 
## 
## 1) ndvi <= 0.8636503; criterion = 1, statistic = 13.193
##   2)*  weights = 0 
## 1) ndvi > 0.8636503
##   3) sarvi <= 1.256241; criterion = 0.999, statistic = 12.003
##     4)*  weights = 0 
##   3) sarvi > 1.256241
##     5) msavi2 <= 0.9451706; criterion = 0.983, statistic = 5.722
##       6) sarvi <= 1.301302; criterion = 1, statistic = 13.936
##         7) a <= 7.166483; criterion = 0.876, statistic = 2.372
##           8) a <= 6.802395; criterion = 0.474, statistic = 0.403
##             9)*  weights = 0 
##           8) a > 6.802395
##             10)*  weights = 0 
##         7) a > 7.166483
##           11)*  weights = 0 
##       6) sarvi > 1.301302
##         12)*  weights = 0 
##     5) msavi2 > 0.9451706
##       13) sarvi <= 1.325933; criterion = 0.983, statistic = 5.649
##         14) a <= 6.802395; criterion = 0.892, statistic = 2.588
##           15)*  weights = 0 
##         14) a > 6.802395
##           16)*  weights = 0 
##       13) sarvi > 1.325933
##         17) msavi2 <= 0.9557341; criterion = 0.961, statistic = 4.268
##           18)*  weights = 0 
##         17) msavi2 > 0.9557341
##           19)*  weights = 0
plot(nt, type="simple", main = "A Tree from the Fraccional Cover Random Forest")

construccion del primer escenario. Incluyendo variable de cobertura fraccional obtenida con ClassLite

#install.packages("devtools")
library(devtools)
#install.packages("reprtree")
#library(reprtree)

cobertura_frac <- randomForest(b4_b ~ dem + a + atb + slope + b4_b + aspect + ndvi + sarvi + tndvi + msavi2 + rdvi,
                           data=train, na.action=na.omit, 
                           importance=TRUE, ntree=5000, mtry = 2, do.trace=100)
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  100 |    688.4    27.58 |
##  200 |    690.2    27.65 |
##  300 |    677.2    27.13 |
##  400 |    683.8    27.40 |
##  500 |    683.6    27.39 |
##  600 |    671.1    26.89 |
##  700 |    680.4    27.26 |
##  800 |    689.6    27.63 |
##  900 |    688.2    27.57 |
## 1000 |    689.5    27.63 |
## 1100 |    686.9    27.52 |
## 1200 |    683.4    27.38 |
## 1300 |    684.8    27.44 |
## 1400 |    685.6    27.47 |
## 1500 |    686.9    27.52 |
## 1600 |    694.4    27.82 |
## 1700 |    698.5    27.99 |
## 1800 |    710.9    28.48 |
## 1900 |    707.9    28.36 |
## 2000 |    709.4    28.42 |
## 2100 |    704.8    28.24 |
## 2200 |    707.6    28.35 |
## 2300 |    704.1    28.21 |
## 2400 |    702.8    28.16 |
## 2500 |      704    28.21 |
## 2600 |    701.7    28.11 |
## 2700 |      702    28.13 |
## 2800 |      703    28.17 |
## 2900 |      704    28.20 |
## 3000 |    705.5    28.27 |
## 3100 |    705.6    28.27 |
## 3200 |    704.8    28.24 |
## 3300 |    699.8    28.04 |
## 3400 |    700.9    28.08 |
## 3500 |    702.6    28.15 |
## 3600 |      702    28.13 |
## 3700 |    705.5    28.27 |
## 3800 |    707.8    28.36 |
## 3900 |    707.1    28.33 |
## 4000 |    708.5    28.39 |
## 4100 |    708.7    28.40 |
## 4200 |    706.3    28.30 |
## 4300 |    705.8    28.28 |
## 4400 |    705.3    28.26 |
## 4500 |    706.3    28.30 |
## 4600 |    707.1    28.33 |
## 4700 |    708.1    28.37 |
## 4800 |    708.2    28.37 |
## 4900 |    709.6    28.43 |
## 5000 |    706.9    28.32 |
#### random forest prediction in the whole study area

r_frac <- predict(topo, fit1, filename="cobertura_fraccional_2014.tif", overwrite=T)

hist(r_frac)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(11,"YlOrRd"))
## Warning in brewer.pal(11, "YlOrRd"): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors
plt <- levelplot(r_frac, margin=F, par.settings=mapTheme, at = do.breaks(c(0,80),10), main="Fraccional Cover 2014")
plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

Construccion escenario sin tener datos de cobertura fraccional.

fit2 <- randomForest(b4_b ~ dem + a + atb + slope + aspect + ndvi + sarvi + tndvi + msavi2 + rdvi 
                     ,   data=train, 
                     ntrees=5000, mtry=5, na.action=na.omit, importance=TRUE)

print(fit2) # view results 
## 
## Call:
##  randomForest(formula = b4_b ~ dem + a + atb + slope + aspect +      ndvi + sarvi + tndvi + msavi2 + rdvi, data = train, ntrees = 5000,      mtry = 5, importance = TRUE, na.action = na.omit) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##           Mean of squared residuals: 865.7611
##                     % Var explained: 65.31
importance(fit2) # importance of each predictor
##           %IncMSE IncNodePurity
## dem    -0.5118891      3597.546
## a       1.7791450      3318.087
## atb    -1.5200742      6190.363
## slope   0.1268016      3359.958
## aspect  2.6304752      5868.394
## ndvi   11.0416428     24483.338
## sarvi  11.6816565     19461.573
## tndvi  13.0306587     31004.003
## msavi2 12.2147747     31785.557
## rdvi    3.1585561      7647.870
par(mar=c(1,1,1,1))

varImpPlot(fit2,type=2)

rf2 <- predict(fit2, test)

summary(rf2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   71.66  109.23  123.04  125.96  131.77  233.05      75
hist(rf2)

compare2a <- cbind (actual=test$b4_b, rf2) 

mean (apply(compare2a, 1, min, na.rm=T)/apply(compare2a, 1, max, na.rm=T)) # calculate accuracy
## [1] 0.9003458
#### random forest prediction in the whole study area

r_frac2 <- predict(topo, fit2, filename="rf_nofraccional_2014.tif", overwrite=T)

hist(r_frac2)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(7,"YlOrRd"))

plt <- levelplot(r_frac2, margin=F, par.settings=mapTheme, at = do.breaks(c(0,80),10), main="Result without Fraccional cover 2014")

plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))

Construccion escenario incluyendo unicamente indices espectrales de vegetacion

fit3 <- randomForest(b4_b ~ ndvi + sarvi + tndvi + msavi2 + rdvi,   data=train, 
                     ntrees=5000, mtry=5, na.action=na.omit, importance=TRUE)
print(fit3) # view results 
## 
## Call:
##  randomForest(formula = b4_b ~ ndvi + sarvi + tndvi + msavi2 +      rdvi, data = train, ntrees = 5000, mtry = 5, importance = TRUE,      na.action = na.omit) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##           Mean of squared residuals: 1477.863
##                     % Var explained: 30
importance(fit3) # importance of each predictor
##          %IncMSE IncNodePurity
## ndvi   16.729251      75038.94
## sarvi  34.231933      98952.13
## tndvi  14.452782      67321.59
## msavi2 14.320714      67931.70
## rdvi   -2.028578      41098.54
varImpPlot(fit3,type=2)

rf3 <- predict(fit3, test)

summary(rf3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   42.65  105.61  120.44  122.52  138.20  230.99       1
hist(rf3)

compare3a <- cbind (actual=test$b4_b, rf3) 

mean (apply(compare3a, 1, min, na.rm=T)/apply(compare3a, 1, max, na.rm=T)) # calculate accuracy
## [1] 0.8096318
#### random forest prediction in the whole study area

r_frac3 <- predict(topo, fit3, filename="rf_fraccional_index.tif", overwrite=T)

hist(r_frac3)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

mapTheme <- rasterTheme(region=brewer.pal(9,"YlGnBu"))

plt <- levelplot(r_frac3, margin=F, par.settings=mapTheme, at = do.breaks(c(0,80),10), main="Result Fraccional cover with vegetation index only")

plt + layer(sp.lines(delimitacion, col="red", lwd=0.5))