# Clean maps for text
# Median fam inc in each year with no labels
# Change 1969-2017


dir.create('output/maps')
dir.create('output/maps/outcomes')


# Download county shapes
cntys = counties(resolution = '20m', cb = T) %>% shift_geometry()

# Drop PR
cntys = filter(cntys, STATEFP != '72') #STATEFP != '02', STATEFP != '15',  # Still works here
cntys$fips = paste(cntys$STATEFP,cntys$COUNTYFP,sep = '')
cntys = rename(cntys, st = STATEFP)


# Get crosswalk from counties to MSAs
# Only including metros not micros
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$msa = 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$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')
msas = select(msas,fips,mettype,msa)

# Merge
cntys = left_join(cntys, msas, by = 'fips') 

# Assign leftover categories
cntys = cntys %>% mutate(msafull = case_when(mettype == 'Metropolitan' ~ msa,
                                          mettype == 'Micropolitan' ~ paste('nm',st,sep = '_'),
                                          is.na(mettype) ~ paste('nm',st,sep = '_'))) %>% select(-msa) %>% rename(msa = msafull)


# Bring in all the output and all the control variables
regout = read_csv('intermediate_data/regression_data/msa_reg_outcomes.csv')
outvars = unique(regout$incsamp)
rdefs = unique(regout$racedef)

popcut = 10000

rcuts = c(-1,.5,.6,.7,10000)
rlabs = c('0-49%','50-59%','60-69%','70% or greater')
regout$bw_rat_med_cats = regout$med_rat %>% cut(rcuts,labels = rlabs ,include.lowest = T)
regout$bw_rat_med_cats_popcut = regout$bw_rat_med_cats
regout[regout$pop_b < popcut, 'bw_rat_med_cats_popcut'] = NA


rkcuts = c(-1,.15,.25,.35,10000)
rklabs = c('Below 15th','At least 15th, below 25th','At least 25th, below 35th','At least 35th')
regout$rgap_cats = regout$rgap %>% cut(rkcuts,labels = rklabs ,include.lowest = T)
regout$rgap_cats_popcut = regout$rgap_cats
regout[regout$pop_b < popcut, 'rgap_cats_popcut'] = NA



for(incs in c('all_faminc_sqrt','emp_emp_inctot','gq_hh_spmuinc_disp_sqrt','gq_hh_spmuinc_sqrt')){
regint = filter(regout, incsamp == incs)


# Make yearly maps
# Dissolve into MSAs

msashape = cntys %>% group_by(msa) %>% summarise()  # Dissolve takes a long time to do


# Loop through rdefs
for(rdef in c('expansive','narrow')){
    
# Loop through years
  for(yr in c(1969,1979,1989,1999,2007,2012,2017)){

    # Get data
    dtaout = filter(regint, year == yr,  racedef == rdef)
    mdta = left_join(msashape, dtaout, by = 'msa')
    
    if(incs %in% c('gq_hh_spmuinc_disp_sqrt','gq_hh_spmuinc_sqrt') & yr < 2010) next()
    
    # Make map with pop cutoff
    pdf(paste('output/maps/outcomes/map_medrat_',incs,'_',yr,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
    gplt = ggplot(mdta, aes(fill = bw_rat_med_cats_popcut)) +
      geom_sf(size = .1, color = 'black') + 
      scale_fill_brewer(name = 'Black-white median family income ratio\n',
                        palette = 'YlGnBu', 
                        labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) 
    dev.off()
    
    if(rdef == 'expansive' & incs == 'all_faminc_sqrt'){
      pdf(paste('output/maps/map_appendix_medrat_',incs,'_',yr,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
      print(gplt)
      dev.off()
    }
    
    # Rank gap
    # Make map with pop cutoff
    pdf(paste('output/maps/outcomes/map_rank_',incs,'_',yr,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
    gplt = ggplot(mdta, aes(fill = rgap_cats_popcut)) +
      geom_sf(size = .1, color = 'black') +
      scale_fill_brewer(name = 'Percentile rank of median Black resident in white family income distribution\n',
                        palette = 'YlGnBu', 
                        labels = c(levels(mdta$rgap_cats_popcut),paste('Black population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5,nrow = 2)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) # This is working if I don't set the CRS at the top
    dev.off()
    
  }
  
  # Make 2-panel map for paper
  if(rdef == 'expansive' & incs == 'all_faminc_sqrt'){
    
    dtaout = filter(regint, year %in% c(1969,2017),  racedef == rdef) %>% select(msa,year,bw_rat_med_cats_popcut)
    dtaout$year = factor(dtaout$year, levels = c(1969, 2017), labels = c('A. 1969', 'B. 2017'))
      mdta = left_join(msashape, dtaout, by = 'msa')
    
    pdf(paste('output/maps/map_paper_panel.pdf',sep = ''),height = 10)
    gplt = ggplot(mdta, aes(fill = bw_rat_med_cats_popcut)) +
      geom_sf(size = .1, color = 'black') + 
      facet_wrap(vars(year),ncol = 1) + 
      scale_fill_brewer(name = 'Black-white median family income ratio\n',
                        palette = 'YlGnBu', 
                        labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) 
    dev.off()
    
    # Black and white for print version
    pdf(paste('output/maps/map_paper_panel_bw.pdf',sep = ''),height = 10)
    gplt = ggplot(mdta, aes(fill = bw_rat_med_cats_popcut)) +
      geom_sf(size = .1, color = 'black') + 
      facet_wrap(vars(year),ncol = 1) + 
      scale_fill_grey(name = 'Black-white median family income ratio\n',na.value = 'white',start = .9,end = .1,
                        # palette = 'Greys', 
                      # breaks = c(.)
                        labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      # scale_fill_brewer(name = 'Black-white median family income ratio\n',
      #                   palette = 'Greys', 
      #                   labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000')),
      #                   guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) 
    dev.off()
    

  }
}





# Now do change
mch = regint %>% filter(year %in% c(1969,2017)) %>% select(year,msa,racedef,pop_b,med_rat,rgap) %>%
  pivot_wider(id_cols = c(msa,racedef),names_from = year,values_from = c(pop_b,med_rat,rgap))

if(incs %in% c('gq_hh_spmuinc_disp_sqrt','gq_hh_spmuinc_sqrt')) next()

mch$change = mch$med_rat_2017 - mch$med_rat_1969
mch$changer = mch$rgap_2017 - mch$rgap_1969

# Place next to abluquerque is missing - replace with ABQ which is MSA 10740
filter(mch, msa == '10740')
mch[mch$msa == '24380', 'change'] = 0.187
mch[mch$msa == '24380', 'changer'] = 0.123
mch[mch$msa == '24380', 'pop_b_1969'] = 6355

# popcut = 10000
chcuts = c(-100000, -.15,-.05,.05,.15,10000)
chlabs = c('Decline over 15 pp','Decline 5-15 pp','Minimal change','Increase 5-15 pp','Increase over 15 pp')
mch$changecat = mch$change %>% cut(chcuts,labels = chlabs ,include.lowest = T)
mch$changecat_popcut = mch$changecat
mch[mch$pop_b_1969 < popcut, 'changecat_popcut'] = NA

chrcuts = c(-100000,-.10, -.05,.05,.1,10000)
chrlabs = c('Decline over 10 percentiles','Decline 5-10 percentiles','Minimal change','Increase 5-10 percentiles','Increase over 10 percentiles')
mch$changercat = mch$changer %>% cut(chrcuts,labels = chrlabs ,include.lowest = T)
mch$changercat_popcut = mch$changercat
mch[mch$pop_b_1969 < popcut, 'changercat_popcut'] = NA

# Make map here
for(rdef in c('expansive')){

mdta = left_join(msashape, mch %>% filter(racedef == rdef), by = 'msa')
  
# Make map with 10k cutoff
pdf(paste('output/maps/outcomes/map_change_',incs,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
gplt = ggplot(mdta, aes(fill = changecat_popcut)) +
  geom_sf(size = .1, color = 'black') + coord_sf(lims_method = 'geometry_bbox',default_crs = NULL) +
  # facet_wrap(.~race_long,ncol = 2) +
  scale_fill_brewer(palette = 'RdYlBu',
                    name = paste('1969-2017 change in Black-white median family income ratio'),
                    labels = c(levels(mdta$changecat_popcut),paste('Black population < 10,000 in 1969')),
                    guide = guide_legend(title.position = 'top',title.hjust = .5)) +
  # scale_fill_brewer(name = '',palette = 'YlGnBu', labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000'))) +
  theme(legend.position = 'bottom')+
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = .5)
  )
print(gplt)
dev.off()

if(rdef == 'expansive' & incs == 'all_faminc_sqrt'){
  pdf(paste('output/maps/map_fig2_ch_',incs,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
 print(gplt)
 dev.off()
}


# Make map with 10k cutoff
pdf(paste('output/maps/outcomes/map_paper_changergap_',incs,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
gplt = ggplot(mdta, aes(fill = changercat_popcut)) +
  geom_sf(size = .1, color = 'black') + coord_sf(lims_method = 'geometry_bbox',default_crs = NULL) +
  # facet_wrap(.~race_long,ncol = 2) +
  scale_fill_brewer(palette = 'RdYlBu',
                    name = paste('1969-2017 change in Black-white median rank gap'),
                    labels = c(levels(mdta$changercat_popcut),paste('Black population < 10,000 in 1969')),
                    guide = guide_legend(title.position = 'top',title.hjust = .5)) +
  # scale_fill_brewer(name = '',palette = 'YlGnBu', labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000'))) +
  theme(legend.position = 'bottom')+
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = .5)
  )
print(gplt)
dev.off()





}

  }
