NDVI-entre-2013-a-2018.R

## vamos realizar a comparação entre dois anos, de 2013 a 2018, 
## para saber se houve ganho ou déficit de vegetação

library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/cided/Documents/R/win-library/3.6/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/cided/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1

library(rgeos)

## rgeos version: 0.4-3, (SVN revision 595)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

library(RColorBrewer)

# Bandas 2013
red2013=raster("Recorte Landsat/2013/banda_red.tif")
nir2013=raster("Recorte Landsat/2013/banda_nir.tif")

# bandas 2018
red2018=raster("Recorte Landsat/2018/banda_red.tif")
nir2018=raster("Recorte Landsat/2018/banda_nir.tif")

# cálculo de NDVI

ndvi2013 <- overlay(nir2013, red2013, fun=function(x,y){(x-y)/(x+y)})
ndvi2018 <- overlay(nir2018, red2018, fun=function(x,y){(x-y)/(x+y)})
par(mfrow=c(1,2))
plot(ndvi2013, col = rev(terrain.colors(15)), main = "NDVI 2013")
plot(ndvi2018, col = rev(terrain.colors(15)), main = "NDVI 2018")

par(mfrow=c(1,1))
plot((ndvi2013-ndvi2018), col = rev(terrain.colors(15)), main = "Vegetação Perdida")

# histogramas para comparar frequência por pixel
par(mfrow=c(1,2))
hist(ndvi2013,
     main="NDVI 2013",
     xlab="ndvi",col="red", ylab="Frequência por Pixels")
hist(ndvi2018,
     main="NDVI 2018",
     xlab="ndvi",col="red", ylab="Frequência por Pixels")

# extrair somente a cobertura vegetal
veg2013 <- calc(ndvi2013, function(x) {
  x[x < 0.3] <- NA
  return(x)
})
veg2018 <- calc(ndvi2018, function(x) {
  x[x < 0.3] <- NA
  return(x)
})
dfveg <- calc((ndvi2013-ndvi2018), function(x) {
  x[x < 0.3] <- NA
  return(x)
})

coresveg=brewer.pal(9,"Greens")
plot(veg2013, main = "Cobertura Vegetal 2013",col = coresveg)
plot(veg2018, main = "Cobertura Vegetal 2018",col = coresveg)

par(mfrow=c(1,1))
plot(dfveg, main = "Cobertura Vegetal",col = coresveg)

# aumentar para 0.6 de NDVI e mudar o estílo das cores
cobertura2013 <- calc(ndvi2013, function(x) {
  x[x < 0.4] <- NA
  return(x)
})
cobertura2018 <- calc(ndvi2018, function(x) {
  x[x < 0.4] <- NA
  return(x)
})
dfcorb <- calc((ndvi2013-ndvi2018), function(x) {
  x[x < 0.4] <- NA
  return(x)
})

par(mfrow=c(1,2))
plot(cobertura2013, col="darkgreen", main="NDVI_0.60 2013")
plot(cobertura2018, col="darkgreen", main="NDVI_0.60 2018")

par(mfrow=c(1,1))
plot(dfcorb, col="darkgreen", main="NDVI_0.60 Diferença")

# Para gravar em Geotiff:
writeRaster(cobertura2013,"NDVI_0.60 2013.tif",overwrite=T, format="GTiff")
writeRaster(cobertura2018,"NDVI_0.60 2018.tif",overwrite=T, format="GTiff")