Peabody/Visuals.r
2026-09-14 12:39:42 -06:00

46 lines
4.7 KiB
R

library(tidyverse)
library(scales)
library(janitor)
library(paletteer)
dir.create("Results", showWarnings = FALSE)
source("Scripts/Load_IMPLAN.r")
source("Scripts/Load_REMI.r")
IMPLAN <- GET_IMPLAN_DATA()
REMI <- GET_REMI_DATA()
DATA <- full_join(IMPLAN,REMI)
REMI <- DATA %>% filter(Source=='REMI') %>% filter(Region=='WY') %>% group_by(Year,Region,Impact) %>% summarize(value=sum(value)) %>% filter(Region=='WY') %>% ungroup
ADJUST <- DATA %>% filter(Source=='IMPLAN') %>% group_by(Year,Impact,Region) %>% summarize(value=sum(value)) %>% ungroup %>% group_by(Year,Impact) %>% summarize(RATIO=max(value)/min(value)-1) %>% ungroup
ADJUST <- rbind(ADJUST,ADJUST %>% filter(Year==2028) %>% mutate(RATIO=0) %>% mutate(Year=2025),ADJUST %>% filter(Year==2028) %>% mutate(RATIO=0) %>% mutate(Year=2026),ADJUST %>% filter(Year==2028) %>% mutate(RATIO=0) %>% mutate(Year=2027))
REMI_DOLLAR <- REMI %>% filter(Impact!='Employment')
TEMP <- REMI_DOLLAR %>% filter(Year==2028,Impact=='Income') %>% mutate(value=0)
REMI_DOLLAR <- rbind(REMI_DOLLAR,rbind(TEMP %>% mutate(Year=2025),TEMP %>% mutate(Year=2026),TEMP %>% mutate(Year=2027)))
REMI_DOLLAR$Impact <- factor(REMI_DOLLAR$Impact,levels=c('Output','Value Added','Income'))
REMI_DOLLAR <- rbind(REMI_DOLLAR,REMI_DOLLAR %>% left_join(ADJUST) %>% mutate(Region='US',value=value*RATIO)%>% select(-RATIO))
REMI_DOLLAR$Region <- factor(ifelse(REMI_DOLLAR$Region=="WY",'Wyoming',"Rest of United States"),levels=rev(c("Wyoming","Rest of United States")))
########################Plot of total impact in US and in Wyoming
REMI_DOLLAR_PLOT <- ggplot(REMI_DOLLAR ,aes(x=Year,y=value/10^9,fill=Region,group=Region))+geom_area()+facet_wrap(~Impact,ncol=1)+ylab("Billion Dollars")+theme_bw()+ theme(plot.title = element_text(size = 22, face = "bold"), axis.title = element_text(size = 18),axis.text = element_text(size = 16),legend.text = element_text(size = 14),strip.text = element_text(size = 14),legend.position="top")+ theme(legend.title=element_blank())+scale_fill_manual(values=c("#3C3B6E","#FFC425"))+scale_x_continuous(breaks=seq(2025,2060,by=5))+scale_y_continuous(breaks=seq(0,5000,by=0.5))
ggsave( filename = "./Results/Dollar Values of Model.png", plot = REMI_DOLLAR_PLOT, width = 8.5, height = 11*(3./4), units = "in", dpi = 300)
#########################Post 2035 Values to show smaller scale
REMI_DOLLAR_PROD_PLOT <- ggplot(REMI_DOLLAR %>% filter(Year>2035) ,aes(x=Year,y=value/10^6,fill=Region,group=Region))+geom_area()+facet_wrap(~Impact,ncol=1)+ylab("Million Dollars")+theme_bw()+ theme(plot.title = element_text(size = 22, face = "bold"), axis.title = element_text(size = 18),axis.text = element_text(size = 16),legend.text = element_text(size = 14),strip.text = element_text(size = 14),legend.position="top")+ theme(legend.title=element_blank())+scale_fill_manual(values=c("#3C3B6E","#FFC425"))+scale_x_continuous(breaks=seq(2025,2060,by=5))+scale_y_continuous(breaks=seq(0,5000,by=100))
ggsave( filename = "./Results/Dollar Values of Model in Production.png", plot = REMI_DOLLAR_PROD_PLOT , width = 8.5, height = 11*(3./4), units = "in", dpi = 300)
#############Dollar values table
NPV <- REMI_DOLLAR %>% mutate(value=value/10^9,Discount_0=value,Discount_2=value/(1+0.02)^(Year-2028),Discount_5=value/(1+0.05)^(Year-2028),Discount_10=value/(1+0.1)^(Year-2028)) %>% group_by(Impact,Region) %>% summarize(Discount_0=sum(Discount_0),Discount_2=sum(Discount_2),Discount_5=sum(Discount_5),Discount_10=sum(Discount_10)) %>% ungroup
NPV[,3:6] <-round(NPV[,3:6],2)
write.csv(NPV,"./Results/NPV_Values.csv",row.names=FALSE)
############Employment
REMI_EMPLOY <- DATA %>% filter(Impact=='Employment',Source=='REMI') %>% group_by(Year,Region,Impact) %>% summarize(value=sum(value)) %>% ungroup
REMI_EMPLOY <- rbind(REMI_EMPLOY, REMI_EMPLOY %>% left_join(ADJUST) %>% mutate(value=value*RATIO,Region='US') %>% select(-RATIO))
REMI_EMPLOY$Region <- factor(ifelse(REMI_EMPLOY$Region=="WY",'Wyoming',"Rest of United States"),levels=rev(c("Wyoming","Rest of United States")))
REMI_EMPLOYMENT_PLOT <- ggplot(REMI_EMPLOY ,aes(x=Year,y=value,fill=Region,group=Region))+geom_area()+ylab("Employment (Full-Time Equivalent)")+theme_bw()+ theme(plot.title = element_text(size = 22, face = "bold"), axis.title = element_text(size = 18),axis.text = element_text(size = 16),legend.text = element_text(size = 14),strip.text = element_text(size = 14),legend.position="top")+ theme(legend.title=element_blank())+scale_fill_manual(values=c("#3C3B6E","#FFC425"))+scale_x_continuous(breaks=seq(2025,2060,by=5))+scale_y_continuous(breaks=c(seq(0,50000,by=1500),13000),labels = label_comma())
ggsave( filename = "./Results/Employment.png", plot = REMI_EMPLOYMENT_PLOT , width = 8.5, height = 11*(2/4), units = "in", dpi = 600)