#
#This program generates deforestation map for 1985, 1995, 2005, 2015
#for a study area and calculates the deforestation rate for each year.
#Load raster library
library(raster)
##read land cover datafiles for 1985, 1995, 2005, and 2015
m1=raster("biomas_1985_rs.tif")
m2=raster("biomas_1995_rs.tif")
m3=raster("biomas_2005_rs.tif")
m4=raster("biomas_2015_rs.tif")
#stack 4 years into one variable biome
biome = stack(m1,m2,m3,m4)
#Define a study area (domain) here
lon1 = -50.0
lon2 = -40.0
lat1 = -20.0
lat2 = -10.0
domain = extent(lon1,lon2,lat1,lat2)
#crop data to the study area defined above
mydata = crop(biome,domain)
#use the function reclassify to simplify the data
#from 31 to 3 cover types: natural, deforested, and water
m = c(-Inf,0,NA, 0.5,13,0, 14,21,1, 22,Inf,2)
rcm = matrix(m, ncol=3, byrow=T)
mydata=reclassify(mydata,rcm)
#find the total number of cells in the study area
tmp = mydata
tmp[tmp >= 0] = 1
#use function cellStats to calculate the number of cells
total_cells = cellStats(tmp,sum)
#find the number of altered (deforested) cells in the study area
tmp = mydata
tmp[tmp == 0] = 0 #set original land cover cells to zero
tmp[tmp == 1] = 1 #deforested cells set to 1
tmp[tmp == 2] = 0 #set water cells to zero
#use function cellStats to calculate the number of cells
deforest_cells= cellStats(tmp,sum)
#calculate deforestation rate (deforested cells/total cells)
deforest_rate = 100*(deforest_cells/total_cells)
#barplot
dates = c("1985","1995","2005","2015")
barplot(deforest_rate ~ dates, col="red",
main="Deforestation rate in study area",
xlab="Year",
ylab="%",
ylim=c(0,60))
box()
