Plot_trade_loss_ExpRes.R

rm(list = ls())

# Libraries
library(tidyverse)
library(viridis)
library(patchwork)
library(hrbrthemes)
library(circlize)
library(chorddiag)  #devtools::install_github("mattflor/chorddiag")
library(reticulate)  # Load the reticulate package for working with Python objects
library(dplyr)
library(readr)

np <- import("numpy")

# Import label data
label_T <- read_csv("./Eora26_2016_bp/labels_T.csv",
                    col_types = cols(
                      Index = col_double(),
                      Country = col_character(),
                      ISO = col_character(),
                      Items = col_character(),
                      Sector = col_character()
                    )) %>%
  mutate(Country_Sector = paste(Country, Sector, sep = "_"))

label_T <- label_T %>%
  slice(1:4914)

name <- c(label_T$Country_Sector, unique(label_T$Country))

RR <- label_T %>%
  distinct(Country) %>%
  nrow()

NN <- label_T %>%
  distinct(Sector) %>%
  nrow()

RRNN <- RR * NN

years <- 2016:2050
TT <- length(years)


# Load no Export restriction order

expResSce <- c('Combined_5y', 'USA_5y', 'CHN_5y', 'IDN_5y', 'IND_5y')
expRes <- expResSce[1]

PATH1 <- "./Results/TechProExpResExtreme/temp/"

bs_cty <- np$load(paste0(PATH1, 'order_0_cty.npy'))
pExtreme <- np$load(paste0(PATH1, expRes, '_pExtreme_585_order.npy'))
nExtreme <- np$load(paste0(PATH1, expRes, '_nExtreme_585_order.npy'))
cExtreme <- np$load(paste0(PATH1, expRes, '_cExtreme_585_order.npy'))


# Total intermediate and final demand
combine_and_sum <- function(data) {
  aa = data[, 1:RR, ] + data[, (RR + 1):ncol(data), ]
  return(apply(aa, c(1, 2), sum))
}

# Intermediate order
inter <- function(data) {
  return(data[, 1:RR, ])
}

# calculate loss of flow
calculate_loss <- function(data_res, data) {
  return(data_res - data)
}

base_tot <- (bs_cty[, 1:RR] + bs_cty[, (RR+1):378]) * TT
pExtreme_tot <- combine_and_sum(pExtreme)
nExtreme_tot <- combine_and_sum(nExtreme)
cExtreme_tot <- combine_and_sum(cExtreme)


# calculate the diff between restriction and baseline
diff_pExtreme <- calculate_loss(pExtreme_tot, base_tot) / 1e6 # convert thousand to billion
diff_nExtreme <- calculate_loss(nExtreme_tot, base_tot) / 1e6
diff_cExtreme <- calculate_loss(cExtreme_tot, base_tot) / 1e6

# diag(diff_pExtreme) <- 0
# diag(diff_nExtreme) <- 0
# diag(diff_cExtreme) <- 0

diff_pExtreme <- as.data.frame(-diff_pExtreme)
diff_nExtreme <- as.data.frame(-diff_nExtreme)
diff_cExtreme <- as.data.frame(-diff_cExtreme)

rownames(diff_pExtreme) <- colnames(diff_pExtreme) <- unique(label_T$Country)
rownames(diff_nExtreme) <- colnames(diff_nExtreme) <- unique(label_T$Country)
rownames(diff_cExtreme) <- colnames(diff_cExtreme) <- unique(label_T$Country)

pExtreme_tot <- as.data.frame(pExtreme_tot)
nExtreme_tot <- as.data.frame(nExtreme_tot)
cExtreme_tot <- as.data.frame(cExtreme_tot)

rownames(pExtreme_tot) <- colnames(pExtreme_tot) <- unique(label_T$Country)
rownames(nExtreme_tot) <- colnames(nExtreme_tot) <- unique(label_T$Country)
rownames(cExtreme_tot) <- colnames(cExtreme_tot) <- unique(label_T$Country)


# 找到 'Dominican Republic' 在rownames和colnames中的索引
row_index <- which(rownames(diff_pExtreme) == 'Dominican Republic')
col_index <- which(colnames(diff_pExtreme) == 'Dominican Republic')

# 将 'Dominican Republic' 替换为 'Dominican Rep.'
rownames(diff_pExtreme)[row_index] <- colnames(diff_pExtreme)[col_index] <- 'Dominican'


# define data
data_long <- diff_pExtreme %>%
  rownames_to_column %>%
  gather(key = 'key', value = 'value', -rowname)

data_long <- data_long %>%
  arrange(desc(value)) %>%
  slice(1:50)

# order of country suffer most inflow and outflow loss
sum_by_key <- aggregate(value ~ key, data = data_long, sum, na.rm = TRUE)
sum_by_rowname <- aggregate(value ~ rowname, data = data_long, sum, na.rm = TRUE)
merged_data <- merge(sum_by_key, sum_by_rowname, by.x = "key", by.y = "rowname", all = TRUE) # value.x: inflow; value_y: outflow
merged_data[is.na(merged_data)] <- 0
final_sum <- aggregate(cbind(value.x, value.y) ~ key, data = merged_data, sum)
final_sum$total_sum <- final_sum$value.x + final_sum$value.y
final_sum <- final_sum[order(-final_sum$total_sum), ]
#print(final_sum[, c("key", "total_sum")])


# color palette
num_colors <- nrow(final_sum)
mycolor <- colorRampPalette(c("#B4DE2CFF", "#FDE725FF"))(num_colors)

country_colors <- data.frame(
  country = unique(final_sum$key),
  color = mycolor[1:length(unique(final_sum$key))]
)

# 创建特定国家与颜色的映射
specific_color_map <- c(
  "China" = "#FF0000",
  "Canada" = "#FF5733",
  "USA" = "#B22234",
  "Mexico" = "#006847",
  "South Korea" = "#003478",
  "Japan" = "#C32148",
  "Hong kong" = "#FFBE00",
  "Germany" = "#FFCC00",
  "UK" = "#00247D",
  "Ireland" = "#169B62",
  "Netherlands" = "#21468B",
  "Singapore" = "#ED2939",
  "Switzerland" = "#FF0000",
  "France" = "#0055A4"
)

# Apply country-specific colors to data
country_colors$color[country_colors$country %in% names(specific_color_map)] <- specific_color_map[country_colors$country[country_colors$country %in% names(specific_color_map)]]

order_cty <- final_sum$key

# parameters
circos.clear()
circos.par(start.degree = 90, 
           gap.degree = 1, 
           gap.after = 5, 
           track.margin = c(-0.15, 0.15), 
           points.overflow.warning = FALSE,
           canvas.xlim = c(-1.15, 1.15),  # horizontal, default is c(-1,1)
           canvas.ylim = c(-1.15, 1.15)   # vertical, , default is c(-1,1)
           )
par(mar = c(0, 0, 2, 0))  

PATH2 <- './FigExpRes_Extreme/'
png(file=paste0(PATH2, expRes, '_trade_loss_pExtreme_585.png'), height=2200, width=2200, bg = "#F4F1EB",res = 600)


# Base plot
chordDiagram(
  x = data_long, 
  grid.col = setNames(country_colors$color, country_colors$country),
  order = order_cty,
  transparency = 0.25,
  directional = 1,
  direction.type = c("arrows", "diffHeight"), 
  diffHeight  = -0.04,
  link.arr.type = "big.arrow",
  annotationTrack = "grid",
  annotationTrackHeight = c(0.05, 0.1)

)
# Add text and axis
circos.trackPlotRegion(
  track.index = 1, 
  bg.border = NA, 
  panel.fun = function(x, y) {
    
    xlim = get.cell.meta.data("xlim")
    sector.index = get.cell.meta.data("sector.index")
    
    
    for (i in seq_along(sector.index)) {
      label = sector.index[i]
      
      circos.text(
        x = mean(xlim), 
        y = 2,  # Adjusted y-coordinate to add spacing between label and plot
        labels = label, 
        facing = "clockwise",
        niceFacing = TRUE,
        adj = c(0, 0.5),
        cex = 0.5
      )
    }
    # Add graduation on axis
    circos.axis(
      h = "top", 
      major.at = seq(from = 0, to = xlim[2], by = 50), 
      minor.ticks = 0.1, 
      major.tick.length = 0.1,
      labels.niceFacing = TRUE,
      labels.cex = 0.3     # 设置刻度线数字的字体大小
    )
    par(cex.axis = 0.2)
    
    title(main = "Extreme Wet scenarios",cex.main = 0.7)
    #title(main = "Extreme Dry scenarios",cex.main = 0.7)
    #title(main = "Compound Extreme Anomalies",cex.main = 0.7)
  }
)
    

dev.off()