# upload necessary packages
library(sf)
library(tmap)
library(raster)
library(pals)
library(dplyr)
library(tidyr)
library(ggplot2)
library(networkD3) # Sankey diagram
library(htmlwidgets)Prepare the landcover data
# identify home directory
normalizePath("~/")
lulc90 <- raster("~/R/Sankey diagram/Sankey diagram/c1990.tif")
lulc00 <- raster("~/R/Sankey diagram/Sankey diagram/c2000.tif")
lulc10 <- raster("~/R/Sankey diagram/Sankey diagram/c2010.tif")
lulc20 <- raster("~/R/Sankey diagram/Sankey diagram/c2020.tif")
#check that rasters are in the same dimensions ---------------------------------
compareRaster(lulc90, lulc00)# TRUE
compareRaster(lulc10, lulc20)# TRUE
# compareRaster(lulc00, lulc10)# FALSE
# compareRaster(lulc00, lulc20)# FALSE
# compareRaster(lulc90, lulc10)# FALSE
# compareRaster(lulc90, lulc20)# FALSE
extent(lulc90)
extent(lulc00)
extent(lulc10)
extent(lulc20)
extent(lulc90) <- extent(lulc20) # match extent, ignore the small differences in extent
extent(lulc00) <- extent(lulc20)
compareRaster(lulc90, lulc00, lulc10, lulc20) # TRUE
#summary of land cover data.
lulc90
lulc00
lulc10
lulc20
# assign land cover names -------------------------------------------------------
lulc <- ratify(lulc90) # access raster attribute
lulc <- levels(lulc)[[1]]
# ID: 0, 1, 2, 3, 4, 5, 6, 7 representing built-up, forest,
# herbaceous, bare land, water, cropland and shrubland, respectively
lulc$class <- c("Built-up", "Forest", "Herbaceous", "Bareland",
"Water", "Cropland", "Shrubland")
# assign new attribute (land use-land cover names)
levels(lulc90) <- lulc
levels(lulc00) <- lulc
levels(lulc10) <- lulc
levels(lulc20) <- lulcPrepare the nodes and connections
# create index
crosstab90_00 <- as.data.frame(crosstab(lulc90, lulc00, useNA=FALSE))
indexing_90 <- (rep(0:6, 7)) # for zero indexing of Sankey diagram
crosstab90_00$c1990 <- indexing_90
crosstab90_00 <- crosstab90_00[order(crosstab90_00$c1990),]
indexing_00 <- (rep(7:13, 7)) # indexing continued from indexing_90, so from number 7
crosstab90_00$c2000 <- indexing_00 #change numbering to indexing_00
#Now the same is done for 2000 and 2010
crosstab00_10 <- as.data.frame(crosstab(lulc00, lulc10, useNA=FALSE))
crosstab00_10$c2000 <- indexing_00 #change numbering to indexing_00
crosstab00_10 <- crosstab00_10[order(crosstab00_10$c2000),] #order them from small to big for easier indexing merge
indexing_10 <- (rep(14:20, 7)) #indexing for 2010
crosstab00_10$c2010 <- indexing_10 #include indexing for 2010
# 2010 and 2020
crosstab10_20 <- as.data.frame(crosstab(lulc10, lulc20, useNA = FALSE))
indexing_20 <- (rep(21:27, 7)) # indexing continued from indexing_10
crosstab10_20$c2020 <- indexing_20 #include indexing for 2020
crosstab10_20 <- crosstab10_20[order(crosstab10_20$c2010),]
crosstab10_20$c2010 <- indexing_10 #change numbering to indexing_10
#rename columns before merging so that both crosstabs have the same names
names(crosstab90_00)[1] = "Source"
names(crosstab90_00)[2] = "Target"
names(crosstab90_00)[3] = "Value"
names(crosstab00_10)[1] = "Source"
names(crosstab00_10)[2] = "Target"
names(crosstab00_10)[3] = "Value"
names(crosstab10_20)[1] = "Source"
names(crosstab10_20)[2] = "Target"
names(crosstab10_20)[3] = "Value"
# view the crosstabs
# dplyr::glimpse(crosstab90_00)
# dplyr::glimpse(crosstab00_10)
# dplyr::glimpse(crosstab10_20)
crosstab_total <- rbind(crosstab90_00, crosstab00_10, crosstab10_20)
# prepare the nodes from 1-28 for the sankey diagram.
# this needs to be a table with one column including the numbers 1 to 21 and a
# second column with each number linking to a land cover type and a year.
# these nodes will become the labels for the sankey diagram
nodes <- c(1:28)
#second column needs to be a combination of landcover and year.
# First copy the landcover types four times
categories <- data.frame(levels(lulc90)) |>
subset(select = -c(ID) )
categories4 <- rbind(categories, categories, categories, categories)
#I then made a list of the years 7 times (for each landcover type)
Year90 <- list(rep(1990, 7))
Year00 <- list(rep(2000, 7))
Year10 <- list(rep(2010, 7))
Year20 <- list(rep(2020, 7))
Year <- rbind.data.frame(Year90, Year00, Year10, Year20)
#the nodes are now combined in three columns. 1: number 1:21, 2: landcover types, 3: years
landcovernodes <- data.frame(nodes, categories4, Year)
#column 2 and 3 are pasted into one, so that landcover and year are in one column.
#the other columns are removed
landcovernodes$name <- paste(landcovernodes$class,
landcovernodes$c.1990..1990..1990..1990..1990..1990..1990..2000..2000..2000..)
landcovernodes <- dplyr::select(landcovernodes, -c(2, 3))
#To ensure the colors correspond with the land cover maps, I made a color file from number 1:7, corresponding with the nodes
my_color <- 'd3.scaleOrdinal() .domain(["1", "2", "3", "4","5","6","7"]) .range(["red", "green", "greenyellow", "gray", "blue", "gold", "olive"])'
#Include a grouping variable to colour the different flows in the diagram
groupa <- list(rep("1", 7))
groupb <- list(rep("2", 7))
groupc <- list(rep("3", 7))
groupd <- list(rep("4", 7))
groupe <- list(rep("5", 7))
groupf <- list(rep("6", 7))
groupg <- list(rep("7", 7))
groups <- rbind.data.frame(groupa, groupb, groupc, groupd, groupe, groupf, groupg)
groups_triple <- rbind.data.frame(groups, groups, groups)
names(groups_triple)[1] <- "group"
crosstab_total$group <- (groups_triple$group)Create the Sankey diagram
#the sankey diagram is made
sankeydiagram <- sankeyNetwork(Links = crosstab_total, Nodes = landcovernodes,
Source = "Source",Target = "Target", Value = "Value",
NodeID = "name", LinkGroup = "group",
colourScale = my_color, fontSize = 25,
nodeWidth = 12,
iterations=0,
margin = list("left" = 100))
# sankeydiagram <- onRender(sankeydiagram,
# '
# function(el, x){
# d3.select(el)
# .selectAll(".node text")
# .attr("text-anchor", "middle")
# .style("writing-mode", "horizontal-rl")
# .style("text-orientation", "upright");
# }
# '
# )
sankeydiagram