world maps per capita

worldmaps.R

@@ -1,6 +1,7 @@
 ##########################################################################
 ##  Worldmaps
 ##########################################################################
+# install.packages("wbstats")
 library(ggplot2)
 library(sf)
 library(rnaturalearth)
@@ -16,7 +17,19 @@ agg <- function(x)
   x <- as.matrix(x) %*% sapply(unique(colnames(x)),"==",colnames(x))
   return(x)
 }
-
+per_capita <- function(countries, data, year = as.integer(format(Sys.Date(), "%Y"))-1){
+  pop_data <- wbstats::wb(indicator = "SP.POP.TOTL", startdate = year, enddate = year)
+  if(is.null(countries)) {
+    data_pc <- data.frame(country = as.character(data$country), 
+                          value = as.numeric(data$value))
+  } else {
+    data_pc <- data.frame(country = as.character(countries), 
+                          value = as.numeric(data))
+  }
+  data_pc$pop <- pop_data$value[match(data_pc$country, pop_data$iso3c)]
+  data_pc$value <- data_pc$value / data_pc$pop
+  return(data_pc[,1:2])
+}
 #----------------------------------------
 # read data
 #----------------------------------------
@@ -51,7 +64,7 @@ index_exio <- data.frame(country = rep(regions_exio$EXIOregion, each=200),
 # select year and allocation
 #----------------------------------------
 year <- 2013
-allocation <- c("mass","price")[2]
+allocation <- c("mass","price")[1]
 ######################################################
 
 #----------------------------------------
@@ -87,30 +100,45 @@ rm(L_fabio, L_link); gc()
 #----------------------------------------
 country <- c("BRA", "IDN")[1]
 product <- c("Cattle", "Soyabeans", "Oil, palm fruit", "Wood")[1]
+unit <- c("(kg per capita)", "(heads per capita)", "(cubic metres per capita)")[if_else(product=="Cattle", 2, if_else(product=="Wood", 3, 1))]
+percapita <- TRUE
 ######################################################
 
 #----------------------------------------
 # calculate footprints
 #----------------------------------------
-if(product=="Wood") element <- colSums(L[index_fabio$ISO == country & index_fabio$item %in% items$Item[items$Com.Group=="Wood"], ])
+if(product=="Wood") element <- colSums(L[index_fabio$ISO == country & index_fabio$item %in% items$Item[items$Com.Group=="Wood"], ]) / 1000
 if(product!="Wood") element <- L[index_fabio$ISO == country & index_fabio$item == product, ]
 
 data <- data.frame(region = colnames(Y),
                    data = colSums(element * Y))
-# share of ROW which is lost (i.e. not plotted)
-sum(data$data[192:197] / sum(data$data))
-world$footprint <- data$data[match(world$adm0_a3_is, data$region)] / sum(data$data)
-
-p <- ggplot(data = world) +
-  geom_sf(aes(fill = footprint), size = 0.05) +
-  labs(fill=paste0("Share of ",country," ",product)) + 
-  scale_fill_viridis_c(direction = -1, na.value = "lightgrey", labels = scales::percent) + 
-  theme(rect = element_blank()) + 
-  coord_sf(crs = "+proj=robin")
+if(percapita) {
+  data <- per_capita(data$region, data$data * 1000, year)
+  world$footprint <- data$value[match(world$adm0_a3_is, data$country)]
+  p <- ggplot(data = world) +
+    geom_sf(aes(fill = footprint), size = 0.05) +
+    labs(fill=paste0(country," ",product,"\n", unit)) + 
+    scale_fill_viridis_c(direction = -1, na.value = "lightgrey") + 
+    theme(rect = element_blank()) + 
+    coord_sf(crs = "+proj=robin")
   # coord_sf(crs = "+proj=moll")
   # coord_sf(crs = "+proj=wintri")
+} else {
+  # share of ROW which is lost (i.e. not plotted)
+  sum(data$data[192:197] / sum(data$data))
+  world$footprint <- data$data[match(world$adm0_a3_is, data$region)] / sum(data$data)
+  p <- ggplot(data = world) +
+    geom_sf(aes(fill = footprint), size = 0.05) +
+    labs(fill=paste0("Share of ",country," ",product)) + 
+    scale_fill_viridis_c(direction = -1, na.value = "lightgrey", labels = scales::percent) + 
+    theme(rect = element_blank()) + 
+    coord_sf(crs = "+proj=robin")
+  # coord_sf(crs = "+proj=moll")
+  # coord_sf(crs = "+proj=wintri")
+}
 
+p
 
-ggsave(filename = paste0("map_",country,"_",product,"_",allocation,".tif"), plot = p, device = "tiff", path = "./output",
-       scale = 1, width = 207, height = 90, units = "mm", dpi = 300)
+ggsave(filename = paste0("map_",country,"_",product,"_",if_else(percapita,"percapita_",""),allocation,".tif"), 
+       plot = p, device = "tiff", path = "./output", scale = 1, width = 207, height = 90, units = "mm", dpi = 300)