pacman::p_load(dplyr,pastecs,data.table,lubridate,magrittr,ggplot2,tidyr,scales,patchwork,foreign,Metrics,stringr,readxl,sf,raster,RColorBrewer)

# Figure 3 ####
LF <- read_excel('***/Nonroad.xlsx',sheet = 'NRCMLF') %>% as.data.table()
LF$功率 <- LF$功率 %>% factor(levels = c('37kW 以下','37-75kW','75-130kW','130kW 以上'),
                          labels = c('Less than 37kW','37-75kW','75-130kW','More than 130kW'))
LF$机型 <- LF$机型 %>% factor(levels = c("挖掘机","推土机","装载机","叉车","压路机"),
                          labels = c("Excavator","Bulldozer","Loader","Forklift","Roller"))

load_F <- ggplot(LF)+
  geom_bar(aes(x = `机型`,y = `负载因子更新后`,fill = `功率`),stat = 'identity',position = 'dodge',color = 'white')+
  geom_hline(yintercept = 0.65,color = 'red')+
  scale_fill_manual(values = c('#3CAE93','#28C1B6','#C8E2BA','#BFBFBF'),name = 'Power Distribution')+
  xlab(NULL)+
  ylab('Load factor')+
  theme_bw()+
  theme(
    text = element_text(face = 'bold',family = 'sans'),
    axis.text = element_text(face = 'bold',family = 'WRYH'),
    axis.title = element_text(face = 'bold',family = 'WRYH'),
    legend.position = 'bottom',
    legend.text = element_text(face = 'bold',family = 'WRYH'),
    legend.title = element_text(face = 'bold',family = 'WRYH'))+ coord_cartesian(ylim = c(0.4, 0.88))
ggsave('***/Load_factor_eng.png',load_F,dpi = 600,height = 4.5,width = 7)

# Figure 4 ####
df_MB <- fread('***/gridded_emis.csv')
CMAQ_cn36_inter <- st_read('***/Intersect_36kmgrid_china_province_map.shp') # 36km grid intersect with china province shpfile 
for (i in 1:3) {
  vari <- list(c('道路移动源','非道路移动源'),'非道路移动源','道路移动源')[i] %>% unlist() %>% print()
  plot_name <- c('MB','nonroad','vehicle')[i]
  Change_MB <- df_MB[MB %in% vari & emis >0,.(emis = sum(emis)), by = c('FID_Grid','poll','Year')] %>% dcast(FID_Grid + poll ~Year, value.var = 'emis') # %>% print()
  Change_MB[is.na(Change_MB)] <- 0
  Change_MB[, E_1511 := (`2015`-`2011`) / 36 / 36 * 1000] # Emission Intensity  kgkm2yr-1
  Change_MB[, E_2015 := (`2020`-`2015`)/ 36 / 36 * 1000]  # Emission Intensity  kgkm2yr-1
  
  Change_MB_NOx <- Change_MB[poll == 'NOx' ,c('FID_Grid','E_1511','E_2015')] %>% melt(.,id.vars = c('FID_Grid'), variable.name = 'case')
  Change_MB_NOx[,re_value := case_when(value <= -800 ~ -800,
                                       value > -800 & value <= -300 ~ -300,
                                       value > -300 & value <= -100 ~ -100,
                                       value > -100 & value <= -50  ~ -50,
                                       value > -50  & value <= -10  ~ -10,
                                       value > -10  & value <= 0    ~ 0,
                                       value > 0    & value <= 10   ~ 10,
                                       value > 10   & value <= 100  ~ 100,
                                       value > 100   ~ 100)]
  # 拼接地理信息
  plot_MB <- left_join(Change_MB_NOx,CMAQ_cn36_inter,by = c('FID_Grid')) %>% as.data.table()
  plot_MB[,re_value := re_value %>% as.factor ]
  
  plot_MB_all <- plot_MB[LongName != "Taiwan",] %>% st_as_sf() 
  
  legend_color <- c("#08306B","#1664AB","#4A97C9","#93C4DE","#CFE1F2","#DFE701", "#F8B949", "#DD691E")
  
  pd <- ggplot()+
    geom_sf(data = plot_MB_all, aes(geometry = geometry,  fill = re_value),colour = NA) + 
    geom_sf(data = border_china_land, aes(geometry = geometry), fill = NA, color="black",size = 0.2)+
    geom_sf(data = Shp_shasha,aes(geometry = geometry))+
    geom_sf(data = nine,aes(geometry = geometry))+
    facet_wrap(vars(case))+
    scale_fill_manual(values = legend_color,name = '')+
    theme_bw()+
    coord_sf(crs = crs_84, ylim = c(16, 54),xlim = c(75, 135))+
    theme( 
      text = element_text(family = "sans",face = 'bold'),
      axis.title=element_text(face="bold",color="black"),
      axis.text = element_text(face="bold",color="black"),
      legend.background = element_blank(),
      legend.text = element_text(family = "sans",face = 'bold'),
      legend.title = element_text(face = "bold",family = 'WRYH'),
      legend.key.size = unit(0.4, 'cm'), 
      legend.key.width = unit(0.4, "cm"),
      plot.title = element_text(face="bold",hjust = 0.5,vjust = 0.1),
      strip.text = element_text(face = "bold",family = 'WRYH'))+ 
    guides(fill = guide_legend(ncol = 1, reverse = TRUE))
  # pd
  ggsave(paste0('***/NOx_change_',plot_name,'.jpg'),pd,width = 11, height = 10,dpi = 600)
  
  # 南海
  pd_Nanhai <- ggplot()+
    geom_sf(data = plot_MB_all, aes(geometry = geometry,  fill = re_value),colour = NA) +
    geom_sf(data = border_china_land, aes(geometry = geometry), fill = NA, color="black",size = 0.2)+
    geom_sf(data = Shp_shasha,aes(geometry = geometry))+
    geom_sf(data = nine,aes(geometry = geometry))+
    facet_wrap(vars(case))+
    scale_fill_manual(values =legend_color)+
    theme_bw()+
    theme( 
      text = element_text(family = "sans",face = 'bold'),
      axis.title=element_text(face="bold",color="black"),
      axis.text = element_text(face="bold",color="black"),
      legend.position = 'none',
      strip.text = element_blank())+
    coord_sf(crs = crs_84) + 
    scale_x_continuous(expand = c(0, 0), limits = c(107, 122), breaks = seq(70, 140, 10)) +
    scale_y_continuous(expand = c(0, 0), limits = c(2, 24), breaks = seq(10, 60, 10)) +
    guides(fill = "none", color = "none") +
    theme_bw() +
    theme(
      panel.grid = element_blank(),
      strip.background = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank()
    )
  # pd
  ggsave(paste0('**/NOx_change_',plot_name,'_nanhai.jpg'),pd_Nanhai,width = 4, height = 4,dpi = 600)
  
  # VOCs
  Change_MB_NOx <- Change_MB[poll == 'VOCs' ,c('FID_Grid','E_1511','E_2015')] %>% melt(.,id.vars = c('FID_Grid'), variable.name = 'case')
  # 确定离散画图的label 
  options(scipen = 200)
  Change_MB_NOx[,re_value := case_when(value <= -200 ~ -200,
                                       value > -200 & value <= -100 ~ -100,
                                       value > -100 & value <= -50  ~ -50,
                                       value > -50  & value <= -10  ~ -10,
                                       value > -10  & value <= -5   ~ -5,
                                       value > -5   & value <= 0    ~ 0,
                                       value > 0    & value <= 5    ~ 5,
                                       value > 5   & value <= 10    ~ 10,
                                       value > 10   ~ 10)]
  # 拼接地理信息
  plot_MB <- left_join(Change_MB_NOx,CMAQ_cn36_inter,by = c('FID_Grid')) %>% as.data.table()
  plot_MB[,re_value := re_value %>% as.factor ]
  plot_MB_all <- plot_MB[LongName != "Taiwan",] %>% st_as_sf() 
  
  pd <- ggplot()+
    geom_sf(data = plot_MB_all, aes(geometry = geometry,  fill = re_value),colour = NA) + 
    geom_sf(data = border_china_land, aes(geometry = geometry), fill = NA, color="black",size = 0.2)+
    geom_sf(data = Shp_shasha,aes(geometry = geometry))+
    geom_sf(data = nine,aes(geometry = geometry))+
    facet_wrap(vars(case))+
    scale_fill_manual(values = legend_color,name = '')+
    theme_bw()+
    coord_sf(crs = crs_84, ylim = c(16, 54),xlim = c(75, 135))+
    theme( 
      text = element_text(family = "sans",face = 'bold'),
      axis.title=element_text(face="bold",color="black"),
      axis.text = element_text(face="bold",color="black"),
      legend.background = element_blank(),
      legend.text = element_text(family = "sans",face = 'bold'),
      legend.title = element_text(face = "bold",family = 'WRYH'),
      legend.key.size = unit(0.4, 'cm'), 
      legend.key.width = unit(0.4, "cm"),
      plot.title = element_text(face="bold",hjust = 0.5,vjust = 0.1),
      strip.text = element_text(face = "bold",family = 'WRYH'))+ 
    guides(fill = guide_legend(ncol = 1, reverse = TRUE))
  # pd
  ggsave(paste0('***VOCs_change_',plot_name,'.jpg'),pd,width = 11, height = 10,dpi = 600)
  
  # 南海
  pd_Nanhai <- ggplot()+
    geom_sf(data = plot_MB_all, aes(geometry = geometry,  fill = re_value),colour = NA) +
    geom_sf(data = border_china_land, aes(geometry = geometry), fill = NA, color="black",size = 0.2)+
    geom_sf(data = Shp_shasha,aes(geometry = geometry))+
    geom_sf(data = nine,aes(geometry = geometry))+
    facet_wrap(vars(case))+
    scale_fill_manual(values = legend_color)+
    theme_bw()+
    theme( 
      text = element_text(family = "sans",face = 'bold'),
      axis.title=element_text(face="bold",color="black"),
      axis.text = element_text(face="bold",color="black"),
      legend.position = 'none',
      strip.text = element_blank())+
    coord_sf(crs = crs_84) +
    scale_x_continuous(expand = c(0, 0), limits = c(107, 122), breaks = seq(70, 140, 10)) +
    scale_y_continuous(expand = c(0, 0), limits = c(2, 24), breaks = seq(10, 60, 10)) +
    guides(fill = "none", color = "none") +
    theme_bw() +
    theme(
      panel.grid = element_blank(),
      strip.background = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank()
    )
  # pd
  ggsave(paste0('***/VOCs_change_',plot_name,'_nanhai.jpg'),pd_Nanhai,width = 4, height = 4,dpi = 600)
}
 # Figure 5 ####
MB_df <- fread('***/All_Mobile_source_emission.csv',encoding = 'UTF-8') # unit tone
MB_df_total$Plot_poll <- factor(MB_df_total$Plot_poll,
                                levels = c('TOGs','VOCs','IVOCs','NOx','CO','CO2','PM25','PM10','NH3'), 
                                labels = c('TOGs','VOCs','IVOCs','NOx','CO','CO2','PM25','PM10','NH3'))
MB_df_total$MB <- MB_df_total$MB %>% factor(.,levels = c('LSPV','MPV','HPV','LSDT','MDT','HDT','Motocycle','NRAM','NRCM','Small_engine','Aircraft','Railway_machine','Ship') %>% rev(),
                                            labels = c('LDV','MDV','HDV','LDT','MDT','HDT','MC','NRAM','NRCM','SGPM','Aircraft','RE','Ship') %>% rev())

# a <- c('#4894B5',"#3A53A4","#65B4D1","#89CD9B","#BEAD4F","#DD691E","#821518",'#B30D00','#D1F6FF','#BFBFBF','#51AD69','#264a5f','#F59467')
p_plot_re <- ggplot(MB_df_total[Plot_poll %in% c('VOCs','IVOCs','TOGs','NOx','CO2','CO','PM25','PM10','NH3') & MB %in% c(
  'LDV','MDV','HDV','LDT','MDT','HDT','NRAM','NRCM','SGPM','Aircraft','RE'),])+
  geom_bar(aes(x = Time,y = Emis,fill = MB), position = 'stack',stat = 'identity',width = 0.8)+
  scale_fill_manual(values = c("#1F78B4","#4894B5","#3A53A4","#DFE701","#89CD9B","#49A871","#8491B4FF",
                               "#F79D1EFF","#DD691E","#BC3C29FF","#FDAF91FF","#1A5354FF","#6A3D9A") %>% rev() ,name = NULL)+
  facet_wrap(Plot_poll~.,scales = 'free',ncol = 3)+
  scale_x_continuous(breaks = seq(2011,2020,3))+
  ylab('Emissions (CO2 Unit: Pg, Others: Tg)')+
  xlab('Time')+
  guides(fill=guide_legend(nrow = 2,byrow = T))+
  theme_bw()+
  theme(
    text = element_text(face = 'bold',family = 'sans'),
    legend.position = 'bottom',
    legend.text = element_text(face = 'bold',family = 'WRYH'),
    axis.title =  element_text(face = 'bold',family = 'WRYH'),
    legend.key.width = unit(0.5,'cm'),
    strip.text = element_text(face = 'bold'),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text = element_text(face = 'bold'))+guides(fill = guide_legend(nrow = 1, reverse = TRUE))
# p_plot_re
ggsave('***/Trend_poll.png',p_plot_re,width = 9,height =8,dpi = 600)

# Figure 6 ####
MB_df <- fread('***/All_Mobile_source.csv') # 单位吨
Normal_process_trend <- rbind(MB_df[MB == 'Vehicle' & Poll %in% c('PM25'),],MB_df[where %in% c('Tailpipe','Evaporation') & Poll %in% c('VOCs') & MB == 'Vehicle',]) 
Normal_process_trend[,Emis:= case_when(where == 'Tailpipe'   ~ Emis ,   # Gg
                                       TRUE ~ Emis * (-1))] # -Gg

Normal_process_trend$S_vtype <- Normal_process_trend$S_vtype %>% factor(.,levels = c("MC","HDV","MDV","LDV","HDT","MDT","LDT"),
                                                                        labels = c("MC","HDV","MDV","LDV","HDT","MDT","LDT"))

Theinvtype <- c('#1F78B4','#4894B5','#3A53A4','#DFE701','#89CD9B','#49A871','#8491B4')

p1 <-   ggplot(Normal_process_trend[Poll == 'VOCs' & where %in% c('Tailpipe' ,'Eva')])+
  geom_bar(aes(Time , Emis, fill = S_vtype),position = 'stack',stat = 'identity',width = 0.75)+
  geom_abline(intercept = 0,slope = 0)+
  scale_fill_manual(values = Theinvtype,name = NULL)+
  scale_x_continuous(breaks = seq(2011,2020,1))+
  scale_y_continuous(breaks = seq(-800,2050,400))+
  ylab('Emissions (Gg)')+
  xlab('Time')+
  theme_bw()+
  theme(
    text = element_text(face = 'bold',family = 'sans'),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = 'None',
    axis.text = element_text(face = 'bold'),
    axis.title = element_text(face = 'bold',family = 'WRYH'))
p2 <- ggplot(Normal_process_trend[Poll == 'PM25'])+
  geom_bar(aes(Time , Emis, fill = S_vtype),position = 'stack',stat = 'identity',width = 0.75)+
  geom_abline(intercept = 0,slope = 0)+
  scale_fill_manual(values = Theinvtype,name = NULL)+
  # scale_fill_manual(values = c(test_R(6,1)),name = NULL)+
  scale_x_continuous(breaks = seq(2011,2020,1))+
  scale_y_continuous(breaks = seq(-100,400,50))+
  guides(fill = guide_legend(nrow = 1))+
  ylab('
       ')+
  xlab('Time')+
  theme_bw()+
  theme(
    text = element_text(face = 'bold',family = 'sans'),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = 'bottom',
    legend.background = element_blank(),
    legend.title = element_text(face = 'bold',family = 'WRYH',hjust = 0.5),
    legend.text = element_text(face = 'bold'),
    axis.text = element_text(face = 'bold'),
    axis.title = element_text(face = 'bold',family = 'WRYH'))
output <- p1 + p2 
output
ggsave('***/Trend_process.png',output,width = 10,height = 4.5,dpi = 600)

# Figure 7 ####
ncl <- c('#CDECFB','#96D4F3','#6AADDC','#4894B5','#49A871','#74C14B','#CFDC57','#F9C34D','#F68233','#D72728','#921519')
test_ncl <- function(num,x){
  z_color <- ncl %>% colorRampPalette(.,bias = x)
  test_ncl <- z_color(num)
}

FIT_compare_city <- fread('***/FIT_compare_VP_city.csv') 
FIT_compare_county <- fread('***/FIT_compare_VP_county.csv') 

p1 <- ggplot(FIT_compare_county)+
  geom_point(aes(x = VP_county_stat  , y = VP_fit,color = factor(Time)))+
  geom_abline(slope = 1)+
  scale_colour_manual(values = test_ncl(10,1.7),name = '')+
  xlab('Statistical Value at the County level   ')+
  ylab('Simulated Value   ')+
  ggtitle(label = '(a)')+
  labs(col = '')+
  theme_bw()+
  theme(
    text = element_text(family = 'sans'),
    panel.grid.major = element_blank(),
    legend.position = 'none',
    plot.title = element_text(face="bold",hjust = 0),
    axis.title = element_text(face="bold",hjust = 0.5,family = 'WRYH'),
    axis.text = element_text(face="bold"),
  ) 
p2 <-ggplot(FIT_compare)+
  geom_point(aes(VP_City / 10000 ,FIT_PVTR / 10000 ,colour = as.factor(Time)),size = 0.7)+
  geom_abline(slope = 1,intercept = 0)+
  scale_colour_manual(values = test_ncl(10,1.7))+
  xlab('Statistical Value at the City level   ')+
  ylab(NULL)+
  labs(col = '')+
  ggtitle(label = '(b)')+
  guides(colour = guide_legend(ncol = 5,byrow = T))+
  theme_bw()+
  theme(
    text = element_text(family = 'sans'),
    legend.position = c(0.35,0.9),
    legend.direction = "horizontal",
    panel.grid.major = element_blank(),
    plot.title = element_text(face="bold",hjust = 0),
    axis.text = element_text(face="bold"),
    axis.title = element_text(face="bold",hjust = 0.5,family = 'WRYH'),
    legend.key.size = unit(0.6, 'cm'), #change legend key size
    legend.key.height = unit(0.3,'cm'), #change legend key height
    legend.key.width = unit(0.1, 'cm'), #change legend key width
    legend.background = element_blank()
  )
figure_VP_compare <- p1 + p2
ggsave('***/Verification_VP.png',figure_VP_compare,dpi = 600,width = 10,height = 5)

# Figure 8 ####
plot_total_vehinum <- fread('***/plot_total_vehinum.csv')
plot_total_vehinum$S_vtype <- plot_total_vehinum$S_vtype %>% factor(.,levels = c('LSPV','MPV','HPV','LSDT','MDT', 'HDT'),
                                                                    labels = c('LDV','MDV','HDV','LDT','MDT', 'HDT'))
p <- ggplot(plot_total_vehinum[S_vtype %in% c('LDV','MDV','HDV','LDT','MDT', 'HDT') & Time > 2010])+
  geom_point(aes(total_n /10000 , Fit_total /10000,  colour = as.factor(Time %>% as.character() %>% as.numeric())),size = .5)+
  geom_abline(slope = 1)+
  facet_wrap(vars(S_vtype), nrow = 2,scales = 'free')+
  scale_colour_manual(values = test_ncl(length(2011:2020),1.7), name = 'Time')+
  xlab('Statistical Value for total VP   ')+
  ylab('Simulated Value for total VP   ')+
  labs(col = 'Time')+
  guides(colour = guide_legend(nrow = 1,byrow = T))+
  theme_bw()+
  theme(text = element_text(family = 'sans',face="bold"),
        legend.position = 'bottom',
        panel.grid = element_blank(),
        legend.background = element_blank(),
        axis.title = element_text(face="bold",family = 'WRYH'),
        legend.title = element_text(face="bold",family = 'WRYH'), #change legend title font size
        legend.key.size = unit(0.1, 'cm'), 
        legend.text = element_text(face="bold",family = 'sans')
  ) #change legend text font size #"bold.italic"
p
ggsave('***/Comparee_Emissionstand_eng.png',p,width = 8,height =5.5,dpi = 600)

# Figure 9 ####
Com_df <- fread('***/Emis_compare_result.csv',encoding = 'UTF-8') 
my.colors <- function(x,b){
  hpwhite      <- rgb(255,255,255,maxColorValue=255)
  hpgrey        <- rgb(205,204,204,maxColorValue=255)
  hpdarkblue    <- rgb( 58, 83,164,maxColorValue=255)
  hpblue        <- rgb( 68,123,190,maxColorValue=255)
  hpturq        <- rgb(112,205,221,maxColorValue=255)
  hpgreen       <- rgb(154,205,112,maxColorValue=255)
  hpyellow      <- rgb(245,235, 21,maxColorValue=255)
  hporange      <- rgb(245,127, 32,maxColorValue=255)
  hpred         <- rgb(237, 34, 36,maxColorValue=255)
  hpdarkred     <- rgb(130, 21, 24,maxColorValue=255)
  Q.colors <- colorRampPalette(c(hpdarkblue,hpturq,hpgreen,hporange,hpdarkred),bias = b)
  my.colors <- Q.colors(x)
  return(my.colors)
}
p <- ggplot(Com_df) +
  geom_point(aes(x = This_study, y = Ref_Inv, color = Ref_plot)) +
  geom_abline(slope = 1) +
  facet_wrap(~Poll, scales = 'free') +
  scale_x_continuous(
    trans = log10_trans(),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  scale_y_continuous(
    trans = log10_trans(),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  xlab('This study (Tone)')+
  ylab('Others study (Tone)')+
  scale_color_manual(values = my.colors(Com_df$Ref %>% unique() %>% length(),1),name = NULL)+
  guides(colour = guide_legend(ncol = 6,byrow = T))+
  theme_bw()+
  theme(text = element_text(face = 'bold',family = 'sans'),
        axis.title = element_text(face = 'bold',family = 'WRYH'),
        strip.text = element_text(face = 'bold',family = 'sans'),
        axis.text = element_text(face = 'bold',family = 'sans'),
        legend.text = element_text(face = 'bold',family = 'sans'),
        legend.position = 'bottom',
        legend.title = element_blank(),
        panel.grid.minor = element_blank(),
        legend.key.height = unit(0.4,'cm'),
        legend.background = element_blank())
# p
ggsave('***/Emis_compare_SD.png',p,width = 10,height = 11,dpi = 600)