Merge branch 'main' of...

Merge branch 'main' of https://code.usgs.gov/wma/vizlab/vulnerability-indicators into add_python_notebook

.gitignore

@@ -4,7 +4,7 @@
 .RData
 .Ruserdata
 .RProfile
-.Renvironment
+.Renviron
 .DS_Store
 1_fetch/out/*
 1_fetch/in/*

00_config.R

+source("0_config/src/sb_cache.R")
+
+initialize_and_cache_sb(
+  sb_username = rstudioapi::showPrompt("Username", "SB Username"),
+  renviron_file = ".Renviron",
+  override_login = FALSE
+)

0_config/src/sb_cache.R

+" Initialize ScienceBase login and cache token in .Renviron
+#"
+#' @param sb_username chr; ScienceBase username. The default prompts the user to
+#'   enter the username in the console.
+#' @param renviron_file chr; path to .Renviron file
+#' @param override_login lgl; should ScienceBase be re-initialized if already
+#'   logged in?
+#'
+#' @return lgl; result of `sbtools::is_logged_in()`
+#'
+initialize_and_cache_sb <- function(sb_username = NULL,
+                                    renviron_file = ".Renviron",
+                                    override_login = FALSE) {
+  # Give error if already logged in and `override_login = FALSE`
+  tryCatch(
+    sbtools:::token_refresh(),
+    warning = function(x) {},
+    error = function(x) FALSE
+  )
+  
+  if (sbtools::is_logged_in()) {
+    if (override_login) {
+      # If logged in and overriding login: warn and re-initialize SB
+      cli::cli_warn(c(
+        "!" = "You are already logged into ScienceBase and re-initialization
+        anyways.",
+        "i" = "If you would not like to re-initialize, abort and re-run with
+        {.arg override_login = FALSE}"
+      ))
+      sbtools::initialize_sciencebase_session(username = sb_username)
+    }
+  } else {
+    # If not logged into SB, initialize SB
+    sbtools::initialize_sciencebase_session(username = sb_username)
+  }
+  
+  # Grab token and update it in .Renviron
+  update_renviron(
+    sb_username = sb_username,
+    sb_token = jsonlite::toJSON(
+      sbtools::current_session()[c("access_token", "refresh_token")]
+    ),
+    renviron_file
+  )
+  
+  # Output if currently logged into SB
+  if (sbtools::is_logged_in()) {
+    return(TRUE)
+  } else {
+    cli::cli_abort(c(
+      "!" = "You are not logged into ScienceBase after attempting it.",
+      "i" = "Try running {.fun initialize_and_cache_sb} again."
+    ))
+  }
+}
+
+#' Create/update .Renviron with ScienceBase token
+#'
+#' @param sb_token chr; token used to initialize a ScienceBase session
+#' @param renviron_file chr; path to .Renviron file
+#'
+#' @return chr; path to .Renviron file
+#'
+update_renviron <- function(sb_username, sb_token, renviron_file) {
+  if (file.exists(renviron_file)) {
+    # If there is an existing .Renviron file...
+    # Read in existing .Renviron file and check if there is an "sb_token" value
+    existing <- readLines(renviron_file)
+    sb_token_idx <- which(startsWith(existing, "sb_token="))
+    sb_username_idx <- which(startsWith(existing, "sb_username="))
+    
+    # If there is already a value for sb_token update it; if not, create one.
+    if (length(sb_token_idx) > 0) {
+      existing[sb_token_idx] <- paste0("sb_token=", sb_token)
+      writeLines(existing, con = renviron_file)
+    } else {
+      write(paste0("sb_token=", sb_token), file = renviron_file, append = TRUE)
+    }
+    
+    # If there is already a value for sb_username update it; if not, create one
+    if (length(sb_username_idx) > 0) {
+      existing[sb_username_idx] <- paste0("sb_username=", sb_username)
+      writeLines(existing, con = renviron_file)
+    } else {
+      write(
+        paste0("sb_username=", sb_username),
+        file = renviron_file, append = TRUE
+      )
+    }
+  } else {
+    # If there isn't an .Renviron file, create one with an "sb_token" value
+    write(
+      paste0("sb_token=", sb_token, "\n", "sb_username=", sb_username),
+      file = renviron_file
+    )
+  }
+  
+  return(renviron_file)
+}

1_fetch.R

@@ -2,19 +2,21 @@ source("1_fetch/src/fetch_utils.R")
 
 p1_targets <- list(
   ##### Vulnerability indicators data #####
+  tar_target(p1_out_data,
+             "1_fetch/out/"),
   tar_target(p1_sb_id,
              '63f79d49d34e4f7eda456572'),
   ##### Vulnerability indicators metadata #####
   tar_target(p1_vul_ind_xml,
-             download_from_sb(sb_id = p1_sb_id,
-                              filename = 'Uncertainty_Summary.xml',
-                              dest_dir = '1_fetch/out'),
+             sb_initialize_and_download(sb_id = p1_sb_id,
+                              names = 'Uncertainty_Summary.xml',
+                              destinations = paste0(p1_out_data, "Uncertainty_Summary.xml")),
              format = 'file'),
   ##### Uncertainty statistics for the indicators #####
   tar_target(p1_unc_stats_csv,
-             download_from_sb(sb_id = p1_sb_id,
-                              filename = 'Uncertainty_Summary.csv',
-                              dest_dir = '1_fetch/out'),
+             sb_initialize_and_download(sb_id = p1_sb_id,
+                              names = 'Uncertainty_Summary.csv',
+                              destinations = paste0(p1_out_data, "Uncertainty_Summary.csv")),
              format = 'file'),
   tar_target(p1_unc_stats,
              readr::read_csv(p1_unc_stats_csv) |> 
@@ -39,5 +41,25 @@ p1_targets <- list(
                )) |> 
                filter(group %in% c('CONUS')) |> 
                rmapshaper::ms_simplify(keep = 0.2) |> 
-               filter(NAME %in% p1_census_states))
-  )
\ No newline at end of file
+               filter(NAME %in% p1_census_states)),
+  # raster data for population density
+  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)
+               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",
+                              names = 'Annual_NLCD_FctImp_2022_CU_C1V0.zip',
+                              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)
+              file_name <- paste0(p1_out_data, "Annual_NLCD_FctImp_2022_CU_C1V0.tif")
+              return(file_name)})
+  )
+  
\ No newline at end of file

1_fetch/src/fetch_utils.R

-#' @title Download file from ScienceBase
-#' @description Download a file from ScienceBase and save it to a specified
-#' directory
-#' @param sb_id id for ScienceBase data release
-#' @param filename filename of item in data release
-#' @param dest_dir directory to which to download the file
-#' @return the filepath of the saved file
-download_from_sb <- function(sb_id, filename, dest_dir) {
-  filepath <- file.path(dest_dir, filename)
-  
-  sbtools::item_file_download(sb_id = sb_id,
-                              names = filename,
-                              destinations = filepath)
-  
-  return(filepath)
-}
\ No newline at end of file
+#' Initialize ScienceBase session and download files
+#'
+#' @param sb_id chr; ScienceBase ID
+#' @param names chr; names of files to download from ScienceBase
+#' @param destinations  chr; write path location for downloaded files
+#' @param renviron_file chr; path to .Renviron file where credentials are cached
+#' @param ... additional arguments passed to `sbtools::item_file_download()`
+#'
+#' @return chr; path to downloaded files
+#'
+sb_initialize_and_download <- function(sb_id, names, destinations,
+                                       renviron_file = ".Renviron", ...) {
+  # Initialize ScienceBase session
+  sb_login_cached(renviron_file = renviron_file)
+  
+  # Download SB files
+  sbtools::item_file_download(
+    sb_id = sb_id,
+    names = names,
+    destinations = destinations,
+    ...
+  )
+  
+  return(destinations)
+}
+
+#' Login to ScienceBase using cached credentials
+#'
+#' @param renviron_file chr; path to .Renviron file
+#'
+#' @return `TRUE` if logged in. Error if not.
+#'
+sb_login_cached <- function(renviron_file) {
+  # If logged in, return TRUE and skip the rest
+  if (sbtools::is_logged_in()) {
+    return(TRUE)
+  }
+  
+  # Try a token refresh
+  tryCatch(
+    sbtools:::token_refresh(),
+    warning = function(x) {},
+    error = function(x) FALSE
+  )
+  
+  if (sbtools::is_logged_in()) {
+    return(TRUE)
+  }
+  
+  # If .Renviron file does not exist, re-initialize
+  if (!file.exists(renviron_file)) {
+    cli::cli_abort(c(
+      "Could not find the specified file: {.file {renviron_file}}.",
+      "i" = "Follow the instructions in {.file README.md} to initalize and cache
+      ScienceBase login credentials."
+    ))
+  }
+  
+  # Read .Renviron file
+  existing <- readLines(renviron_file)
+  sb_token_idx <- which(startsWith(existing, "sb_token="))
+  sb_username_idx <- which(startsWith(existing, "sb_username="))
+  
+  # If SB credentials not found, throw error
+  if (any(length(sb_token_idx) == 0, length(sb_username_idx) == 0)) {
+    cli::cli_abort(c(
+      "Could not find the username or token in the specified file:
+      {.file {renviron_file}}.",
+      "i" = "Follow the instructions in {.file README.md} to re-initalize and
+      cache ScienceBase login credentials."
+    ))
+  }
+  
+  # Get ScienceBase credentials
+  sb_token <- stringr::str_remove(existing[sb_token_idx], "^sb_token=")
+  sb_username <- stringr::str_remove(existing[sb_username_idx], "^sb_username=")
+  
+  # Initialize ScienceBase session with cached credentials
+  sbtools::initialize_sciencebase_session(
+    username = sb_username,
+    token_text = sb_token
+  )
+  
+  if (sbtools::is_logged_in()) {
+    return(TRUE)
+  } else {
+    cli::cli_abort(c(
+      "Could not login to ScienceBase using cached credentials.",
+      "i" = "Follow the instructions in {.file README.md} to re-initalize and
+      cache ScienceBase login credentials."
+    ))
+  }
+}

2_process.R

@@ -70,11 +70,12 @@ p2_targets <- list(
                # distinct(determinant, .keep_all = TRUE) |> 
                readr::write_csv('public/determinant_uncertainty.csv')
              ),
-tar_target(p2_unc_determinant_json,
-           read_csv(p2_unc_agg_summary_csv) |>
-             toJSON(pretty = TRUE) |>
-             write("public/determinant_uncertainty.json")
-           ),
+# commented out for now so we don't overwrite spanish names
+#tar_target(p2_unc_determinant_json,
+#           read_csv(p2_unc_agg_summary_csv) |>
+#             toJSON(pretty = TRUE) |>
+#             write("public/determinant_uncertainty.json")
+#           ),
 tar_target(`p2_unc_agg_summary_ind_csv`,
            p2_unc_agg_ind_summary |>
              left_join(p2_top_trend_ind_stats) |>
@@ -241,5 +242,29 @@ tar_target(p2_census_acs5sub_disability_data,
                            survey_var = "acs5", 
                            percent_rename = FALSE),
            pattern = map(p2_census_acs5_disability_layers),
-           iteration = "list")
+           iteration = "list"),
+# prep for raster data processing and plotting
+tar_target(p2_conus_sf,
+           fetch_conus_sf(states = p1_census_states)),
+tar_target(p2_conus_sf_proj,
+           p2_conus_sf |>  
+             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 = p1_census_states) |> 
+             rmapshaper::ms_simplify(keep = 0.2) |>  
+             st_transform('EPSG:5070')),
+# process population density raster data
+tar_target(p2_pop_density_processed,
+           process_pop_dens(in_raster = p1_pop_density_raster_tif, 
+                            conus_sf = p2_conus_sf, conus_proj = p2_conus_sf_proj,
+                            outfile_path = "2_process/out/pop_density.tif"),
+           format = "file"),
+# process impervious surfaces raster data
+tar_target(p2_imp_surf_processed,
+           process_imp_surf(in_raster = p1_imp_surf_tif, conus_proj = p2_conus_sf_proj,
+                            outfile_path = "2_process/out/imp_surfaces.tif"),
+           format = "file")
+
 )
\ No newline at end of file

2_process/src/data_utils.R

@@ -47,4 +47,73 @@ process_perc <- function(tot_pop, tot_var){
   return(joined_perc_df)}
 prep_tree_data <- function(data) {
   
+}
+
+## data processing for population density raster map
+#' @title Fetch conus for raster processing 
+#' @param states, chr string with long name of states of interest to filter by
+fetch_conus_sf <- function(states){
+  
+  conus_sf <- tigris::states(cb = TRUE) |> 
+    dplyr::filter(NAME %in% states
+    ) |> 
+    rmapshaper::ms_simplify(keep = 0.2) 
+}
+
+#' @title Process population density raster data for plotting
+#' @param in_raster character string - .tif file path of raster data for specified variable of interest
+#' @param conus_sf, sf of conus states outline
+#' @param conus_proj, projected sf of conus states outline
+#' @param outfile_path, outfile path for processed tifs 
+process_pop_dens <- function(in_raster, conus_sf, conus_proj, outfile_path){
+  
+  raw_data <- rast(in_raster)
+
+  # crop population data to area of interest 
+  terra::window(raw_data) <- terra::ext(conus_sf) 
+  
+  # project population data
+  conus_sf_rast <- rast(conus_proj, resolution = c(1000, 1000)) 
+  pop_usa_proj <- project(raw_data, conus_sf_rast)
+  # match boundaries of population data to conus data
+  pop_usa_mask <- terra::mask(pop_usa_proj, terra::vect(conus_proj))
+  
+  # change to tibble using tidy terra for processing into log scale 
+  usa_dat <- as_tibble(pop_usa_mask, xy = TRUE)
+  usa_dat$pop <- ifelse(usa_dat$gpw_v4_population_count_rev11_2020_30_sec < 1 
+                        & usa_dat$gpw_v4_population_count_rev11_2020_30_sec > 0, 0.1, 
+                        usa_dat$gpw_v4_population_count_rev11_2020_30_sec)
+  usa_dat$pop_log10 <- log10(usa_dat$pop)
+  usa_dat <- usa_dat |> 
+    dplyr::select(-gpw_v4_population_count_rev11_2020_30_sec, -pop)
+  # convert back into raster
+  usa_dat_rast <- as_spatraster(usa_dat) 
+  
+  writeRaster(usa_dat_rast, filename = outfile_path, overwrite=TRUE)
+  
+  return(outfile_path)
+}
+
+
+#' @title Process impervous surfaces raster data for plotting
+#' @param in_raster character string - .tif file path of raster data for specified variable of interest
+#' @param conus_proj, projected sf of conus states outline
+#' @param outfile_path, outfile path for processed tifs 
+process_imp_surf <- function(in_raster, conus_proj, outfile_path){
+  
+  imp_surf_raw <- rast(in_raster)
+  
+  # crop population data to area of interest
+  terra::window(imp_surf_raw) <- terra::ext(conus_proj)
+  
+  # project population data
+  conus_sf_rast <- rast(conus_proj, resolution = c(1000, 1000)) 
+  imp_surf_proj <- project(imp_surf_raw, conus_sf_rast)
+  
+  # match boundaries of population data to conus data
+  imp_surf_mask <- terra::mask(imp_surf_proj, terra::vect(conus_proj))
+  
+  writeRaster(imp_surf_mask, filename = outfile_path, overwrite=TRUE)
+  
+  return(outfile_path)
 }
\ No newline at end of file

3_visualize.R

@@ -93,139 +93,183 @@ p3_targets <- list(
     format = "file"
   ),
   tar_target(
-  p3_perc_disable_png_en,
-  plot_census_map(
-    census_data = p2_census_acs5sub_disability_data[[1]],
-    percent_leg = TRUE,
-    dollar_leg = FALSE,
-    lim_vals = c(0, 50),
-    break_vals = c(0, 10, 20, 30, 40, 50),
-    var = 'estimate',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/perc_disabled_census_2022_en.png",
-    leg_title = "Percent disabled, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$demographic_characteristics,
-    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,
-    low_ramp_col = "#eef0ff"
+    p3_perc_disable_png_en,
+    plot_census_map(
+      census_data = p2_census_acs5sub_disability_data[[1]],
+      percent_leg = TRUE,
+      dollar_leg = FALSE,
+      lim_vals = c(0, 50),
+      break_vals = c(0, 10, 20, 30, 40, 50),
+      var = 'estimate',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/perc_disabled_census_2022_en.png",
+      leg_title = "Percent disabled, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$demographic_characteristics,
+      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,
+      low_ramp_col = "#eef0ff"
     ),
-  format = "file"
-),
-# Spanish version's of maps -----------------------------------------------
-tar_target(
-  p3_med_income_png_es,
-  plot_census_map(
-    census_data = p2_perc_census_acs5_layers_sf[[2]],
-    lim_vals = c(0, 155000),
-    percent_leg = FALSE,
-    dollar_leg = TRUE,
-    var = 'estimate',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/med_income_census_2022_es.png",
-    leg_title = "Media de ingresos por hogar, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$socioeconomic_status,
-    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,
-    low_ramp_col = "#fef1f1"
+    format = "file"
   ),
-  format = "file"
-),
-tar_target(
-  p3_perc_latino_png_es,
-  plot_census_map(
-    census_data = p2_perc_census_acs5_layers_sf[[4]],
-    percent_leg = TRUE,
-    dollar_leg = FALSE,
-    lim_vals = c(0, 100),
-    break_vals = c(0, 25, 50, 75, 100),
-    var = 'percent',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/perc_hispanic_census_2022_es.png",
-    leg_title = "Porcentaje de Hispanos, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$demographic_characteristics,
-    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,
-    low_ramp_col = "#eef0ff"
+  # population density raster plot
+  tar_target(
+    p3_pop_density_png,
+    plot_raster(in_raster = p2_pop_density_processed, 
+                conus_sf = p2_conus_sf_proj,
+                conus_inner = p2_conus_inner,
+                counties_sf = p2_conus_counties_sf,
+                viz_config_df = p0_viz_config_df,
+                viz_config_pal = "#38837a",
+                rast_type = "pop_density",
+                low_ramp_col = "#E8F4F4",
+                base_conus_fill = "#E8F4F4", 
+                base_conus_color = "#E8F4F4", 
+                border_color = p0_viz_config_df$counties_outline_col, 
+                leg_title = expression(paste(bold("Population per 1 "), bold(km^2), bold(", 2020"))),
+                font_size  = p0_viz_config_df$font_size_desktop,
+                barheight = 1,
+                barwidth = 20,
+                outfile_path = "3_visualize/out/pop_density_rast_2020_en.png",
+                width = p0_viz_config_df$width_desktop,
+                height = p0_viz_config_df$height_desktop),
+    format = "file"),
+  # impervious surfaces raster plot
+  tar_target(
+    p3_imp_surf_png,
+    plot_raster(in_raster = p2_imp_surf_processed,
+                conus_sf = p2_conus_sf_proj,
+                conus_inner = p2_conus_inner,
+                counties_sf = p2_conus_counties_sf,
+                viz_config_df = p0_viz_config_df,
+                viz_config_pal = "#38837a",
+                rast_type = "imp_surfaces",
+                low_ramp_col = "white",
+                base_conus_fill = NA, 
+                base_conus_color = "#E8F4F4", 
+                border_color = "gray90", 
+                leg_title = "Percent impervious surface, 2022",
+                font_size  = p0_viz_config_df$font_size_desktop,
+                barheight = 1,
+                barwidth = 20,
+                outfile_path = "3_visualize/out/imp_surface_rast_2022_en.png",
+                width = p0_viz_config_df$width_desktop,
+                height = p0_viz_config_df$height_desktop),
+    format = "file"),
+  # Spanish version's of maps -----------------------------------------------
+  tar_target(
+    p3_med_income_png_es,
+    plot_census_map(
+      census_data = p2_perc_census_acs5_layers_sf[[2]],
+      lim_vals = c(0, 155000),
+      percent_leg = FALSE,
+      dollar_leg = TRUE,
+      var = 'estimate',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/med_income_census_2022_es.png",
+      leg_title = "Media de ingresos por hogar, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$socioeconomic_status,
+      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,
+      low_ramp_col = "#fef1f1"
+    ),
+    format = "file"
   ),
-  format = "file"
-),
-tar_target(
-  p3_avg_household_size_png_es,
-  plot_census_map(
-    census_data = p2_census_acs5sub_household_data[[1]],
-    percent_leg = FALSE,
-    dollar_leg = FALSE,
-    lim_vals = c(1, 5),
-    var = 'estimate',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/avg_household_size_2022_es.png",
-    leg_title = "Tamaño promedio de los hogares, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$demographic_characteristics,
-    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,
-    low_ramp_col = "#eef0ff"
+  tar_target(
+    p3_perc_latino_png_es,
+    plot_census_map(
+      census_data = p2_perc_census_acs5_layers_sf[[4]],
+      percent_leg = TRUE,
+      dollar_leg = FALSE,
+      lim_vals = c(0, 100),
+      break_vals = c(0, 25, 50, 75, 100),
+      var = 'percent',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/perc_hispanic_census_2022_es.png",
+      leg_title = "Porcentaje de Hispanos, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$demographic_characteristics,
+      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,
+      low_ramp_col = "#eef0ff"
+    ),
+    format = "file"
   ),
-  format = "file"
-),
-tar_target(
-  p3_median_rent_png_es,
-  plot_census_map(
-    census_data = p2_census_acs5sub_household_data[[2]],
-    percent_leg = FALSE,
-    dollar_leg = TRUE, 
-    lim_vals = c(0, 3000),
-    var = 'estimate',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/median_rent_2022_es.png",
-    leg_title = "Alquiler bruto medio, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$land_tenure,
-    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,
-    low_ramp_col = "#ebfaf8"
+  tar_target(
+    p3_avg_household_size_png_es,
+    plot_census_map(
+      census_data = p2_census_acs5sub_household_data[[1]],
+      percent_leg = FALSE,
+      dollar_leg = FALSE,
+      lim_vals = c(1, 5),
+      var = 'estimate',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/avg_household_size_2022_es.png",
+      leg_title = "Tamaño promedio de los hogares, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$demographic_characteristics,
+      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,
+      low_ramp_col = "#eef0ff"
+    ),
+    format = "file"
   ),
-  format = "file"
-),
-tar_target(
-  p3_perc_disable_png_es,
-  plot_census_map(
-    census_data = p2_census_acs5sub_disability_data[[1]],
-    percent_leg = TRUE,
-    dollar_leg = FALSE,
-    lim_vals = c(0, 50),
-    break_vals = c(0, 10, 20, 30, 40, 50),
-    var = 'estimate',
-    conus_sf = p1_conus_sf,
-    outfile_path = "3_visualize/out/perc_disabled_census_2022_es.png",
-    leg_title = "Porcentaje de discapacitados, 2022",
-    viz_config_df = p0_viz_config_df,
-    viz_config_pal = p0_viz_config_pal$demographic_characteristics,
-    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,
-    low_ramp_col = "#eef0ff"
+  tar_target(
+    p3_median_rent_png_es,
+    plot_census_map(
+      census_data = p2_census_acs5sub_household_data[[2]],
+      percent_leg = FALSE,
+      dollar_leg = TRUE, 
+      lim_vals = c(0, 3000),
+      var = 'estimate',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/median_rent_2022_es.png",
+      leg_title = "Alquiler bruto medio, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$land_tenure,
+      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,
+      low_ramp_col = "#ebfaf8"
+    ),
+    format = "file"
   ),
-  format = "file"
-)
+  tar_target(
+    p3_perc_disable_png_es,
+    plot_census_map(
+      census_data = p2_census_acs5sub_disability_data[[1]],
+      percent_leg = TRUE,
+      dollar_leg = FALSE,
+      lim_vals = c(0, 50),
+      break_vals = c(0, 10, 20, 30, 40, 50),
+      var = 'estimate',
+      conus_sf = p1_conus_sf,
+      outfile_path = "3_visualize/out/perc_disabled_census_2022_es.png",
+      leg_title = "Porcentaje de discapacitados, 2022",
+      viz_config_df = p0_viz_config_df,
+      viz_config_pal = p0_viz_config_pal$demographic_characteristics,
+      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,
+      low_ramp_col = "#eef0ff"
+    ),
+    format = "file"
+  )
 )
\ No newline at end of file

3_visualize/src/plot_utils.R

@@ -17,7 +17,7 @@
 #' @param break_vals, set legend breaks 
 #' @param low_ramp_col, assign color to low end of color ramp for legend 
 plot_census_map <- function(census_data, conus_sf, leg_title, outfile_path, var, 
-                            percent_leg, viz_config_df, viz_config_pal, width, height,
+                            percent_leg, viz_config_df, viz_config_pal, width, height, 
                             font_size, barwidth, barheight, lim_vals, break_vals,
                             dollar_leg, low_ramp_col){
   
@@ -105,4 +105,132 @@ plot_census_map <- function(census_data, conus_sf, leg_title, outfile_path, var,
 ggsave(outfile_path, final_map, width = width, height = height, dpi = viz_config_df$dpi, bg = viz_config_df$bg_col, units = "in")
   
   return(outfile_path)
+}
+
+
+# Map raster data 
+#'
+#' @param in_raster character string - .tif file path of raster data for specified variable of interest 
+#' @param conus_sf, sf of conus states outline
+#' @param conus_inner, sf of conus states outline not including US borders
+#' @param counties_sf, sf of counties outline
+#' @param viz_config_df `data.frame` width, height, counties outline color, conus outline color, background color, font nam, and font size 
+#' @param viz_config_pal `data.frame` assign colors for postively and negatively correlated dimensions for census maps
+#' @param rast_type character string - choice of "pop_density" or "imp_surfaces" to indicate which map is being made
+#' @param low_ramp_col, color for low end of color ramp 
+#' @param base_conus_fill, base fill used for conus background
+#' @param base_conus_color, color used for state outlines
+#' @param border_color, color used for inner state outlines and county outlines
+#' @param leg_title, title to be used for the legend
+#' @param font_size, set font size 
+#' @param barheight, set colorbar bar height
+#' @param barwidth, set colorbar bar width 
+#' @param outfile_path, outfile path for pngs 
+#' @param width, set figure width dimension
+#' @param height, set figure height dimension
+plot_raster <- function(in_raster, conus_sf, conus_inner, counties_sf, 
+                        viz_config_df, viz_config_pal, rast_type, low_ramp_col, 
+                        base_conus_fill, base_conus_color, border_color, 
+                        leg_title, font_size, barheight, barwidth, outfile_path, 
+                        width, height){
+  
+  raster_data <- rast(in_raster)
+  
+  # plot
+  font_legend <- viz_config_df$load_font
+  font_add_google(font_legend)
+  showtext_opts(dpi = 300, regular.wt = 200, bold.wt = 700)
+  showtext_auto(enable = TRUE)
+  
+  rast_map <- ggplot() +
+    geom_sf(data = conus_sf, fill = base_conus_fill, color = base_conus_color, linewidth = 0.9) +
+    geom_spatraster(data = raster_data)+ 
+    geom_sf(data = conus_inner, fill = NA, color = border_color, linewidth = 0.5) +
+    geom_sf(data = counties_sf, fill = NA, color = border_color, linewidth = 0.1) +
+    labs(fill = leg_title) + 
+    theme(plot.background = element_rect(fill = "white", color = "white"),
+          panel.background = element_rect(fill = "white", color = "white"),
+          text = element_text(family = viz_config_df$font_legend, size = font_size),
+          axis.title = element_blank(), 
+          axis.text = element_blank(), 
+          axis.ticks = element_blank(), 
+          panel.grid = element_blank(),
+          legend.background = element_blank(),
+          legend.direction = "horizontal",
+          legend.margin = margin(t = 5, b = 5),
+          legend.position = "bottom",
+          legend.title.align = 0.5) +
+    guides(fill = guide_colorbar(
+      title.position = "top",
+      title.theme = element_text(face = 'bold', family = font_legend, size = font_size),
+      direction = "horizontal",
+      position = "bottom",
+      barwidth = barwidth, 
+      barheight = barheight 
+    ))
+  
+  if(rast_type == "pop_density"){
+    rast_map <- rast_map + 
+      scale_fill_gradientn(
+        colors = colorRampPalette(c(low_ramp_col, viz_config_pal))(100), 
+        name = leg_title,
+        labels = c(1, 10, 100, '1k', '10k', '100k'),
+        na.value=NA
+      )
+   
+  } else if(rast_type == "imp_surfaces") {
+    rast_map <- rast_map +
+      scale_fill_gradientn(
+        colors = colorRampPalette(c(low_ramp_col, viz_config_pal))(100), 
+        name = leg_title,
+        limits = c(0, 100),
+        breaks = c(0, 25, 50, 75, 100),
+        labels = c("0%", "25%", "50%", "75%", "100%"),
+        na.value=NA
+      )
+  }
+  
+  background_color = "white"
+  plot_margin =  0.025
+  
+  # cowplot to get map sizes larger
+  canvas <- grid::rectGrob(
+    x = 0, y = 0,
+    width = 6, height = 6,
+    gp = grid::gpar(fill = background_color, alpha = 1, col = background_color
+    )
+  )
+  rast_map_legend <-  get_plot_component(rast_map, 'guide-box', return_all = TRUE)
+  # compose final plot
+  final_map <- ggdraw(ylim = c(0,1),
+                      xlim = c(0,1)) +
+    # White background
+    draw_grob(canvas,
+              x = 0, y = 1,
+              height = 6, width = 6,
+              hjust = 0, vjust = 1) +
+    # Add main plot
+    draw_plot(rast_map + theme(legend.position="none"),
+              x = -0.01,
+              y = 0.08,
+              height = 0.98,
+              width = (1-plot_margin)*1.03) 
+    # Add legend 
+  
+  final_map <- final_map +
+    if(rast_type == "pop_density"){
+      draw_plot(rast_map_legend[[3]],
+                x = 0.48,
+                y = 0.03,
+                height = 0.09 ,
+                width = 0.1 - plot_margin) 
+    } else if(rast_type == "imp_surfaces") {
+      draw_plot(rast_map_legend[[3]],
+                x = 0.48,
+                y = 0.02,
+                height = 0.09 ,
+                width = 0.1 - plot_margin) 
+    }
+  
+  ggsave(outfile_path, final_map, width = width, height = height, dpi = viz_config_df$dpi, bg = viz_config_df$bg_col, units = "in")
 }
\ No newline at end of file

README.md

@@ -16,7 +16,7 @@ This project:
 
 ## Running the targets pipeline in R
 
-Clone the repo. In RStudio, run `library(targets)` and `tar_make()`. This will require `tidycensus` credentialing to run, see instructions below.
+Clone the repo. In RStudio, run `library(targets)` and `tar_make()`. This will require `tidycensus`, `sciencebase`, and `NASA Earthdata` credentialing to run, see instructions below.
 
 
 ## `tidycensus` set up
@@ -36,6 +36,28 @@ census_api_key("YOUR API KEY GOES HERE", install = TRUE)
 readRenviron("~/.Renviron")
 ```
 
+## ScienceBase login
+To access controlled ScienceBase items, an additional step must be completed to authenticate a ScienceBase session before the targets pipeline is run.
+
+1.  Ensure `.gitignore` includes a line for ".Renviron" to prevent credentials from being committed the the repository by git.
+2.  Run the entire contents of `00_config.R` either (1) interactively (i.e., line by line) or (2) by running `source("00_config.R")`.
+3.  You will be prompted to enter the token text by R/Rstudio. A browser window should pop-up directing you to <https://sciencebase.usgs.gov/manager/> to login.
+4.  Once logged in, copy your API token: click the account drop down button (person shaped icon) and click **Copy API Token**.
+5.  Paste the API token into the RStudio popup.
+
+Once the above steps are complete, `initialize_and_cache_sb()` will cache your username and token into the .Renviron file. From here, the pipeline should be able to use those cached credentials to re-initialize and refresh ScienceBase sessions in different R sessions for up to 10 hours.
+
+## NASA Earthdata
+Population density data must be manually downloaded from <https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-count-rev11/data-download>. You will be prompted to login to a NASA Earthdata account. If you do not have one, you will need to create one. Once logged in, use these selections:
+
+- Temporal: Single Year
+- FileFormat: GeoTiff
+- Resolution: 30 second
+- File Select: Year 2020
+
+Download the zip in the `1_process/in/` folder.
+
+
 ## Building the website locally
 
 Clone the repo. In the directory, run `npm install` to install the required modules. Once the dependencies have been installed, run `npm run dev` to run locally from your browser.

_targets.R

@@ -15,7 +15,10 @@ tar_option_set(packages = c('tidyverse',
                             'readr',
                             'cowplot',
                             'rmapshaper',
-                            'jsonlite'))
+                            'jsonlite',
+                            'terra',
+                            'tidyterra',
+                            'raster'))
 
 # Phase target makefiles
 source("0_config.R")