#
#This program generates the 1985 deforestation map for a study area
#defined by longitude and latitude coordinates (lon1,lon2,lat1,lat2)
#
#Load libraries
library(raster)
library(rasterVis)
##read the file biomas_1985_rs.tif with Brazil land cover for 1985
#file contains 31 land cover type
bio1985=raster("biomas_1985_rs.tif")
#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(bio1985,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)
#use ratify to convert from continuous to categorical data
mydata = ratify(mydata)
rat = levels(mydata)[[1]]
rat$landcover = c("Natural","Altered","Water")
levels(mydata) = rat
#set colors for the map
#http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
cc=c("darkolivegreen4","bisque3","blue")
#generate land cover map
p1 = levelplot(mydata,margin=F, colorkey=T,
par.settings=custom.theme(region=cc),
scales = list(cex=1, tck=-0.5),
xlab="",ylab="",main="Deforestation map")
#plot map
print(p1)
