#This program calculates the deforestation rate for for 1985, 1995, 2005, 2015, generates a barplot and adds a linear regression to the plot.

#It uses a different method to calculate the number of forested and deforested points

#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 = -60.0
lon2 = -50.0
lat1 = -15.0
lat2 = -5.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)

#another method to count the  numbers of natural and altered cells
natural = freq(mydata,value=0)  #forested
altered = freq(mydata,value=1)  #deforested
water   = freq(mydata,value=2)  #water

#another way to calculate the total number of cells in the rasters
total_cells = natural+altered+water

#calculate the deforestation rate for each year
deforest_rate = 100*(altered/total_cells)
deforest_rate
## biomas_1985_rs biomas_1995_rs biomas_2005_rs biomas_2015_rs 
##       7.444583      12.000347      21.585764      24.033889
#barplot
years =   c("1985","1995","2005","2015")
barplot(deforest_rate ~ years, col="red", 
        main="Deforestation rate in study area", 
        xlab="Year",
        ylab="%",
        ylim=c(0,40))
box()

#use function lm to find the linear regression for the deforest_rate 
yy=seq(1,4)  #number of years (1,2,3,4)
reg=lm(deforest_rate~yy)  #linear regression

#add linear regression line to barplot
abline(reg,col="darkred",lwd=2)

#print the statistical summary of the linear regression on screen
reg_model = summary(reg)  
print(reg_model)
## 
## Call:
## lm(formula = deforest_rate ~ yy)
## 
## Residuals:
## biomas_1985_rs biomas_1995_rs biomas_2005_rs biomas_2015_rs 
##        0.08144       -1.29813        2.35195       -1.13526 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   1.4278     2.5267   0.565   0.6289  
## yy            5.9353     0.9226   6.433   0.0233 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.063 on 2 degrees of freedom
## Multiple R-squared:  0.9539, Adjusted R-squared:  0.9309 
## F-statistic: 41.39 on 1 and 2 DF,  p-value: 0.02332