Analisis de cobertura fraccional
Leonardo Hurtado
26 de noviembre de 2017
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: splibrary(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewerlibrary(tmap)
setwd("H:\\percepcion avanzada\\final")
dir()## [1] "codigo original.R"
## [2] "codigo_final.Rmd"
## [3] "codigo_final_files"
## [4] "codigofinal.Rmd"
## [5] "delimitacion.dbf"
## [6] "delimitacion.prj"
## [7] "delimitacion.sbn"
## [8] "delimitacion.sbx"
## [9] "delimitacion.shp"
## [10] "delimitacion.shx"
## [11] "dem"
## [12] "indices"
## [13] "LC80060582014035LGN00_refl.tfw"
## [14] "LC80060582014035LGN00_refl.tif"
## [15] "LC80060582014035LGN00_refl.tif.aux.xml"
## [16] "LC80060582014035LGN00_refl.tif.ovr"
## [17] "msavi2.tfw"
## [18] "msavi2.tif"
## [19] "msavi2.tif.aux.xml"
## [20] "msavi2.tif.ovr"
## [21] "muestreo.dbf"
## [22] "muestreo.prj"
## [23] "muestreo.sbn"
## [24] "muestreo.sbx"
## [25] "muestreo.shp"
## [26] "muestreo.shp.LEONARDO.6700.5924.sr.lock"
## [27] "muestreo.shx"
## [28] "ndvi.tfw"
## [29] "ndvi.tif"
## [30] "ndvi.tif.aux.xml"
## [31] "ndvi.tif.ovr"
## [32] "rdvi.tfw"
## [33] "rdvi.tif"
## [34] "rdvi.tif.aux.xml"
## [35] "rdvi.tif.ovr"
## [36] "sarvi.tfw"
## [37] "sarvi.tif"
## [38] "sarvi.tif.aux.xml"
## [39] "sarvi.tif.ovr"
## [40] "tndvi.tfw"
## [41] "tndvi.tif"
## [42] "tndvi.tif.aux.xml"
## [43] "tndvi.tif.ovr"
## [44] "tps20171119t185412_srtm30m.tfw"
## [45] "tps20171119t185412_srtm30m.tif"
## [46] "tps20171119t185412_srtm30m.tif.aux.xml"
## [47] "tps20171119t185412_srtm30m.tif.ovr"
## [48] "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 : H:\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 1mapTheme <- 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 : H:\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 : H:\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 : H:\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 : H:\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 : H:\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 <- aspectdeclaracion 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
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 colorsplt <- 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))