# Prep employment by industry data from Censuses

# Bring in MSAs
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$msa = str_pad(msas$CBSA.Code,width = 5,pad = '0')
msas$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')

msas$st = str_pad(msas$FIPS.State.Code,width = 2,pad = '0')
msas$cnty = str_pad(msas$FIPS.County.Code, width = 3, pad = '0')
msas = filter(msas, mettype != '', st != '72')
msas = select(msas, st,cnty,mettype, msa)


# Bring in historical county crosswalk
county_hist = read_csv('raw_data/county_hist/eglp_county_crosswalk_endyr_2010.csv')

# Bring in 1970 industry data
ind70 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds98_1970_county.csv'))
ind70 = mutate(ind70, 
               cemp_public  = C09014, 
               cemp_manuf = C09004 + C09005, 
               cemp_prof = C09013, 
               cemp_fire = C09009,
               cemp_tot = C09001 + C09002 + C09003 + C09004 + C09005 + C09006 + C09007 +
                 + C09008 + C09009 + C09010 + C09011 + C09012 + C09013 + C09014 , 
               Year = 1970, 
               NHGISST = str_pad(STATEA,width = 3, side = 'right',pad = '0'),
               NHGISCTY = str_pad(COUNTYA, width = 4, side = 'right',pad = '0'))

ind70l = left_join(ind70, county_hist, by = c('Year','NHGISST','NHGISCTY'))
stopifnot(sum(is.na(ind70l$weight)) == 1) # Dillingham AK not merging, only has 500 employees
ind70l[ind70l$STATEA  == '02' & ind70l$COUNTYA == '070','NHGISST_2010'] = 20
ind70l[ind70l$STATEA  == '02' & ind70l$COUNTYA == '070','NHGISCTY_2010'] = '070'
ind70l[ind70l$STATEA  == '02' & ind70l$COUNTYA == '070','weight'] = 1
stopifnot(sum(is.na(ind70l$weight)) == 0) # Dillingham AK not merging, only has 500 employees
stopifnot(abs(sum(ind70l$weight) -nrow(ind70)) < 1)

ind70l = mutate(ind70l, 
                st = str_pad(NHGISST_2010,3,pad = '0') %>% substr(1,2),
                cnty = NHGISCTY_2010,
                cemp_public = weight * cemp_public,
                cemp_manuf = weight * cemp_manuf,
                cemp_prof = weight * cemp_prof,
                cemp_fire = weight * cemp_fire,
                cemp_tot = weight * cemp_tot,
                year = 1969)
ind70cl = ind70l %>% group_by(st,cnty,year) %>% summarise(across(contains('cemp'),~sum(.,na.rm = T)))
stopifnot(abs(sum(ind70cl$cemp_tot) -sum(ind70$cemp_tot) ) <1)



# 1980 data
ind80 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds107_1980_county.csv'))
ind80 = mutate(ind80, 
               cemp_public  = DIA015, 
               cemp_manuf = DIA003 + DIA004, 
               cemp_prof = DIA012 + DIA013 + DIA014, 
               cemp_fire = DIA009,
               cemp_tot = DIA001 + DIA002 + DIA003 + DIA004 + DIA005 + DIA006 + DIA007 +
                 + DIA008 + DIA009 + DIA010 + DIA011 + DIA012 + DIA013 + DIA014 + DIA015 ,
               NHGISST = str_pad(STATEA,width = 3, side = 'right',pad = '0'),
               NHGISCTY = str_pad(COUNTYA, width = 4, side = 'right',pad = '0'))

ind80l = left_join(ind80, county_hist %>% filter(Year == 1980), by = c('NHGISST','NHGISCTY'))
stopifnot(sum(is.na(ind80l$weight)) == 0) 
stopifnot(abs(sum(ind80l$weight) - nrow(ind80)) < 1)

ind80l = mutate(ind80l, 
                st = str_pad(NHGISST_2010,3,pad = '0') %>% substr(1,2),
                cnty = NHGISCTY_2010,
                cemp_public = weight * cemp_public,
                cemp_manuf = weight * cemp_manuf,
                cemp_prof = weight * cemp_prof,
                cemp_fire = weight * cemp_fire,
                cemp_tot = weight * cemp_tot,
                year = 1979)
ind80cl = ind80l %>% group_by(st,cnty,year) %>% summarise(across(contains('cemp'),~sum(.,na.rm = T)))
stopifnot(abs(sum(ind80cl$cemp_tot) - sum(ind80$cemp_tot) ) < 1)



# 1990 data
ind90 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds123_1990_county.csv'))
ind90 = mutate(ind90, 
               cemp_public  = E4P017, 
               cemp_manuf = E4P004 + E4P005, 
               cemp_prof = E4P016 + E4P015 + E4P014, 
               cemp_fire = E4P010,
               cemp_tot = E4P001 + E4P002 + E4P003 + E4P004 + E4P005 + E4P006 + E4P007 +
               + E4P008 + E4P009 + E4P010 + E4P011 + E4P012 + E4P013 + E4P014 + E4P015 +
                 + E4P016 + E4P017,
               NHGISST = str_pad(STATEA,width = 3, side = 'right',pad = '0'),
               NHGISCTY = str_pad(COUNTYA, width = 4, side = 'right',pad = '0'))

ind90l = left_join(ind90, county_hist %>% filter(Year == 1990), by = c('NHGISST','NHGISCTY'))
stopifnot(sum(is.na(ind90l$weight)) == 0) 
stopifnot(abs(sum(ind90l$weight) - nrow(ind90)) <1)

ind90l = mutate(ind90l, 
                st = str_pad(NHGISST_2010,3,pad = '0') %>% substr(1,2),
                cnty = NHGISCTY_2010,
                cemp_public = weight * cemp_public,
                cemp_manuf = weight * cemp_manuf,
                cemp_prof = weight * cemp_prof,
                cemp_fire = weight * cemp_fire,
                cemp_tot = weight * cemp_tot,
                year = 1989)
ind90cl = ind90l %>% group_by(st,cnty,year) %>% summarise(across(contains('cemp'),~sum(.,na.rm = T)))
stopifnot(abs(sum(ind90cl$cemp_tot) - sum(ind90$cemp_tot) ) < 1)



# 2000 data
ind00 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds151_2000_county.csv'))
ind00 = mutate(ind00, 
               cemp_public  = GMH026 + GMH013, 
               cemp_manuf = GMH003 + GMH016, 
               cemp_prof = GMH022 + GMH023 + GMH009 + GMH010, 
               cemp_fire = GMH008 + GMH021,
               cemp_tot = GMH001 + GMH002 + GMH003 + GMH004 + GMH005 + GMH006 + GMH007 +
                 + GMH008 + GMH009 + GMH010 + GMH011 + GMH012 + GMH013 + GMH014 + GMH015 +
                 + GMH016 + GMH017 + GMH018 + GMH019 + GMH020 + GMH021 + GMH022 + GMH023 +
                 + GMH024 + GMH025 + GMH026 ,
               NHGISST = str_pad(STATEA,width = 3, side = 'right',pad = '0'),
               NHGISCTY = str_pad(COUNTYA, width = 4, side = 'right',pad = '0'))

ind00l = left_join(ind00, county_hist %>% filter(Year == 2000), by = c('NHGISST','NHGISCTY'))
stopifnot(sum(is.na(ind00l$weight)) == 0) 
stopifnot(abs(sum(ind00l$weight) - nrow(ind00)) < 1)

ind00l = mutate(ind00l, 
                st = str_pad(NHGISST_2010,3,pad = '0') %>% substr(1,2),
                cnty = NHGISCTY_2010,
                cemp_public = weight * cemp_public,
                cemp_manuf = weight * cemp_manuf,
                cemp_prof = weight * cemp_prof,
                cemp_fire = weight * cemp_fire,
                cemp_tot = weight * cemp_tot,
                year = 1999)
ind00cl = ind00l %>% group_by(st,cnty,year) %>% summarise(across(contains('cemp'),~sum(.,na.rm = T)))
stopifnot(abs(sum(ind00cl$cemp_tot) - sum(ind00$cemp_tot) ) < 1)


# 2009 
ind07 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds195_20095_county.csv'))
ind07 = mutate(ind07, 
               cemp_public  = RP5E028 + RP5E055, 
               cemp_manuf = RP5E007 + RP5E034, 
               cemp_prof = RP5E044 + RP5E048 + RP5E021 + RP5E017, 
               cemp_fire = RP5E014 + RP5E041,
               cemp_tot = RP5E001, year = 2007) %>% select(st = STATEA, cnty = COUNTYA,
                                                           year, contains('cemp'))
 


# 2014 
ind12 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds206_20145_county.csv'))
ind12 = mutate(ind12, 
               cemp_public  = ABGTE028 + ABGTE055, 
               cemp_manuf = ABGTE007 + ABGTE034, 
               cemp_prof = ABGTE044 + ABGTE048 + ABGTE021 + ABGTE017, 
               cemp_fire = ABGTE014 + ABGTE041,
               cemp_tot = ABGTE001, year = 2012) %>% select(st = STATEA, cnty = COUNTYA,
                                                            year, contains('cemp'))

# 2019 
ind17 = read_csv(paste0('raw_data/nhgis_industries/',nhgis_ind_name,'_ds244_20195_county.csv'))
ind17 = mutate(ind17, 
               cemp_public  = ALZHE028 + ALZHE055, 
               cemp_manuf = ALZHE007 + ALZHE034, 
               cemp_prof = ALZHE044 + ALZHE048 + ALZHE021 + ALZHE017, 
               cemp_fire = ALZHE014 + ALZHE041,
               cemp_tot = ALZHE001, year = 2017) %>% select(st = STATEA, cnty = COUNTYA,
                                                            year, contains('cemp'))

# Now combine all of them and add on MSAs
inddta = bind_rows(ind70cl, ind80cl,ind90cl, ind00cl, ind07,ind12, ind17)

# All counties have positive employment, so don't have to worry about zeros
stopifnot(inddta$cemp_tot > 0)

inddta = left_join(inddta, msas,by = c('st','cnty'))

inddta = inddta %>% mutate(msafull = case_when(mettype == 'Metropolitan' ~ msa,
                                                      mettype == 'Micropolitan' ~ paste('nm',st,sep = '_'),
                                                      is.na(mettype) ~ paste('nm',st,sep = '_'))) %>% select(-msa)

# Now calculate metro percentages by msa type


  
  # Get totals by MSA and year
  drom = inddta %>% group_by(msafull,year) %>% summarise(
    cemp_public = sum(cemp_public, na.rm = T),
    cemp_manuf = sum(cemp_manuf, na.rm = T),
    cemp_prof = sum(cemp_prof, na.rm = T),
    cemp_fire = sum(cemp_fire, na.rm = T),
    cemp_tot = sum(cemp_tot, na.rm = T),
    
    cemp_pct_public = cemp_public / cemp_tot,
    cemp_pct_manuf = cemp_manuf / cemp_tot,
    cemp_pct_prof = cemp_prof / cemp_tot,
    cemp_pct_fire = cemp_fire / cemp_tot,
    
  )
  stopifnot(!is.na(drom$cemp_pct_public))
  
  stopifnot(drom$cemp_pct_public >= 0 & drom$cemp_pct_public <= 1)
  stopifnot(drom$cemp_pct_manuf >= 0 & drom$cemp_pct_manuf <= 1)
  stopifnot(drom$cemp_pct_prof >= 0 & drom$cemp_pct_prof <= 1)
  stopifnot(drom$cemp_pct_fire >= 0 & drom$cemp_pct_fire <= 1)
  
drom = drom %>% rename(msa = msafull)


write_csv(drom, 'intermediate_data/predictors/msa_census_emp_all.csv')


