######################################################################
# Historical land use and land-use change in Great Britain 1930s-2007
# Code to bring in data and produce figures. Data are described in separate readme file.
# Alistair Auffret, August 2023
######################################################################

## Setting up.

# Assign colours for plotting. Colours follow those used in the historical Land Utilisation Survey maps
ara.col<-"tan4"
gra.col<-"darkseagreen"
urb.col<-"firebrick2"
wdl.col<-"forestgreen"
wat.col<-"deepskyblue1"
imp.gra.col<-"gold"

# Write function for colour transparency for use in the figures
transp<-function(col,alpha){
  trans.col<-adjustcolor(col,alpha.f = alpha)
  return(trans.col)
}


## 1. Summary of land use change
LU.GB<-read.csv("GB_LandUseChange_Data.csv") # bring in data

## Historical extents
sum(LU.GB$LUS.ara*LU.GB$hab.pix) / sum(LU.GB$hab.pix) # arable covered 22%
sum(LU.GB$LUS.gra*LU.GB$hab.pix) / sum(LU.GB$hab.pix) # grassland covered 65%
sum(LU.GB$LUS.urb*LU.GB$hab.pix) / sum(LU.GB$hab.pix) # urban covered 4%
sum(LU.GB$LUS.wdl*LU.GB$hab.pix) / sum(LU.GB$hab.pix) # woodland covered 6%

## Change at hectad level, seen in Fig 1.
boxplot(LU.GB[,c("chg.urb", "chg.wdl", "chg.ara", "chg.gra","chg.imp.gra")], col="white",names=c("Urban", "Woodland", "Arable", "Open", "Improved grassland"), outline=FALSE, ylim=c(-1,1),staplewex=0, border="white", frame=F,cex.axis=1.2,cex.lab=1.4, ylab="Fractional change")

points(cbind(jitter(rep(1,nrow(LU.GB)),1,0.4),LU.GB$chg.urb),col=transp(urb.col,0.2), pch=16,lwd=2)
points(cbind(jitter(rep(2,nrow(LU.GB)),1,0.4),LU.GB$chg.wdl),col=transp(wdl.col,0.2), pch=16,lwd=2)
points(cbind(jitter(rep(3,nrow(LU.GB)),1,0.4),LU.GB$chg.ara),col=transp(ara.col,0.2), pch=16,lwd=2)
points(cbind(jitter(rep(4,nrow(LU.GB)),1,0.4),LU.GB$chg.gra),col=transp(gra.col,0.2), pch=16,lwd=2)
points(cbind(jitter(rep(5,nrow(LU.GB)),1,0.4),LU.GB$chg.imp.gra),col=transp(imp.gra.col,0.2), pch=16,lwd=2)
abline(0,0, lty=5,lwd=3)

segments(seq(0.7,5.7,1),c(median(LU.GB$chg.urb),median(LU.GB$chg.wdl),median(LU.GB$chg.ara),median(LU.GB$chg.gra),median(LU.GB$chg.imp.gra)),seq(1.3,6.3,1), lwd=2.5)

points(1:5,c(mean(LU.GB$chg.urb),mean(LU.GB$chg.wdl),mean(LU.GB$chg.ara),mean(LU.GB$chg.gra),mean(LU.GB$chg.imp.gra)), pch=21, bg="white", cex=1.5, lwd=2)



## 2. Fate of lowland grassland and meadow
low.gra<-read.csv("GB_LandUseChange_LowlandGrasslandChange.csv")

# create summary table of fate of pixels classified as lowland meadow and pasture in historical maps. 
low.gra.perc<-as.matrix(100*sort(colSums(low.gra[,c("LCM.ara", "LCM.gra", "LCM.urb", "LCM.wdl", "LCM.wat", "LCM.imp.gra")]))/sum(low.gra[,c("LCM.ara", "LCM.gra", "LCM.urb", "LCM.wdl", "LCM.wat", "LCM.imp.gra")]))
low.gra.perc # Only 10% is currently in grassland category according to modern land cover maps.

# create barplot seen in Fig 1
barplot(low.gra.perc, horiz=FALSE, col=transp(c(wat.col,urb.col, wdl.col, gra.col,ara.col,imp.gra.col),0.8),legend.text = FALSE, axes=TRUE, cex.axis = 2.5)


## 3. Maps of land use change and historical land use
# Import data and packages
library(terra) # for working with raster data
library(RColorBrewer) # for colour pallettes used for land-use change figures
library(viridis) # for colour pallette used for land conversion figure
LU.GB.ras<-rast("GB_LandUseChange_Raster.tif" )

# Specify colour pallettes
pal.size=7
green.palette <- brewer.pal(n = pal.size, name = "Greens") # Green pallette for historical land use
pur.palette <- brewer.pal(n = pal.size, name = "PRGn") # Purple to green pallette for land use change
vir.fun<-colorRampPalette(viridis(pal.size)) # viridis pallette (green to dark blue) for land conversion
vir.palette<-vir.fun(pal.size)

# Plot
# Change in land use, seen in Fig 1.
plot(LU.GB.ras[[c("chg.ara", "chg.gra", "chg.urb", "chg.wdl", "chg.imp.gra")]],col=pur.palette, range=c(-1,1)) #

# Land conversion, seen in Fig 2.
plot(LU.GB.ras[["land.conv"]], col=vir.palette,range=c(0,1))

# Historical land cover (no figure in paper)
plot(LU.GB.ras[[c("lus.ara", "lus.gra", "lus.urb", "lus.wdl", "lus.wat")]],col=green.palette, range=c(0,1))