library(dplyr)
library(ggplot2)
library(tidyr)
library(scales)

#load baseline emission data
load('chc_ei_210521.Rdata')

# baseline emission summary
ei0 <- data.frame(summarise(group_by(ei,calendar_year,location,county,airbasin,vessel_type,vessel_type_s,engine_type,model_year,tier,horsepower_bin),
                average_hp = weighted.mean(average_hp,pop,na.rm=T),pop=sum(pop),tot_annual_hours=sum(tot_annual_hours),
                final_hp_hour=sum(final_hp_hour),fuel_exh_gpy=sum(fuel_exh_gpy),nox_tpd=sum(nox_tpd),pm10_tpd=sum(pm10_tpd),
                co_tpd=sum(co_tpd),hc_tpd=sum(hc_tpd),co2_tpd=sum(co2_tpd)))
ei0 <- transform(ei0,id=seq.int(nrow(ei0)),retrofit='no',age=calendar_year-model_year)

eil <- gather(ei,'fields','values',which(colnames(ei)=='pop'):which(colnames(ei)=='pm2_5_tpd'))
eil$scenario <- 'Base'
em_base_sum <- data.frame(summarise(group_by(eil,calendar_year,county,airbasin,vessel_type,vessel_type_s,engine_type,model_year,tier,horsepower_bin,average_hp,fields,scenario),
                                    values=sum(values)))

# ****************** Excursion rule: new excursion vessels become hybrid in 2025, 30% zero-emission  ****************** 
chc_excursion_electric <- ei0 %>% filter(model_year>=2025 & vessel_type=='Excursion') %>% 
  mutate(final_hp_hour = 0.3*final_hp_hour,pop = 0.3*pop,tier = 'Electric',model_year_n = model_year,low_hphr_pct=0)

# Excursion rule, 70% diesel emissions
ei01 <- mutate(ei0, final_hp_hour = ifelse(model_year>=2025 & vessel_type=='Excursion', 0.7*final_hp_hour, final_hp_hour),
               pop = ifelse(model_year>=2025 & vessel_type=='Excursion', 0.7*pop, pop))

# ***************** Add new low-use percentage, low-use vessels do not change in rules, low use emissions calculated *****************
low_use <- read.csv('low_use_pct_dca.csv')
ei1 <- merge(ei01,low_use,c('vessel_type','tier'),all.x = T)
ei1[is.na(ei1)] <- 0
# add placeholder for ei1
ei1 <- transform(ei1,model_year_n=model_year)

#1. Add new low-use percentages for all engines, emissions unchanged
low_use_em <- function(df){
  em_low <- gather(df,'fields','values',which(colnames(df)=='pop'):which(colnames(df)=='co2_tpd'))
  em_low <- mutate(em_low,values_low=values*low_hphr_pct,low_use_n=1)
  em_low1 <- mutate(em_low,values=values_low,values_low=NULL)
  # eng_low_fin <- data.frame(summarise(group_by(em_low,calendar_year,airbasin,vessel_type,engine_type,tier,horsepower_bin,low_use_n,fields),
  #                                    values=sum(values_low)))
  return(em_low1)
}
em_low_fin <- low_use_em(ei1)

# ***************** Short-run ferries go full zero emission *****************
chc_shortFerry <- subset(ei1, vessel_type == 'Ferry-Short Run')
short_ferry_rule <- function(df){
  df <- transform(df,model_year_n=model_year)
  # assign model year 0 for ZEV engine
  s1 <- df %>% filter(calendar_year>=2026) %>% mutate(model_year_n = calendar_year,tier='Electric')
  eng_fin <- rbind(s1,subset(df,!id %in% s1$id))
  return(eng_fin)
}
chc_shortFerry_rule <- short_ferry_rule(chc_shortFerry)

#****Commercial fishing only goes to Tier 3, not Tier 4. Tier 4 doesn't apply to commercial fishing.that's why it is assigned 2022 MY.
chc_comFish <- subset(ei1, vessel_type == 'Commercial Fishing')

comFish_rule <- function(df){
  c1 <- df %>% filter(model_year<=1987 & calendar_year >= 2030) %>% mutate(model_year_n = 2030, age = calendar_year-model_year_n)
  c2 <- df %>% filter(model_year %in% seq(1988,1997,1) & calendar_year >= 2031) %>% mutate(model_year_n = 2031, age = calendar_year-model_year_n)
  c3 <- df %>% filter(model_year %in% seq(1998,2006,1) & calendar_year >= 2032) %>% mutate(model_year_n = 2032, age = calendar_year-model_year_n)
  df_c <- rbind(c1,c2,c3)
  df_uc <- subset(df,!id %in% df_c$id)
  df_uc <- transform(df_uc, model_year_n =model_year) 
  df_fin <- rbind(df_c,df_uc)
  return(df_fin)
}
chc_comFish_rule <- comFish_rule(chc_comFish)

#********************* Update Tier 0 to 1 engines to Tier 3/4 according to the schedule, and then add DPF at a later time
chc_regular <- subset(ei1, !vessel_type %in% c('Commercial Fishing','Ferry-Short Run'))

t01_rule <- function(df){
  ## update tier 0/1 engines
  p1 <- df %>% filter(model_year<=1993 & calendar_year >= 2023) %>% mutate(model_year_n = 2023, age = calendar_year-model_year_n)
  p2 <- df %>% filter(model_year %in% seq(1994,2001,1) & calendar_year >= 2024) %>% mutate(model_year_n = 2024, age = calendar_year-model_year_n)
  p3 <- df %>% filter(model_year %in% seq(2002,2006,1) & calendar_year >= 2025) %>% mutate(model_year_n = 2025, age = calendar_year-model_year_n)
  chc_tier01_c <- rbind(p1,p2,p3) # Tier 0/1 that has been changed
  
  # add dpf
  dpf_year <- read.csv('case1_tier01dpf.csv')
  p1.1 <- chc_tier01_c %>% filter(vessel_type %in% dpf_year[dpf_year$rule_year==2029,'vessel_type'] & calendar_year >= 2029) %>% mutate(retrofit='yes')
  p1.2 <- chc_tier01_c %>% filter(vessel_type %in% dpf_year[dpf_year$rule_year==2030,'vessel_type'] & calendar_year >= 2030) %>% mutate(retrofit='yes')
  p1.3 <- chc_tier01_c %>% filter(vessel_type %in% dpf_year[dpf_year$rule_year==2031,'vessel_type'] & calendar_year >= 2031) %>% mutate(retrofit='yes')
  chc_tier01_c_dpf <- rbind(p1.1,p1.2,p1.3)
  chc_tier01_c_npdf <- subset(chc_tier01_c,!id %in% chc_tier01_c_dpf$id)
  chc_tier01_c1 <- rbind(chc_tier01_c_dpf,chc_tier01_c_npdf)
  
  # combine unchanged with changed
  chc_tier01_uc <- subset(df,tier %in% c(0,1) & !id %in% chc_tier01_c1$id)
  chc_tier01_fin <- rbind(chc_tier01_c1,chc_tier01_uc)
  return(chc_tier01_fin)
}
chc_t01 <- subset(chc_regular,tier %in% c(0,1))
chc_t01_rule <- t01_rule(chc_regular)

all.equal(sum(chc_t01$pop),sum(chc_t01_rule$pop)) # check

#************************** Tier 2/3/4 engines **************************
chc_t234_t3 <- subset(chc_regular,!tier %in% c(0,1) & horsepower_bin < 800)

#a. engines hp < 800 hp go to Tier 3+DPF
t3dpf_rule <- function(df){
  # case 2: Tier 3 and DPF
  c2_tier3dpf <- read.csv('case2_tier3dpf_new.csv')
  chc_c2 <- merge(df,c2_tier3dpf[,-4],c('vessel_type','horsepower_bin','model_year'))
  # Tier 2 go to Tier 3 and DPF
  chc_c2_t2 <- chc_c2 %>% filter(tier==2 & calendar_year>=rule_year) %>% mutate(model_year_n = rule_year, retrofit='yes')
  # Tier 3 add DPF
  chc_c2_t3 <- chc_c2 %>% filter(tier==3 & calendar_year>=rule_year) %>% mutate(retrofit='yes')
  chc_c2_changed <- rbind(chc_c2_t2,chc_c2_t3)
  chc_c2_fin <- rbind(chc_c2_changed,subset(chc_c2,! id %in% chc_c2_changed$id))
  return(chc_c2_fin)
}
chc_t234_t3rule <- t3dpf_rule(chc_t234_t3)

all.equal(sum(chc_t234_t3$pop),sum(chc_t234_t3rule$pop))

#b. engines hp>800 hp go to Tier 4/ Tier 4+DPF
chc_t234_t4 <- subset(chc_regular,!tier %in% c(0,1) & horsepower_bin >= 800)

t4dpf_rule <- function(df){
  c4_tier4dpf <- read.csv('case4_tier4dpf_new.csv')
  tier4dpf_pct <- read.csv('low_use_tier4_dca.csv')
  c4_tier4dpf <- merge(c4_tier4dpf,tier4dpf_pct,'vessel_type',all.x = T)

  chc_c4 <- merge(df,c4_tier4dpf[,-4],c('vessel_type','horsepower_bin','model_year'))
  # tier 2/3 go to tier 4 + DPF
  chc_c4_t23 <- chc_c4 %>% filter(tier %in% c(2,3) & calendar_year>=rule_year) %>% mutate(retrofit='yes',model_year_n=rule_year)
  # tier 4 add DPF
  chc_c4_t4 <- chc_c4 %>% filter(tier==4 & calendar_year>=rule_year) %>% mutate(retrofit='yes')
  chc_c4_changed <- rbind(chc_c4_t23,chc_c4_t4)
  chc_c4_fin <- rbind(chc_c4_changed,subset(chc_c4,!id %in% chc_c4_changed$id))
  return(chc_c4_fin)
}
chc_t234_t4rule <- t4dpf_rule(chc_t234_t4)

all.equal(sum(chc_t234_t4$pop),sum(chc_t234_t4rule$pop))

# *********************************** combine engines under rule ***********************************
chc_p1 <- rbind(chc_excursion_electric,chc_shortFerry_rule,chc_t01_rule,chc_t234_t3rule[,colnames(chc_t01_rule)]) #chc_comFish_rule separate calculation
colnames(chc_shortFerry_rule)
# tier4-dpf fraction is 1 for all those vessels
chc_p1$tier4_dpf_pct <- 1
# combine all cases
chc_regular_rule <- rbind(chc_p1,chc_t234_t4rule[,colnames(chc_p1)])
all.equal(sum(chc_regular_rule$final_hp_hour)+sum(chc_comFish_rule$final_hp_hour),sum(ei0$final_hp_hour))
str(chc_regular_rule)

# ************************************ Calculate emissions ************************************
# update tier and ef for 800 hp engines to tier 4
ef <- read.csv('C:/Users/liliu/OneDrive - California Air Resources Board/CHC/chc_model/input_emission_factors.csv')
ef_t <- ef %>% filter(horsepower_bin==800 & model_year>2022) %>% mutate(tier=4,CO=0.5039562,HC=0.03675035,NOx=1.042157,PM=0.02655537)
ef_new <- rbind(ef_t,subset(ef,!(horsepower_bin==800 & model_year>2022)))
fcf <- read.csv('C:/Users/liliu/OneDrive - California Air Resources Board/CHC/chc_model/input_fuel_correct.csv')
deterioation <- read.csv('C:/Users/liliu/OneDrive - California Air Resources Board/CHC/chc_model/input_deterioration.csv')

ef_new <- mutate(ef_new,model_year_n=model_year,model_year=NULL,engine_type_id=NULL)
ef <- mutate(ef,model_year_n=model_year,model_year=NULL,engine_type_id=NULL)
fcf <- mutate(fcf,model_year_n=model_year,model_year=NULL)

## update emission factors for Tier 4/3
cal_emissions <- function(df,EF){
  # if electric, fuel and emissions are 0
  df_electric <- subset(df, tier == 'Electric')
  df_electric <- mutate(df_electric,fuel_exh_gpy=0, nox_tpd=0, pm10_tpd=0, hc_tpd=0, co_tpd=0, co2_tpd=0)
  
  df_diesel <- subset(df, tier != 'Electric')
  chc_em_1 <- merge(df_diesel,select(EF,-tier),c('model_year_n','engine_type','horsepower_bin'),all.x = T)
  chc_em_1 <- merge(chc_em_1,fcf,c('model_year_n','calendar_year'))
  chc_em_1 <- merge(chc_em_1,select(deterioation,-c(engine_type_id,vessel_type_id)),c('vessel_type','engine_type','horsepower_bin'),all.x = T)
  # emission g/year to ton/day, calculate non-low-use emissions
  chc_em_2 <- mutate(chc_em_1,pop = pop*(1-low_hphr_pct),tot_annual_hours = tot_annual_hours*(1-low_hphr_pct),
                     final_hp_hour=final_hp_hour*(1-low_hphr_pct), fuel_exh_gpy=final_hp_hour*bsfc_lbsperhphr/7.0936325,
                     nox_tpd = final_hp_hour*NOx*nox_fcf*1.102*10^(-6)/365, pm10_tpd = final_hp_hour*PM*pm_fcf*1.102*10^(-6)/365,
                     hc_tpd = final_hp_hour*HC*hc_fcf*1.102*10^(-6)/365, co_tpd = final_hp_hour*CO*1.102*10^(-6)/365,
                     co2_tpd = fuel_exh_gpy*10.21*1.102/1000/365)
  # retrofit, add DPF penalty 4%
  chc_em_3 <- mutate(chc_em_2,pm10_tpd = ifelse(retrofit=='yes',pm10_tpd*(1-0.85*tier4_dpf_pct),pm10_tpd),
                     co2_tpd = ifelse(retrofit=='yes',co2_tpd*(1+0.04*tier4_dpf_pct),co2_tpd))
  # deterioration
  chc_em_4 <- mutate(chc_em_3,nox_tpd = ifelse(age<useful_life,nox_tpd*(1+nox_det*age),nox_tpd*(1+nox_det*useful_life)),
                     pm10_tpd = ifelse(age<useful_life,pm10_tpd*(1+pm_det*age),pm10_tpd*(1+pm_det*useful_life)),
                     hc_tpd = ifelse(age<useful_life,hc_tpd*(1+hc_det*age),hc_tpd*(1+hc_det*useful_life)),
                     co_tpd = ifelse(age<useful_life,co_tpd*(1+co_det*age),co_tpd*(1+co_det*useful_life)))

  # combine
  chc_all <- rbind(df_electric,chc_em_4[,colnames(df_electric)])
  chc_long <- gather(chc_all,'fields','values',which(colnames(chc_all)=='pop'):which(colnames(chc_all)=='co2_tpd'))
  chc_long$low_use_n=0
  return(chc_long)
}

em_regular_fin <- cal_emissions(chc_regular_rule,ef_new)
em_com_fish <- cal_emissions(chc_comFish_rule,ef)
# make column names the same
em_regular_fin <- mutate(em_regular_fin,tier4_dpf_pct=NULL)

#add back low-use vessel emissions
em_rule0 <- rbind(em_regular_fin,em_com_fish,em_low_fin)
sum(subset(em_rule0,fields=='final_hp_hour')$values)==sum(ei0$final_hp_hour)

# summarise at county level
em_rule0_sum <- data.frame(summarise(group_by(em_rule0,calendar_year,county,airbasin,vessel_type,vessel_type_s,engine_type,model_year_n,tier,horsepower_bin,average_hp,fields),
                            values = sum(values)))
colnames(em_rule0_sum)[which(colnames(em_rule0_sum) == 'model_year_n')] <- 'model_year'

#*********************************** apply compliance extension to rule***********************************
em_rule0_ves <- summarise(group_by(em_rule0_sum,vessel_type,calendar_year,fields),values=sum(values))
em_pre_ext <- subset(em_rule0_ves,fields %in% c('co_tpd','hc_tpd','nox_tpd','pm10_tpd'))

# add extension ratio
ext_ratio <- read.csv('compliance_extension.csv')
ext_template <- merge(em_pre_ext,ext_ratio,c('vessel_type','calendar_year'))

# apply a percentage reduction every year for each vessel type and field

comp_ext <- function(ves_type,field){
  dat <- ext_template[ext_template$vessel_type==ves_type & ext_template$fields==field,]
  lm_reg <- lm(values~calendar_year,data=dat)
  # em_new <- ext_template[0,]
  em_new <- data.frame(calendar_year=seq(dat[1,2],dat[2,2],1),values=0,vessel_type=ves_type,fields=field)
  em_new$values <- predict(lm_reg,em_new)
  return(em_new)
}

em_ext <- ext_template[0,]
for (ves_type in unique(ext_template$vessel_type)){
  for (field in unique(ext_template$fields)) {
    em_new <- comp_ext(ves_type,field)
    em_ext <- rbind(em_ext,em_new)
  }
}

scale_ratio <- merge(em_rule0_ves,em_ext,c('vessel_type','calendar_year','fields'))
scale_ratio <- mutate(scale_ratio,ratio = values.y/values.x)

# apply extension
em_rule0_ext <- merge(em_rule0_sum,scale_ratio[,c(1:3,6)],c('vessel_type','calendar_year','fields'),all.x = T)
em_rule0_ext[is.na(em_rule0_ext)] <- 1
em_rule0_ext <- mutate(em_rule0_ext,values_ext = values*ratio)
em_rule0_ext_fin <- mutate(em_rule0_ext,values=values_ext,values_ext=NULL,ratio=NULL)
em_rule0_fin <- mutate(em_rule0_ext,values_ext=NULL,ratio=NULL)

#*********************************  MORE RULEs: aux engines, renewable fuel, GHG/co2 ZEV penetration ********************************* 
# ** all engines need to use renewable diesel since 2023.** 
more_rules <- function(df,rule_name){
  #1. all aux engines 0.012 reduction in fuel,emissions
  df1 <- df %>% mutate(values= ifelse((engine_type == "aux" & calendar_year>=2023 & fields!='pop'), values*0.998,values))
  
  #2. Renewable diesel:NOx 10% PM 30% reduction for Tier 0-Tier 3 engines;No reduction for Tier 4, fuel use up 5%
  # add renewable fule penetration factors
  df1$renew_fraction <- 0
  df1[df1$calendar_year==2023,'renew_fraction'] <- 0.982
  df1[df1$calendar_year==2024,'renew_fraction'] <- 0.988
  df1[df1$calendar_year==2025,'renew_fraction'] <- 0.994
  df1[df1$calendar_year>=2026,'renew_fraction'] <- 1
  # calculate emissions
  df2 <- df1 %>% mutate(values= ifelse((tier != 4 & fields=='nox_tpd' & calendar_year>=2023), values*0.9*renew_fraction,values))
  df3 <- df2 %>% mutate(values= ifelse((fields=='pm10_tpd'& calendar_year>=2023), values*0.7*renew_fraction,values))
  df4 <- df3 %>% mutate(values= ifelse((fields=='fuel_exh_gpy'& calendar_year>=2023), values*1.05*renew_fraction, values))
  df5 <- df4 %>% mutate(values= ifelse((fields=='co_tpd'& calendar_year>=2023), values*0.9*renew_fraction, values))
  df_fin <- df5 %>% mutate(values= ifelse((fields=='hc_tpd'& calendar_year>=2023), values*0.95*renew_fraction, values))

  # 3. ZE penetration, GHG reduction only
  zev_input <- read.csv('zev_fraction.csv',check.names = F)
  zev_pct <- gather(zev_input,'vessel_type','zev_pct',2:6)
  eng_tot <- merge(df_fin,zev_pct,c('calendar_year','vessel_type'),all.x = T)
  eng_tot$zev_pct[is.na(eng_tot$zev_pct)] <- 0
  # apply ZE reductions to only CO2, since the ZE credits (7 year provision on existing engines would cause increase in existing emissions)
  em_zev <- mutate(eng_tot,values = ifelse(fields == 'co2_tpd',values*(1-zev_pct),values),scenario = rule_name,renew_fraction=NULL)
  return(em_zev)
}

#*********************************#*********************************#*********************************#*********************************
em_rule <- more_rules(em_rule0_fin,'Rule')
em_rule_ext <- more_rules(em_rule0_ext_fin,'Rule_Extension')

sum(em_base_sum[em_base_sum$fields=='co2_tpd','values'])
sum(em_rule[em_rule$fields=='co2_tpd','values'])

# ****************************************************************** comparison ******************************************************************
# compare emissions
em_rule_ext$zev_pct <- NULL
comp_em <- rbind(em_base_sum,em_rule_ext)
comp_tot <- data.frame(summarise(group_by(comp_em,scenario,calendar_year,fields),values=sum(values)))
comp_ab <- data.frame(summarise(group_by(comp_em,scenario,airbasin,calendar_year,fields),values=sum(values)))
comp_ves <- data.frame(summarise(group_by(comp_em1,scenario,calendar_year,vessel_type,fields),values=sum(values)))

ggplot(comp_tot,aes(calendar_year,values,color=scenario))+geom_line()+facet_wrap(~fields,scales='free_y')+expand_limits(y = 0)+theme_bw()+
  labs(y='Emissions (tpd)',x='Calendar Year',title='Statewide CHC emissions') + theme(legend.position = c(0.85, 0.2),legend.title = element_blank())
ggsave('comp_chc_emi_210628.tiff',dpi=600,width = 7,height=5)

ggplot(subset(comp_ves,fields=='nox_tpd'),aes(calendar_year,values,color=scenario))+geom_line()+facet_wrap(~vessel_type,scales='free_y')+
  expand_limits(y = 0)+theme_bw()+ labs(y='NOx Emissions (tpd)',x='Calendar Year',title='Statewide CHC NOx emissions by vessel type')+ 
  theme(legend.position = c(0.8, 0.1),legend.title = element_blank())
ggsave('comp_chc_ves_nox.tiff',width = 9,height=5.5)

ggplot(subset(comp_ves,fields=='pm10_tpd'),aes(calendar_year,values,color=scenario))+geom_line()+facet_wrap(~vessel_type,scales='free_y')+
  expand_limits(y = 0)+theme_bw()+ labs(y='PM10 Emissions (tpd)',x='Calendar Year',title='Statewide CHC PM emissions by vessel type')+ 
  theme(legend.position = c(0.8, 0.1),legend.title = element_blank())
ggsave('comp_chc_ves_pm.tiff',width = 9,height=5.5)


