Peabody/Visuals.r
2026-08-20 16:56:32 -06:00

218 lines
16 KiB
R

library(tidyverse)
library(janitor)
#install.packages("paletteer")
library(scales)
library(paletteer)
DATA_DIR <- 'Exports/'
OUTPUT_DIR <- './Results/Tables_and_Figures/'
dir.create(OUTPUT_DIR,recursive=TRUE,showWarnings=FALSE)
######################Employments
EMPLOY <- read_csv(paste0(DATA_DIR,"Direct, Indirect, and Induced -Employment.csv" ),skip=5) %>% pivot_longer(c(-Category,-Units),names_to="Year") %>% mutate(Year=parse_number(Year))
EMPLOY_NUM <- EMPLOY %>% filter(Units!="Proportion") %>% select(-Units,Employment=value) %>% filter(Year>=2027) %>% mutate(Employment=Employment*1000)
EMPLOY_NUM$Category <- gsub(" Employment","",EMPLOY_NUM$Category)
EMPLOY_NUM$Category <- factor(EMPLOY_NUM$Category,levels=rev(c("Direct","Indirect","Induced","Total")))
#EMPLOY_NUM %>% filter(Employment=='Direct')
#EMPLOY_NUM
#EMPLOY_NUM %>% filter(Category=='Total',Employment!=0)%>% mutate(Period=ifelse(Year>2033,"Operating","Construction")) %>% group_by(Period) %>% summarize(Employment=mean(Employment))
#EMPLOY_NUM %>% filter(Category=='Total')%>% mutate(Period=ifelse(Year>2033,"Operating","Construction")) %>% group_by(Period) %>% summarize(Employment=mean(Employment))
#EMPLOY_NUM %>% filter(Category=='Total') %>% print(n=100)
COLORS <- rev(paletteer_d("fishualize::Alosa_fallax",n=4,direction=1)[-2])
MAX_VAL <- round(max(EMPLOY_NUM$Employment))
#MAX_VAL
COOL_DOWN <- as.numeric(EMPLOY_NUM[5,3])
#EMPLOY_NUM
#EMPLOY_NUM %>% filter(Category=='Total') %>% arrange(Year)%>% print(n=100)
JOB_PLOT <-ggplot(EMPLOY_NUM %>% filter(Category!="Total"),aes(x=Year,y=Employment))+geom_area(position = "stack",aes(fill=Category,group=Category))+theme_bw()+scale_x_continuous(breaks=c(2027,seq(2030,2060,by=5)))+scale_fill_manual(values=COLORS)+geom_line(data=EMPLOY_NUM %>% filter(Category=="Total"),linewidth=1)+geom_text(data=df,aes(label = text),size=6, vjust = "inward", hjust = "inward")+theme(text = element_text(size = 20),legend.position = "top")+scale_y_continuous(labels = scales::comma,breaks=seq(-250,10000,by=500))
png(paste0(OUTPUT_DIR,"Job_plot.png" ) , units = "in", width = 11, height = 8, res = 600)
JOB_PLOT
dev.off()
#JOB_PLOT
###################Output graph
OUTPUT <- read_csv(paste0(DATA_DIR,"Direct, Indirect, and Induced -Output.csv" ),skip=5) %>% pivot_longer(c(-Category,-Units),names_to="Year") %>% mutate(Year=parse_number(Year))
OUTPUT_NUM <- OUTPUT %>% select(-Units,Output=value) %>% filter(Year>=2028) %>% mutate(Output=Output*1000)
OUTPUT_NUM
OUTPUT_NUM$Category <- gsub(" Output","",OUTPUT_NUM$Category)
OUTPUT_NUM$Category <- factor(OUTPUT_NUM$Category,levels=rev(c("Direct","Indirect","Induced","Total")))
OUTPUT_PLOT <- ggplot(OUTPUT_NUM %>% filter(Category!="Total"),aes(x=Year,y=Output))+geom_area(position = "stack",aes(fill=Category,group=Category))+theme_bw()+scale_x_continuous(breaks=c(2027,seq(2030,2060,by=5)))+scale_fill_manual(values=COLORS)+geom_line(data=OUTPUT_NUM%>% filter(Category=="Total"),linewidth=1)+theme(text = element_text(size = 20),legend.position = "top")+scale_y_continuous(labels = scales::comma,breaks=seq(0,2000,by=500))+ylab("Economic Output (Million USD)")
OUTPUT_PLOT
png(paste0(OUTPUT_DIR,"Output_plot.png"), units = "in", width = 10, height = 8, res = 600)
OUTPUT_PLOT
dev.off()
####GDP and Output
GDP <- read_csv(paste0(DATA_DIR,"Gross Domestic Product - By Region - GDP by Region.csv" ),skip=5) %>% pivot_longer(c(-Region,-Units),names_to="Year",values_to="GDP") %>% mutate(Year=parse_number(Year),GDP=GDP*1000) %>% select(-Units)
GDP <- GDP %>% filter(Region !='All Regions')
GDP$Region <- gsub(" County","",GDP$Region)
KEY_REGIONS <- GDP %>% group_by(Region) %>% summarize(GDP = median(GDP)) %>% arrange(desc(GDP)) %>% filter(GDP>0.5) %>% pull(Region)
KEY_REGIONS
OTHER_GDP <- GDP %>% filter(!(Region %in% KEY_REGIONS)) %>% group_by(Year) %>% summarize(Region="Other Counties",GDP=sum(GDP)) %>% ungroup
GDP <- rbind(GDP %>% filter(Region %in% KEY_REGIONS),OTHER_GDP)
GDP$Region <- factor(GDP$Region,levels=rev(c("Other Counties",KEY_REGIONS)))
GDP <- GDP %>% filter(Year>=2027)
sum(GDP$GDP)
OUTPUT_NUM %>% pull(Output) %>% sum
GDP <- GDP %>% group_by(Year) %>% summarize(GDP=sum(GDP)) %>% filter(Year>2027) %>% arrange(Year) %>% pull(GDP)
sum(GDP)
length(GDP)
sum(GDP*(1.09)^-(0:32))
GDP_SUB_REGION_PLOT <- ggplot(GDP %>% filter(Region!='Lincoln'),aes(x=Year,y=GDP,fill=Region))+geom_area(position='stack')+scale_x_continuous(breaks=c(2027,seq(2030,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,1000,by=100))+theme_bw()+theme(text = element_text(size = 20),legend.position = "top")+ylab("GDP (Million USD)")+scale_fill_manual(values=paletteer_d("PNWColors::Shuksan2"))
GDP_SUB_REGION_PLOT
png(paste0(OUTPUT_DIR,"GDP_Other_Counties.png"), units = "in", width = 10, height = 8, res = 600)
GDP_SUB_REGION_PLOT
dev.off()
GDP_REGION_SUMMARY <- GDP %>% mutate(Region=factor(ifelse(Region=='Lincoln','Lincoln','Rest of Wyoming'),levels=rev(c('Lincoln','Rest of Wyoming')))) %>% group_by(Region,Year) %>% summarize(GDP=sum(GDP)) %>% ungroup
GDP_PLOT <- ggplot(GDP_REGION_SUMMARY ,aes(x=Year,y=GDP,fill=Region))+geom_area(position='stack')+scale_x_continuous(breaks=c(2029,seq(2035,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,1000,by=25))+theme_bw()+theme(text = element_text(size = 20),legend.position = "top")+ylab("GDP (Million USD)")+scale_fill_manual(values=c(paletteer_d("PNWColors::Shuksan2")[c(5)],paletteer::paletteer_d("calecopal::lupinus")[2]))
#paletteer::paletteer_d("calecopal::lupinus")[3]
png(paste0(OUTPUT_DIR,"GDP_Region_Plot.png"), units = "in", width = 10, height = 8, res = 600)
GDP_PLOT
dev.off()
###################
INDUSTRY_JOBS <- read_csv(paste0(DATA_DIR,'Employment- By Industry - Employment by Industry.csv'),skip=5) %>% pivot_longer(c(-Industry,-Units),names_to="Year",values_to="Jobs") %>%mutate(Year=parse_number(Year)) %>% filter(Industry!='All Industries') %>% select(-Units) %>% filter(Year>=2029)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Data processing, hosting, and related services; Other information services',"Data processing",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Professional, scientific, and technical services',"Technical services",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Administrative and support services',"Administrative services",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='State and Local Government',"Government",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Food services and drinking places',"Restaurants",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Other transportation equipment manufacturing',"Equipment manufacturing",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Repair and maintenance',"Maintenance",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Utilities',"Energy Production",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS$Industry <- ifelse(INDUSTRY_JOBS$Industry=='Ambulatory health care services',"Hospitals",INDUSTRY_JOBS$Industry)
INDUSTRY_JOBS <- INDUSTRY_JOBS %>% group_by(Year) %>% mutate(Rank=rank(-Jobs)) %>% mutate(Rank=ifelse(Rank>10,11,Rank),Industry=ifelse(Rank==11,"Other",Industry)) %>% group_by(Year,Industry,Rank) %>% summarize(Jobs=sum(Jobs)) %>% ungroup(Industry,Rank) %>% arrange(Year,Rank) %>% ungroup
#INDUSTRY_JOBS %>% pull(Industry) %>% unique
ORDER <- c(INDUSTRY_JOBS %>% filter(Industry!='Other') %>% group_by(Industry) %>% summarize(MEAN=mean(Jobs)) %>% arrange(desc(MEAN)) %>% pull(Industry) %>% unique,"Other")
INDUSTRY_JOBS$Industry <- factor(INDUSTRY_JOBS$Industry,levels=ORDER)
#INDUSTRY_JOBS
JOB_TYPE_PLOT <- ggplot(INDUSTRY_JOBS ,aes(x=Year,y=Jobs,fill=Industry))+geom_bar(stat='identity')+ paletteer::scale_fill_paletteer_d("colorBlindness::Blue2DarkRed18Steps",name="")+theme_bw()+theme(text = element_text(size = 20),legend.position = "top")+guides(fill=guide_legend(nrow=4,byrow=TRUE)) +scale_x_continuous(breaks=c(2029,seq(2035,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,3000,by=500))
png(paste0(OUTPUT_DIR,"Job_Dist_Plot.png"), units = "in", width = 11, height = 8, res = 600)
JOB_TYPE_PLOT
dev.off()
JOB_FACET_DATA <- INDUSTRY_JOBS %>% mutate(Period=ifelse(Year>2033,"Operating","Construction")) %>% group_by(Period,Industry) %>% summarize(Jobs=mean(Jobs)) %>% ungroup
JOB_FACET <- ggplot(JOB_FACET_DATA,aes(x=Period,y=Jobs,fill=factor(Period)))+geom_bar(stat='identity',width=0.99)+facet_wrap(.~Industry,nrow=5)+theme_bw()+labs(x = NULL)+theme(text=element_text(size=20),legend.position="top",palette.colour.discrete=c("darkslategray1","darkorchid1" ),axis.text.x = element_blank())+scale_fill_discrete(name = "") +scale_y_continuous(labels = scales::comma,breaks=seq(0,1000,by=100))
png(paste0(OUTPUT_DIR,"Job_Facet_Plot.png"), units = "in", width = 11, height = 11, res = 600)
JOB_FACET
dev.off()
##################GDP
read_csv(paste0(DATA_DIR,'Gross Domestic Product - By Component - GDP Components.csv'),skip=5)
GDP_TOTAL <- read_csv(paste0(DATA_DIR,'Gross Domestic Product - By Component - GDP Components.csv'),skip=5)%>% select(-Units) %>% pivot_longer(c(-Category),names_to="Year") %>% mutate(value=value*1000,Year=parse_number(Year)) %>% filter(Year>=2029)
GDP_TOTAL <- GDP_TOTAL %>% filter(Category %in% c('Gross Domestic Product (GDP)','Consumption','Investment','Change in Private Inventories','Net Trade','Government Spending','Exogenous Final Demand'))
GDP_TOTAL <- GDP_TOTAL %>% filter(Category=='Gross Domestic Product (GDP)') %>% select(-Category) %>% rename(GDP=value)
GDP_PLOT <- ggplot(GDP_TOTAL ,aes(x=Year,y=GDP))+geom_line(linewidth=1,color='slateblue')+theme_bw()+theme(text=element_text(size=20),legend.position="top")+ylab("Wyoming GDP (Million USD)")+scale_x_continuous(breaks=c(seq(2030,2060,by=2)))+scale_y_continuous(breaks=seq(0,500,by=25))
png(paste0(OUTPUT_DIR,"GDP_Time_Plot.png"), units = "in", width = 11, height = 8, res = 600)
GDP_PLOT
dev.off()
################Personal Income
PERSONAL_INCOME <- read_csv(paste0(DATA_DIR,'Personal Income - By Region - Personal Income by Region.csv'),skip=5)%>% select(-Units) %>% pivot_longer(c(-Region),names_to="Year") %>% mutate(value=value*1000,Year=parse_number(Year)) %>% filter(Year>=2029) %>% rename(Income=value)
TOTAL_PERSONAL_INCOME <- PERSONAL_INCOME %>% filter(Region=='All Regions') %>% select(-Region)
#TOTAL_PERSONAL_INCOME
TAXES <- read_csv("Results/Tax_PI/Revenues.csv",skip=5) %>% clean_names() %>% filter(!is.na(revenue))
colnames(TAXES) <- gsub("fy","",colnames(TAXES))
TAXES <- TAXES %>% pivot_longer(c(-revenue,-units),names_to='year',values_to='Taxes') %>% mutate(Taxes=Taxes*1000) %>% select(-units)
TAXES$revenue <- gsub('State Sales & Use Taxes - ','Sales Tax: ', TAXES$revenue)
TAXES <- TAXES %>% filter(year>2028,!is.na(revenue))
TAXES <- TAXES %>% filter(Taxes>0)
TAXES$revenue[!grepl("Sales",TAXES$revenue)] <- 'Other Taxes'
TAXES <- TAXES %>% rename(Revenue='Taxes','Tax'=revenue,'Year'=year) %>% select(Year,Tax,Revenue)
TAXES <- TAXES %>% group_by(Year,Tax) %>% summarize(Revenue=sum(Revenue))
TOTAL_TAXES <- TAXES %>% group_by(Year) %>% summarize('Million (USD)'=sum(Revenue)/10^6,Metric='Tax Revenue') %>% mutate(Year=as.numeric(Year))
#TOTAL_TAXES
DOLLAR_VALUES <- rbind(OUTPUT_NUM %>% filter(Category=='Total') %>% select(-Category) %>% mutate(Metric='Economic Output') %>% rename('Million (USD)'=Output),
GDP_TOTAL %>% mutate(Metric='GDP') %>% rename('Million (USD)'=GDP),
TOTAL_TAXES,
TOTAL_PERSONAL_INCOME %>% mutate(Metric='Wages') %>% rename('Million (USD)'=Income))
FACET_INDICATORS_PLOT <- ggplot(DOLLAR_VALUES, aes(x=Year,y=`Million (USD)`,fill=Metric))+geom_area()+facet_wrap(.~Metric,nrow=2)+theme_bw()+scale_x_continuous(breaks=c(2029,seq(2035,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,3000,by=50))+theme(legend.position = "top",text=element_text(size=20))+scale_fill_manual(values=paletteer::paletteer_d("calecopal::lupinus") )
#DOLLAR_VALUES
DIRECT_TAX <- c(10.408,22.355,24.647,58.221,55.835,53.449,51.062,48.676,46.290,43.904,41.518,39.132,37.318,35.766,34.214,32.662,31.110,29.557,28.005,26.453,24.901,23.349,21.797,20.244,18.692,17.140,15.588,14.036,12.484,10.931,9.379)
DIRECT_TAX <- cbind(2030:2060,DIRECT_TAX ) %>% as_tibble
colnames(DIRECT_TAX) <- c("Year","Tax")
#DOLLAR_VALUES
DOLLAR_VALUES[DOLLAR_VALUES$Metric=='Tax Revenue','Million (USD)'] <- DOLLAR_VALUES[DOLLAR_VALUES$Metric=='Tax Revenue','Million (USD)'] + c(0,DIRECT_TAX$Tax)
TAX_PLOT <- ggplot(DOLLAR_VALUES %>% filter(Metric %in% c('Tax Revenue')), aes(x=Year,y=`Million (USD)`,fill=Metric))+geom_area()+theme_bw()+scale_x_continuous(breaks=c(2029,seq(2035,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,3000,by=5))+scale_fill_manual(values=paletteer::paletteer_d("calecopal::lupinus")[4],name="" )+theme(legend.position = "top",text=element_text(size=20))+ guides(fill = "none")
#TAX_PLOT
png(paste0(OUTPUT_DIR,"Total_Tax_Plot.png"), units = "in", width = 11, height = 8, res = 600)
TAX_PLOT
dev.off()
FACET_INDICATORS_PLOT <- ggplot(DOLLAR_VALUES , aes(x=Year,y=`Million (USD)`,fill=Metric))+geom_area()+facet_wrap(.~Metric,nrow=2)+theme_bw()+scale_x_continuous(breaks=c(2029,seq(2035,2060,by=5)))+scale_y_continuous(labels = scales::comma,breaks=seq(0,3000,by=50))+scale_fill_manual(values=paletteer::paletteer_d("calecopal::lupinus") )+theme(legend.position = "top",text=element_text(size=20))
#FACET_INDICATORS_PLOT
png(paste0(OUTPUT_DIR,"Facet_Indicator_Plot.png"), units = "in", width = 11, height = 10, res = 600)
FACET_INDICATORS_PLOT
dev.off()
#############
Period_Indicator_Summary <- rbind(OUTPUT_NUM %>% filter(Category=='Total') %>% mutate(Period=ifelse(Year<=2033,"Construction","Operation")) %>% group_by(Period) %>% summarize(Metric='Economic Output','Average'=mean(Output),Max=max(Output)),GDP %>% mutate(Period=ifelse(Year<=2033,"Construction","Operation")) %>% group_by(Period) %>% summarize(Metric='GDP',Average=mean(GDP),Max=max(GDP)),EMPLOY_NUM %>% mutate(Period=ifelse(Year<=2033,"Construction","Operation")) %>% group_by(Period) %>% summarize(Metric='Employment',Average=mean(Employment),Max=max(Employment)))
Period_Indicator_Summary[,c(3:4)] <- round(Period_Indicator_Summary[,c(3:4)],0)
#Period_Indicator_Summary
write_csv(Period_Indicator_Summary,paste0(OUTPUT_DIR,"Economic_Summary_Indicator.csv" ),col_names=TRUE)
GDP_SUMMARY <- GDP %>% group_by(Year) %>% summarize(GDP=sum(GDP),NPV_3=(GDP/((1+0.03)^(Year-2029))),NPV_5=(GDP/((1+0.05)^(Year-2029))),NPV_7=(GDP/((1+0.07)^(Year-2029)))) %>% ungroup
NPV_3 <- GDP_SUMMARY %>% pull(NPV_3) %>% sum
NPV_5 <- GDP_SUMMARY %>% pull(NPV_5) %>% sum
NPV_7 <- GDP_SUMMARY %>% pull(NPV_7) %>% sum
NPV_SUMMARY <- round(t(c(NPV_3,NPV_5,NPV_7)),0) %>% as_tibble
#NPV_SUMMARY
colnames(NPV_SUMMARY) <- c('3%','5%','7%')
write_csv(NPV_SUMMARY,paste0(OUTPUT_DIR,"GDP_Net_Present_Value.csv" ))
#######################
TAX <- DOLLAR_VALUES %>% filter(Metric=='Tax Revenue') %>% select(Year,Tax=`Million (USD)`)
#TAX %>% mutate(Period=ifelse(Year>2033,"Operating","Construction")) %>% group_by(Period) %>% summarize(Mean_Tax=mean(Tax),Max_Tax=max(Tax))
TAX_SUMMARY <- TAX %>% group_by(Year) %>% summarize(Tax=sum(Tax),NPV_3=(Tax/((1+0.03)^(Year-2029))),NPV_5=(Tax/((1+0.05)^(Year-2029))),NPV_7=(Tax/((1+0.07)^(Year-2029)))) %>% ungroup
NPV_3 <- TAX_SUMMARY %>% pull(NPV_3) %>% sum
NPV_5 <- TAX_SUMMARY %>% pull(NPV_5) %>% sum
NPV_7 <- TAX_SUMMARY %>% pull(NPV_7) %>% sum
NPV_TAX_SUMMARY <- round(t(c(NPV_3,NPV_5,NPV_7)),0) %>% as_tibble
#NPV_TAX_SUMMARY
colnames(NPV_TAX_SUMMARY) <- c('3%','5%','7%')
write_csv(NPV_TAX_SUMMARY,paste0(OUTPUT_DIR,"Tax_Net_Present_Value.csv" ))