diff --git a/1_fetch.R b/1_fetch.R
index 810731704e67f75d30dad0981b75bdeff9930948..2e43036cf171b6487a61023d191e52415b2e854f 100644
--- a/1_fetch.R
+++ b/1_fetch.R
@@ -46,8 +46,10 @@ p1_targets <- list(
   tar_target(p1_pop_density_raster_zip,
              '1_fetch/in/gpw-v4-population-count-rev11_2020_30_sec_tif.zip'),
   tar_target(p1_pop_density_raster_tif,
-             unzip(p1_pop_density_raster_zip, 'gpw_v4_population_count_rev11_2020_30_sec.tif',
-                   exdir = p1_out_data)),
+             {unzip(p1_pop_density_raster_zip, 'gpw_v4_population_count_rev11_2020_30_sec.tif',
+                   exdir = p1_out_data)
+               file_name <- paste0(p1_out_data, "gpw_v4_population_count_rev11_2020_30_sec.tif")
+               return(file_name)}),
   # raster data for impervious surfaces
   tar_target(p1_imp_surf_zip,
              sb_initialize_and_download(sb_id = "664e0da6d34e702fe8744579",
@@ -55,7 +57,9 @@ p1_targets <- list(
                               destinations = paste0(p1_out_data, "Annual_NLCD_FctImp_2022_CU_C1V0.zip")),
              format = 'file'),
   tar_target(p1_imp_surf_tif,
-             unzip(p1_imp_surf_zip, 'Annual_NLCD_FctImp_2022_CU_C1V0.tif',
-                   exdir = p1_out_data))
+             {unzip(p1_imp_surf_zip, 'Annual_NLCD_FctImp_2022_CU_C1V0.tif',
+                   exdir = p1_out_data)
+              file_name <- paste0(p1_out_data, "Annual_NLCD_FctImp_2022_CU_C1V0.tif")
+              return(file_name)})
   )
   
\ No newline at end of file
diff --git a/2_process.R b/2_process.R
index bc448f702b3d4cbd18d96d68d31b6e9028a6c1ce..af6c3b31d5e1f57debbdc9243244af461fe026e9 100644
--- a/2_process.R
+++ b/2_process.R
@@ -251,6 +251,14 @@ tar_target(p2_conus_sf_proj,
              st_transform(p1_proj)),
 tar_target(p2_conus_inner,
            rmapshaper::ms_innerlines(p2_conus_sf_proj)),
+tar_target(p2_conus_counties_sf,
+           counties_sf <- tigris::counties(cb = TRUE, state = c('Washington', 'Oregon', 'California', 'Idaho', 'Nevada',
+                                                                'Utah', 'Arizona', 'Montana', 'Wyoming', 'Colorado',
+                                                                'New Mexico', 'North Dakota', 'South Dakota', 'Nebraska', 'Kansas',
+                                                                'Oklahoma', 'Texas', 'Minnesota', 'Iowa', 'Missouri',
+                                                                'Arkansas', 'Louisiana')) |> 
+             rmapshaper::ms_simplify(keep = 0.2) |>  
+             st_transform('EPSG:5070')),
 tar_target(p2_pop_density_processed,
            process_pop_dens_raster(in_raster = p1_pop_density_raster_tif, #proj = p1_proj, 
                                    conus = p2_conus_sf, conus_proj = p2_conus_sf_proj,
diff --git a/3_visualize.R b/3_visualize.R
index c8c14b1d9c7af23d7726db1ee62be5a9916d868f..77b879518f703d1af5fdceb9a322e99e447d6808 100644
--- a/3_visualize.R
+++ b/3_visualize.R
@@ -121,7 +121,7 @@ tar_target(
   plot_raster(in_raster = p2_pop_density_processed, conus_sf = p2_conus_sf_proj, viz_config_df = p0_viz_config_df,
               conus_inner = p2_conus_inner, viz_config_pal = p0_viz_config_pal$living_conditions,
               width = p0_viz_config_df$width_desktop, height = p0_viz_config_df$height_desktop, 
-              font_size  = p0_viz_config_df$font_size_desktop, barwidth = 20, barheight = 1,
+              font_size  = p0_viz_config_df$font_size_desktop, barwidth = 20, barheight = 1, counties_sf = p2_conus_counties_sf,
               outfile_path = "3_visualize/out/pop_density_rast_2020_en.png", low_ramp_col = "#E8F4F4"),
   format = "file"),
 # impervious surfaces raster plot
@@ -131,8 +131,8 @@ tar_target(
                        conus_col = "gray90", conus_inner = p2_conus_inner, viz_config_pal = p0_viz_config_pal$living_conditions,
                        width = p0_viz_config_df$width_desktop, height = p0_viz_config_df$height_desktop, 
                        font_size  = p0_viz_config_df$font_size_desktop, barwidth = 20, barheight = 1,
-                       outfile_path = "3_visualize/out/imp_surface_rast_2022_en.png", 
-                       low_ramp_col = "white" )
+                       outfile_path = "3_visualize/out/imp_surface_rast_2022_en.png", counties_sf = p2_conus_counties_sf,
+                       low_ramp_col = "white")
 ),
 # Spanish version's of maps -----------------------------------------------
 tar_target(
diff --git a/3_visualize/src/plot_utils.R b/3_visualize/src/plot_utils.R
index 9ac4f016241a9d30883b33c0bdf383887ce6d69a..626d1aab1d56642b979e04e663e198011c6ff98a 100644
--- a/3_visualize/src/plot_utils.R
+++ b/3_visualize/src/plot_utils.R
@@ -115,7 +115,7 @@ ggsave(outfile_path, final_map, width = width, height = height, dpi = viz_config
 
 plot_raster <- function(in_raster, conus_sf, conus_inner, outfile_path, 
                         viz_config_df, viz_config_pal, width, height,
-                        font_size, barwidth, barheight,
+                        font_size, barwidth, barheight, counties_sf,
                         low_ramp_col){
   # in_raster = p2_pop_density_processed, conus_sf = p2_conus_sf_proj, conus_inner = p2_conus_inner, viz_config_pal = p0_viz_config_pal$living_conditions
   # width = p0_viz_config_df$width_desktop, height = p0_viz_config_df$height_desktop, font_size  = p0_viz_config_df$font_size_desktop, barwidth = 20, barheight = 1,
@@ -129,27 +129,21 @@ plot_raster <- function(in_raster, conus_sf, conus_inner, outfile_path,
   showtext_auto(enable = TRUE)
   
   # doesn't work with cowplot::get_plot_component()
-  legend_label <- expression("Population per 1 km"^2) #bquote('Population per 1'~km^2)
+  legend_label <- expression(bold("Population per 1 km"^2)) #bquote('Population per 1'~km^2)
   
   (pop_map <- ggplot() +
       #geom_spatraster(data = pop_usa_proj)+ 
       geom_sf(data = conus_sf, fill = low_ramp_col, color = low_ramp_col, linewidth = 0.5) +
       geom_spatraster(data = raster_data)+ #x = x, y = y, 
       geom_sf(data = conus_inner, fill = NA, color = viz_config_df$counties_outline_col, linewidth = 0.5) +
+      geom_sf(data = counties_sf, fill = NA, color = viz_config_df$counties_outline_col, linewidth = 0.1) +
       scale_fill_gradientn(
         colors = colorRampPalette(c(low_ramp_col, viz_config_pal))(100), 
         name = legend_label,
-        #limits = lim_vals,  # c(0, 100),
-        #breaks = break_vals, # c(0, 25, 50, 75, 100),
-        #labels = if (dollar_leg) scales::label_dollar() else function(x) paste0(x, "%"),
         labels = c(1, 10, 100, '1k', '10k', '100k'),
         na.value=NA
       )+
-      #scale_fill_viridis_c(na.value = NA, option = 'mako', breaks = c(0, 1, 2, 3, 4, 5), direction = -1,
-                           # leave a little room for the NAs and Inf (which are 0s)
-                           # to be darker than the -1 values
-       #                    labels = c(1, 10, 100, '1k', '10k', '100k'), begin = 0.05) +
-      labs(fill = legend_label) + #"Population per 1 km^2"
+      labs(fill = legend_label) + 
       theme(plot.background = element_rect(fill = "white", color = "white"),
             panel.background = element_rect(fill = "white", color = "white"),
             axis.title = element_blank(), 
@@ -168,13 +162,6 @@ plot_raster <- function(in_raster, conus_sf, conus_inner, outfile_path,
         barwidth = barwidth, 
         barheight = barheight 
       ))
-    # guides(fill = guide_colorbar(
-    #    direction = "horizontal",
-    #    title = "Population (per km2)",
-    #    barwidth = 15, 
-    #    title.position = "top",
-    #    title.theme = element_text(color = "white"),
-    #    label.theme = element_text(color = "white")))
   )
   
   background_color = "white"
@@ -216,7 +203,7 @@ plot_raster <- function(in_raster, conus_sf, conus_inner, outfile_path,
 
 plot_imp_surf_raster <- function(in_raster, conus_sf, conus_inner, conus_col, outfile_path, 
                         viz_config_df, viz_config_pal, width, height,
-                        font_size, barwidth, barheight,
+                        font_size, barwidth, barheight, counties_sf,
                         low_ramp_col){
   # in_raster = p2_imp_surf_processed, conus_sf = p2_conus_sf_proj, conus_col = "gray90", viz_config_pal = p0_viz_config_pal$living_conditions
   # width = p0_viz_config_df$width_desktop, height = p0_viz_config_df$height_desktop, font_size  = p0_viz_config_df$font_size_desktop, barwidth = 20, barheight = 1,
@@ -231,8 +218,9 @@ plot_imp_surf_raster <- function(in_raster, conus_sf, conus_inner, conus_col, ou
   
   (pop_map <- ggplot() +
       geom_sf(data = conus_sf, fill = NA, color = conus_col, linewidth = 0.9) +
-      geom_spatraster(data = raster_data)+ #x = x, y = y, 
+      geom_spatraster(data = raster_data)+ 
       geom_sf(data = conus_inner, fill = NA, color = conus_col, linewidth = 0.5) +
+      geom_sf(data = counties_sf, fill = NA, color = conus_col, linewidth = 0.1) +
       scale_fill_gradientn(
         colors = colorRampPalette(c(low_ramp_col, viz_config_pal))(100), 
         name = "Percent impervious surface, 2022",