knitr::opts_chunk$set(echo=T, warning=F, comment=NA, message=F, eval=F)
Load rISIMIP package
library(rISIMIP)
# Choose one of totals, 5crops, 15crops lu_type <- "totals"
# Specify path of file directory filedir <- "/media/matt/Data/Documents/Wissenschaft/Data"
The rISIMIP package already includes landuse data for 1995 and 2080 of the totals landuse data, the code of how this data is retrieved can be found in the data-raw folder of the Github repository of this package.
data("landuse-totals_2005soc")
Combine all crop data
# Get summarised landuse data d <- data(package="rISIMIP") landuse <- d$results[,"Item"][grep(x=d$results[,"Item"], pattern="landuse-totals_")] # Read files into dataframe crops <- lapply(landuse, function(x){ data <- get(data(list=x)) data$year <- strsplit(strsplit(basename(x), split="_")[[1]][4], split="[.]")[[1]][1] data$scenario <- strsplit(basename(x), split="_")[[1]][1] data$model <- strsplit(basename(x), split="_")[[1]][2] return(data) }) crops <- do.call(plyr::rbind.fill, crops) #Calculate area of each cell in km2 data(landseamask_generic, package="rISIMIP") isimip_area <- raster::area(landseamask_generic, na.rm=TRUE) # km2 isimip_area <- as.data.frame(raster::rasterToPoints(isimip_area)) colnames(isimip_area) <- c("x", "y", "area") sum(isimip_area$area) # Add area to crops dataframe library(dplyr) crops <- left_join(crops, isimip_area) # Remove 2005soc data crops <- crops %>% filter(scenario != "2005soc")
Create summary bar plot
# Remove duplicate of 1995 data #crops <- crops %>% filter(year != 1995 | scenario != "rcp60") # Calculate remaining area not considered by land-use types crops$other <- 1 - (crops %>% dplyr::select(-c(x,y,area, year, scenario, model, matches("total"))) %>% rowSums(na.rm=TRUE)) # Calculate total area per category and turn into long format crops_sum <- crops %>% mutate_at(vars(-c(x,y,area,year,scenario,model)), funs(. * area)) %>% group_by(year, scenario, model) %>% summarise_at(vars(-c(x,y,area)), sum) %>% dplyr::select(-matches("total")) %>% filter(year < 2090) %>% tidyr::gather(var, value, -c(year, scenario, model)) %>% tidyr::unite(time_rcp, year, sep=" ")# Adapt unite accordingly # Calculate mean across models crops_mean <- crops_sum %>% group_by(time_rcp, var, scenario) %>% dplyr::summarise(mean=mean(value, na.rm=T)) crops_mean <- tidyr::drop_na(crops_mean) crops_mean$perc <- crops_mean$mean/sum(isimip_area$area)*100 library(ggplot2) colnames(crops_mean) <- c("time_rcp", "var", "scenario", "area", "perc") crops_mean$var <- factor(crops_mean$var) ggplot(data=crops_mean, aes(x=time_rcp, y=area/1000000, fill=var)) + geom_bar(stat="identity", position="stack") + facet_wrap(~scenario, ncol=1) + #geom_text(aes(label=round(perc, digits=2))) + scale_fill_discrete(name="Land-use type") + #scale_x_discrete(name="", labels=c("1995", "2050 \n RCP2.6", "2050 \n RCP6.0", "2080 \n RCP2.6", "2080 \n RCP6.0")) + scale_y_continuous(name="Area (km² x 1,000,000)", limits=c(0,150), expand=c(0,0)) + theme_bw() + theme(legend.title=element_text(face="bold"), axis.title.y=element_text(face="bold"), rect = element_rect(fill = "transparent", colour=NA))
Dominant land-use type
# Calculate dominant landuse per grid cell dominant_lu <- crops %>% filter(year != 1995 | scenario != "rcp60") %>% select(-c(area, matches("total"))) %>% tidyr::gather(var, value, -c(x,y,year,scenario, model)) %>% group_by(x,y,year,scenario, model) %>% summarise(lu=which.max(value)) %>% as.data.frame() labels_lu <- colnames(crops)[!colnames(crops) %in% c("x","y","year","scenario","model", "area")] labels_lu <- labels_lu[unique(dominant_lu$lu)] dominant_lu$lu <- factor(dominant_lu$lu, labels=labels_lu) dominant_lu <- dominant_lu %>% filter(year %in% c("1995", "2080")) %>% filter(scenario %in% c("rcp26", "rcp60")) %>% tidyr::unite(time_rcp, year, scenario, sep=" ") dominant_lu$time_rcp <- factor(dominant_lu$time_rcp, labels=c("1995", "2080 RCP2.6", "2080 RCP6.0")) # Create Map data(outline, package="ggmap2") ggplot() + geom_tile(data=dominant_lu, aes(x=x, y=y, fill=factor(lu))) + facet_grid(model ~ time_rcp) + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + scale_fill_discrete(name="Dominant landuse type", labels=labels_lu) + theme_bw() + theme(strip.background= element_blank()) + scale_x_continuous(name=expression(paste("Longitude (",degree,")")), expand=c(0.01,0.01), breaks=c(-120, -60, 0, 60, 120, 180)) + scale_y_continuous(name=expression(paste("Latitude (",degree,")")), expand=c(0.01,0.01), breaks=c(-60, -40, -20, 0, 20, 40, 60,80)) + coord_quickmap(xlim=c(-180,180), ylim=c(-60,80))
Land-use coverage
#Calculate mean across GCMs landuse <- crops %>% select(-c(area, cropland_total)) %>% tidyr::gather(var, value, -c(x,y,year,scenario,model)) %>% group_by(x,y,year,scenario, var) %>% summarise(value=mean(value,na.rm=T)) landuse$var <- factor(landuse$var) # Specify colour ramp colourtheme <- colorRampPalette(c("white", "#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))(255) # Create Map of 1995 landuse_1995 <- landuse[landuse$year == 1995 & landuse$scenario == "rcp26",] landuse_1995 <- landuse_1995[landuse_1995$var %in% c("cropland_irrigated", "cropland_rainfed", "pastures"),] ggplot() + geom_tile(data=landuse_1995, aes(x=x, y=y, fill=value*100)) + facet_wrap(var ~ ., ncol=1) + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + scale_fill_gradientn(name="% Cover", colours=colourtheme, na.value="transparent", limits=c(0,100)) + theme_bw() + theme(strip.background= element_blank()) + scale_x_continuous(name=expression(paste("Longitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-180, -90, 0, 90, 180)) + scale_y_continuous(name=expression(paste("Latitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-60, -40, -20, 0, 20, 40, 60, 80)) + coord_cartesian(xlim=c(-180,180), ylim=c(-60,84)) #Plot map for individual future years data("outline", package="ggmap2") landuse_rcp26 <- landuse[landuse$year == 2080 & landuse$scenario == "rcp26",] ggplot() + geom_tile(data=landuse_rcp26, aes(x=x, y=y, fill=value*100)) + facet_wrap(. ~ var) + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + scale_fill_gradientn(name="% Cover", colours=colourtheme, na.value="transparent", limits=c(0,100)) + theme_bw() + theme(strip.background= element_blank()) + scale_x_continuous(name=expression(paste("Longitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-180, -90, 0, 90, 180)) + scale_y_continuous(name=expression(paste("Latitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-60, -40, -20, 0, 20, 40, 60,80)) + coord_cartesian(xlim=c(-180,180), ylim=c(-60,84)) landuse_rcp60 <- landuse[landuse$year == 2080 & landuse$scenario == "rcp60",] ggplot() + geom_tile(data=landuse_rcp60, aes(x=x, y=y, fill=value*100)) + facet_wrap(. ~ var) + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + scale_fill_gradientn(name="% Cover", colours=colourtheme, na.value="transparent", limits=c(0,100)) + theme_bw() + theme(strip.background= element_blank()) + scale_x_continuous(name=expression(paste("Longitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-180, -90, 0, 90, 180)) + scale_y_continuous(name=expression(paste("Latitude (",degree,")")), expand=c(0.05,0.05), breaks=c(-60, -40, -20, 0, 20, 40, 60,80)) + coord_cartesian(xlim=c(-180,180), ylim=c(-60,84))
Land-use change
# Change land-use categories landuse_all <- crops %>% dplyr::select(-c(area, matches("total"))) landuse_all$biofuel_cropland <- landuse_all$biofuel_cropland_irrigated + landuse_all$biofuel_cropland_rainfed landuse_all$biofuel_cropland[is.na(landuse_all$biofuel_cropland)] <- 0 landuse_all$cropland <- landuse_all$cropland_rainfed + landuse_all$cropland_irrigated # Change format of land-use data landuse_fut <- landuse_all %>% select(-c(cropland_irrigated, cropland_rainfed, biofuel_cropland_irrigated, biofuel_cropland_rainfed)) %>% tidyr::gather(var, value, -c(x,y,year, scenario, model)) landuse_fut <- tidyr::spread(landuse_fut, year, value) landuse_fut <- landuse_fut %>% tidyr::replace_na(list(`1995`=0, `2080`=0)) # Calculate delta landuse delta_landuse <- landuse_fut %>% group_by(x,y,scenario,var) %>% dplyr::summarise_at(vars(`1995`:`2080`), mean, na.rm=TRUE) delta_landuse <- delta_landuse %>% mutate_at(vars(`2080`), funs(. - `1995`)) delta_landuse <- tidyr::gather(delta_landuse, year, value, -c(x, y, scenario, var)) delta_landuse$year <- as.numeric(delta_landuse$year) #Subset data for plotting landuse <- delta_landuse[delta_landuse$year == 2080,] landuse <- landuse %>% tidyr::unite(year, scenario, col="time_rcp") landuse$time_rcp <- factor(landuse$time_rcp, labels=c("2080 RCP2.6", "2080 RCP6.0")) landuse$var <- factor(landuse$var, labels=c("Biofuel cropland", "Cropland", "Other", "Pastures")) ggplot() + geom_tile(data=landuse, aes(x=x, y=y, fill=value*100)) + facet_grid(var ~ time_rcp, scale="free_y", switch="y") + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + scale_fill_gradientn(name="%", colours=rev(colorRampPalette( c("#00007F", "blue", "#007FFF", "cyan", "white", "yellow", "#FF7F00", "red", "#7F0000"))(255)), breaks=c(-90,-60,-30,0,30,60,90), values=c(1,0.7,0.52,0.5,0.48,0.3,0), limits=c(-90,90), na.value="transparent") + theme_classic() + theme(axis.title = element_blank(),axis.line = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(), panel.grid = element_blank(), strip.background= element_blank(), strip.placement = "outside", strip.text = element_text(size=10, face="bold"), rect = element_rect(fill = "transparent")) + coord_quickmap(xlim=c(-180,180), ylim=c(-60,85), expand=FALSE)
Land-use change in Europe
#' getData from the raster package downloads automatically GADM data per country deu <- raster::getData(name="GADM", country="DEU", level=1, path=paste0(filedir, "/GADM")) # Level1 gives regional boundaries #' Subset data by Bayern bavaria <- deu[deu$NAME_1 == "Bayern",] #Subset data for plotting landuse <- delta_landuse[delta_landuse$year == 2080,] landuse <- landuse %>% tidyr::unite(year, scenario, col="time_rcp") landuse$time_rcp <- factor(landuse$time_rcp, labels=c("2080 RCP2.6", "2080 RCP6.0")) landuse$var <- factor(landuse$var, labels=c("Biofuel cropland", "Cropland", "Other", "Pastures")) ggplot() + geom_tile(data=landuse, aes(x=x, y=y, fill=value*100)) + facet_grid(var ~ time_rcp, scale="free_y", switch="y") + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + geom_polygon(data=bavaria, aes(x=long, y=lat, group=group), fill="transparent", colour="black", lty="dashed") + scale_fill_gradientn(name="%", colours=rev(colorRampPalette( c("#00007F", "blue", "#007FFF", "cyan", "white", "yellow", "#FF7F00", "red", "#7F0000"))(255)), breaks=c(-100,-75,-50,-25,0,25,50,75,100), values=c(1,0.7,0.52,0.5,0.48,0.3,0), limits=c(-100,100), na.value="transparent") + theme_classic() + theme(axis.title = element_blank(),axis.line = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(), panel.grid = element_blank(), strip.background= element_blank(), strip.placement = "outside", strip.text = element_text(size=10, face="bold"), rect = element_rect(fill = "transparent")) + coord_quickmap(xlim=c(-11.5,40.5), ylim=c(35.5,70), expand=FALSE)
# Calculate threshold value (threshold <- delta_landuse %>% group_by(var, year) %>% filter(value > 0) %>% summarise(threshold = quantile(value, probs=0.75))) # Subset data by threshold top25_lu <- left_join(delta_landuse, threshold) top25_lu <- top25_lu %>% group_by(x,y,scenario, var, year) %>% filter(value > threshold) %>% select(-threshold) top25_lu$value <- 1 # Plot top25_lu %>% filter(year==2080) %>% ggplot(aes(x=x,y=y)) + geom_tile() + geom_polygon(data=outline, aes(x=long,y=lat, group=group), fill="transparent", colour="black") + facet_grid(var~scenario, switch="y") + labs(x="", y="") + theme_bw() + theme(strip.placement = "outside", strip.background = element_blank())
# List data files data_total <- list.files(paste0(filedir, "/ISIMIP2b/InputData/landuse"), pattern="totals", recursive=T, full.names=T) data_urban <- list.files(paste0(filedir, "/ISIMIP2b/InputData/landuse"), pattern="urbanareas", recursive=T, full.names=T) # Re-create isimip area data(landseamask_generic, package="rISIMIP") isimip_area <- raster::area(landseamask_generic, na.rm=TRUE) # km2 # Calculate total area of each land use class varnames <- c("cropland_total", "pastures", "cropland_irrigated", "cropland_rainfed") data_crop_areas <- lapply(varnames, function(x){ data_total_area_sum <- lapply(data_total, function(y){ data <- raster::stack(y, varname=x) data_total_area <- data*isimip_area data_total_area_sum <- as.data.frame(raster::cellStats(data_total_area, stat='sum', na.rm=TRUE)) colnames(data_total_area_sum)[1] <- "area" nc <- ncdf4::nc_open(y) timeref <- lubridate::year(strsplit(nc$dim[["time"]]$units, " ")[[1]][3]) if(ncdf4::ncvar_get(nc, nc$dim$time)[1] == 0){ data_total_area_sum$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) } else if (ncdf4::ncvar_get(nc, nc$dim$time)[1] != 0){ data_total_area_sum$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) - 1 } ncdf4::nc_close(nc) data_total_area_sum$scenario <- strsplit(basename(y), split="_")[[1]][1] data_total_area_sum$var <- x return(data_total_area_sum) }) data_crop_areas <- do.call("rbind", data_total_area_sum) return(data_crop_areas) }) data_crop_areas <- do.call("rbind", data_crop_areas) varnames <- c("biofuel_cropland_irrigated", "biofuel_cropland_rainfed") data_bio_crop_areas <- lapply(varnames, function(x){ data_total_area_sum <- lapply(data_total[c(4,6,7)], function(y){ data <- raster::stack(y, varname=x) data_total_area <- data*isimip_area data_total_area_sum <- as.data.frame(raster::cellStats(data_total_area, stat='sum', na.rm=TRUE)) colnames(data_total_area_sum)[1] <- "area" nc <- ncdf4::nc_open(y) timeref <- lubridate::year(strsplit(nc$dim[["time"]]$units, " ")[[1]][3]) if(ncdf4::ncvar_get(nc, nc$dim$time)[1] == 0){ data_total_area_sum$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) } else if (ncdf4::ncvar_get(nc, nc$dim$time)[1] != 0){ data_total_area_sum$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) - 1 } ncdf4::nc_close(nc) data_total_area_sum$scenario <- strsplit(basename(y), split="_")[[1]][1] data_total_area_sum$var <- x return(data_total_area_sum) }) data_crop_areas <- do.call("rbind", data_total_area_sum) return(data_crop_areas) }) data_bio_crop_areas <- do.call("rbind", data_bio_crop_areas) data_crop_areas <- dplyr::bind_rows(data_crop_areas, data_bio_crop_areas) # Turn into raster files data_urban_areas <- lapply(data_urban, function(x){ data <- raster::stack(x) data_urban_ts <- as.data.frame(raster::cellStats(data*isimip_area, stat='sum', na.rm=TRUE)) nc <- ncdf4::nc_open(x) timeref <- lubridate::year(strsplit(nc$dim[["time"]]$units, " ")[[1]][3]) if(ncdf4::ncvar_get(nc, nc$dim$time)[1] == 0){ data_urban_ts$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) } else if (ncdf4::ncvar_get(nc, nc$dim$time)[1] != 0){ data_urban_ts$year <- timeref + ncdf4::ncvar_get(nc, nc$dim$time) - 1 } ncdf4::nc_close(nc) data_urban_ts$var <- "urbanareas" colnames(data_urban_ts)[1] <- "area" data_urban_ts$scenario <- strsplit(basename(x), split="_")[[1]][1] return(data_urban_ts) }) data_urban_areas <- do.call("rbind", data_urban_areas)
Plot time series
data_ts <- dplyr::bind_rows(data_crop_areas, data_urban_areas) data_ts$scenario[data_ts$scenario == "2100rcp26soc"] <- "rcp26" data_ts$scenario[data_ts$scenario == "rcp26soc"] <- "rcp26" data_ts$scenario[data_ts$scenario == "rcp60soc"] <- "rcp60" data_ind <- data_ts[data_ts$var != "cropland_total",] ggplot(data = data_ind, aes(x = year, y = area/1000000, colour = factor(var), linetype=factor(scenario))) + scale_colour_discrete(name="Landuse type", labels=c("Biofuel cropland irrigated", "Biofuel cropland rainfed", "Cropland irrigated", "Cropland rainfed", "Pastures", "Urban areas")) + labs(x= "Year", y="Area (km² x 1,000,000)") + scale_linetype_discrete(name="Scenario") + geom_line() + theme_bw() + theme(strip.background= element_blank()) data_total <- data_ts[!data_ts$var %in% c("biofuel_cropland_irrigated", "biofuel_cropland_rainfed", "cropland_irrigated", "cropland_rainfed"),] ggplot(data = data_total, aes(x = year, y = area/1000000, colour = factor(var), linetype=factor(scenario))) + scale_colour_discrete(name="Landuse type", labels=c("Cropland total", "Pastures", "Urban areas")) + scale_linetype_discrete(name="Scenario") + labs(x= "Year", y="Area (km² x 1,000,000)") + geom_line() + theme_bw() + scale_x_continuous(breaks=c(1700,1800,1900,2000,2100,2200)) + theme(strip.background= element_blank())
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.