#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