# Get PUMA to county allocations for 1970, 1980, 1990, 2000, 2010
# Match with 2018 MSA delineations
# Keeping one residual observation per state

dir.create('output/crosswalks')

#Bring in 2020 CBSA lookup
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$cbsa = str_pad(msas$CBSA.Code,width = 5,pad = '0')
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$fips = paste(msas$st,msas$cnty,sep = '')

msas = filter(msas, !is.na(st))

msas$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')
msas = mutate(msas, met2 = case_when(mettype == 'Micropolitan' ~ 'Non-metropolitan', mettype == 'Metropolitan' ~ 'Metropolitan'))

msas = select(msas,fips,cbsa,mettype,met2)


# Bring in 1970 County groups
# 1970 CGs only have one number set bc they cross state boundaries sometimes
# I split CGs by state so that I can calculate statistics by state - same process as with MSAs where you make more than 1 observation and weight them

cg7 = read_csv('intermediate_data/1970_crosswalk/xwalk_1970_cg_cnty_byrace.csv')
cg7$fips = paste(cg7$st, cg7$cnty,sep = '')

# Merge on MSAs
cg7msa = left_join(cg7,msas,by = 'fips')
cg7msa$mettype = replace_na(cg7msa$mettype,'Non-core')
cg7msa$met2 = replace_na(cg7msa$met2,'Non-metropolitan')

# Three different kinds of msa definition 
# CBSA + non-CBSA by state

cg7msa = cg7msa %>% mutate(
  
  msa = case_when(
    met2 == 'Metropolitan' ~ cbsa,
    met2 == 'Non-metropolitan' ~ paste('nm',st,sep = '_')
  )
  
)

# Note for 1970 have to do by race
pm7 = cg7msa %>% group_by(race,stpuma = cg,st,msa,met2) %>% summarise(
  afactor = sum(cwt),pop = sum(pop)
  )


# Bring in 1980 County groups
# Unlike PUMAs, county groups don't cross county boundaries
# But still need to make afactors in case only some CG counties are in a given MSA
cg8 = read.csv('raw_data/geocorr/cg98stat.csv')
cg8$st = str_pad(cg8$State.FIPS.Code,2,pad = '0')
cg8 = filter(cg8, st !=72)
cg8 = filter(cg8, !(st == 56 & County.FIPS.code == '.'))
cg8$cnty = str_pad(cg8$County.FIPS.code, 3, pad = '0')
cg8$fips = paste(cg8$st, cg8$cnty,sep = '')
cg8$cg = str_pad(cg8$X1980.County.Group,2,pad = '0')
cg8$stpuma = paste(cg8$st, cg8$cg,sep = '_')
cg8$cgpop = as.numeric(as.character(cg8$County.Group.Population))
cg8$cpop = as.numeric(as.character(cg8$County.Population))

cg8$afactor = cg8$cpop / cg8$cgpop
stopifnot(cg8$afactor <= 1)

# Merge on msas
cg8msa = left_join(cg8,msas,by = 'fips')
cg8msa$mettype = replace_na(cg8msa$mettype,'Non-core')
cg8msa$met2 = replace_na(cg8msa$met2,'Non-metropolitan')


cg8msa = cg8msa %>% mutate(
  msa = case_when(
    met2 == 'Metropolitan' ~ cbsa,
    met2 == 'Non-metropolitan' ~ paste('nm',st,sep = '_')
  )
  
)

pm8 = cg8msa %>% group_by(stpuma,st,msa,met2) %>% summarise(afactor = sum(afactor),pop = sum(cpop))


# Bring in 1990 PUMAs
puma9 = read.csv('raw_data/geocorr/geocorr1990_puma_to_county.csv')
puma9$st = str_pad(puma9$STATE.CODE..FIPS.,width = 2,pad = '0')
puma9$puma = str_pad(puma9$PUMA.code.from.1990.5...A..Sample,width = 5,pad = '0')
puma9$stpuma = paste(puma9$st,puma9$puma,sep = '_')
puma9$fips = str_pad(puma9$county, width = 5, pad = '0')

# Merge on msas
puma9msa = left_join(puma9,msas, by = 'fips')

puma9msa$mettype = replace_na(puma9msa$mettype,'Non-core')
puma9msa$met2 = replace_na(puma9msa$met2,'Non-metropolitan')


puma9msa = puma9msa %>% mutate(

  msa = case_when(
    met2 == 'Metropolitan' ~ cbsa,
    met2 == 'Non-metropolitan' ~ paste('nm',st,sep = '_')
  )
  
)

pm9 = puma9msa %>% group_by(stpuma,st,msa,met2) %>% summarise(afactor = sum(apuma.to.county.allocation.factor),pop = sum(X1990.POP))

sum(pm9$afactor > 1) #34 instances, but all just barely, looks like rounding
stopifnot(pm9$afactor <1.004)
pm9[pm9$afactor >1,'afactor'] = 1

## Now do 2000
puma0 = read.csv('raw_data/geocorr/geocorr2000_puma_to_county.csv',skip = 1)
puma0$st = str_pad(puma0$FIPS.state,width = 2,pad = '0')
puma0$puma = str_pad(puma0$puma5,width = 5,pad = '0')
puma0$stpuma = paste(puma0$st,puma0$puma,sep = '_')
puma0$fips = str_pad(puma0$county, width = 5, pad = '0')

puma0msa = left_join(puma0, msas, by = 'fips')

puma0msa$mettype = replace_na(puma0msa$mettype,'Non-core')
puma0msa$met2 = replace_na(puma0msa$met2,'Non-metropolitan')


puma0msa = puma0msa %>% mutate(
  msa = case_when(
    met2 == 'Metropolitan' ~ cbsa,
    met2 == 'Non-metropolitan' ~ paste('nm',st,sep = '_')
  )
  
)

pm0 = puma0msa %>% group_by(stpuma,st,msa,met2) %>% summarise(afactor = sum(puma5.to.county.allocation.factor),pop = sum(Total.Pop..2000.census))

# Check that afactor is <= 1
sum(pm0$afactor > 1) #54 instances, but all just barely, looks like rounding
stopifnot(pm0$afactor <1.003)
pm0[pm0$afactor >1,'afactor'] = 1


# Now 2010
puma1 = read.csv('raw_data/geocorr/geocorr2014_puma_to_county.csv',skip = 1)

puma1$st = str_pad(puma1$State.code,width = 2,pad = '0')
puma1$puma = str_pad(puma1$puma12,width = 5,pad = '0')
puma1$stpuma = paste(puma1$st,puma1$puma,sep = '_')
puma1$fips = str_pad(puma1$County.code, width = 5, pad = '0')

puma1msa = left_join(puma1, msas, by = 'fips')

puma1msa$mettype = replace_na(puma1msa$mettype,'Non-core')
puma1msa$met2 = replace_na(puma1msa$met2,'Non-metropolitan')


puma1msa = puma1msa %>% mutate(
  msa = case_when(
    met2 == 'Metropolitan' ~ cbsa,
    met2 == 'Non-metropolitan' ~ paste('nm',st,sep = '_')
  )
  
)

pm1 = puma1msa %>% group_by(stpuma,st,msa,met2) %>% summarise(afactor = sum(puma12.to.county.allocation.factor),pop = sum(Total.population..2010.))

sum(pm1$afactor > 1) #33 instances, but all just barely, looks like rounding
stopifnot(pm1$afactor <1.004)
pm1[pm1$afactor >1,'afactor'] = 1


#Check that populations sum
stopifnot(sum(pm7$pop) == sum(cg7$pop))
stopifnot(sum(pm8$pop) == sum(cg8$cpop))
stopifnot(sum(pm9$pop) == sum(puma9$X1990.POP))
stopifnot(sum(pm0$pop) == sum(puma0$Total.Pop..2000.census))
stopifnot(sum(pm1$pop) == sum(puma1$Total.population..2010.))


# Add years
pm7$year = 1970
pm8$year = 1980
pm9$year = 1990
pm0$year = 2000
pm1$year = 2010

# Add race for other years
pm8$race = 'all'
pm9$race = 'all'
pm0$race = 'all'
pm1$race = 'all'


#output/2021 crosswalks
write_csv(pm7, 'output/crosswalks/puma_msa_1970_byrace.csv')
write_csv(pm8, 'output/crosswalks/puma_msa_1980_byrace.csv')
write_csv(pm9, 'output/crosswalks/puma_msa_1990_byrace.csv')
write_csv(pm0, 'output/crosswalks/puma_msa_2000_byrace.csv')
write_csv(pm1, 'output/crosswalks/puma_msa_2010_byrace.csv')

pmall = rbind(pm7,pm8,pm9,pm0,pm1)
write_csv(pmall, 'output/crosswalks/puma_msa_allyrs_byrace.csv')

