diff --git a/FrEDI/R/aggregate_impacts.R b/FrEDI/R/aggregate_impacts.R
index 3a6fe4c9..1bd66d68 100644
--- a/FrEDI/R/aggregate_impacts.R
+++ b/FrEDI/R/aggregate_impacts.R
@@ -47,9 +47,12 @@
 ### This function aggregates outputs produced by temperature binning
 aggregate_impacts <- function(
     data,             ### Data frame of outputs from temperature binning
-    aggLevels   = c("national", "modelaverage", "impactyear", "impacttype"),  ### Levels of aggregation
+    aggLevels   = c("national", "modelaverage", "impacttype", "impactyear"),  ### Levels of aggregation
     columns     = c("annual_impacts"), ### Columns to aggregate
-    groupByCols = c("sector", "variant", "impactYear", "impactType", "model_type", "model", "region"),
+    groupByCols = c("sector", "variant", "impactType", "impactYear") |>
+      c("region", "state", "postal") |>
+      c("model_type", "model") |>
+      c("includeaggregate", "sectorprimary"),
     silent      = TRUE
 ){
   ###### Defaults ######
@@ -59,7 +62,7 @@ aggregate_impacts <- function(
   yearCol0     <- c("year")
   # summaryCols  <- columns; rm(columns)
   ####### By State  ######
-  byState      <- c("state") %in% groupByCols
+  byState      <- c("state") %in% (data |> names())
   if(byState){stateCols0 <- c("state", "postal")} else{stateCols0 <- c()}
   popCol0      <- byState |> ifelse("state_pop", "reg_pop")
   # byState |> print()
@@ -76,7 +79,7 @@ aggregate_impacts <- function(
   ###### Aggregation Levels  ######
   ### Types of summarization to do: default
   # aggList0     <- c("national", "modelaverage", "impactyear", "impacttype", "all")
-  aggList0     <- c("national", "modelaverage", "impactyear", "impacttype")
+  aggList0     <- c("national", "modelaverage", "impacttype", "impactyear")
   null_aggLvls <- aggLevels |> is.null()
   aggLevels    <- aggLevels |> tolower()
   aggNone      <- "none" %in% aggLevels
@@ -145,7 +148,8 @@ aggregate_impacts <- function(
 
   ###### Grouping Columns  ######
   ### Use default group by columns if none specified...otherwise, check which are present
-  groupByCols0 <- c("sector", "variant", "impactType", "impactYear", "region") |> c(stateCols0)
+  groupByCols0 <- c("sector", "variant", "impactType", "impactYear")
+  groupByCols0 <- groupByCols0 |> c("region") |> c(stateCols0)
   groupByCols0 <- groupByCols0 |> c("model_type", "model")
   groupByCols0 <- groupByCols0 |> c("sectorprimary", "includeaggregate")
   if(groupByCols |> is.null()){groupByCols <- groupByCols0}
@@ -170,7 +174,7 @@ aggregate_impacts <- function(
     scalarCols  <- c("physScalar", "physAdj", "damageAdj", "econScalar", "econAdj", "econMultiplier") |> paste("Name")
     # scalarCols  <- c("physScalar", "physAdj", "damageAdj", "econScalar", "econAdj", "econMultiplier")
     # scalarCols  <- scalarCols0 |> map(~.x |> paste0(c(scalarSuffix0))) |> unlist()
-    dropCols    <- c("physicalmeasure") |> c(scalarCols)
+    dropCols    <- c("physicalmeasure") |> c(scalarCols) |> c(popCol0, "national_pop")
     isDropCol   <- groupByCols %in% dropCols
     hasDropCols <- isDropCol |> any()
     ### If hasDropCols
@@ -217,7 +221,7 @@ aggregate_impacts <- function(
   summaryCols <- summaryCols[isPresent0]
   rm(summaryCols0, nullSumCols, isPresent0, hasNaCols0)
 
-   ### Drop some columns from summary columns
+  ### Drop some columns from summary columns
   if(aggImpTypes){
     scalarCols  <- c("physScalar", "physAdj", "damageAdj", "econScalar", "econAdj", "econMultiplier") |> paste("Value")
     scalarCols  <- scalarCols |> c("physScalar", "econScalar", "physEconScalar")
@@ -229,7 +233,7 @@ aggregate_impacts <- function(
     if(hasDropCols){
       ### Drop levels
       summaryCols  <- summaryCols |> (function(y){y[!(y %in% dropCols)]})()
-    }
+    } ### End if(hasDropCols)
 
     ### If, message user
     if(hasDropCols & msgUser){
@@ -257,7 +261,7 @@ aggregate_impacts <- function(
     if(hasDropCols){
       ### Drop levels
       summaryCols  <- summaryCols |> (function(y){y[!(y %in% dropCols)]})()
-      } ### End if(hasDropCols)
+    } ### End if(hasDropCols)
 
     ###
     if(msgUser){
@@ -294,7 +298,8 @@ aggregate_impacts <- function(
   numCols0      <- numCols0    |> c("gdp_usd", "national_pop", "gdp_percap")
   numCols0      <- numCols0    |> c("driverValue") |> c(popCol0) |> c(yearCol0)
   numCols0      <- numCols0    |> unique()
-  data          <- data |> mutate_at(c(chrCols0), as.character)
+  # data          <- data |> mutate_at(c(chrCols0), as.character)
+  data          <- data |> mutate_at(c(numCols0), as.character)
   data          <- data |> mutate_at(c(numCols0), as.numeric  )
   rm(mutate0)
 
@@ -305,6 +310,7 @@ aggregate_impacts <- function(
   regPopCols    <- c("year", "region") |> c(stateCols0) |> c(popCol0) |> unique()
   natPopCols    <- c("year", "region") |> c("national_pop")
   driverCols    <- c("year", "model_type", "driverType", "driverUnit", "driverValue")
+
   ### Get names in names
   names0        <- data |> names()
   baseCols      <- baseCols   |> (function(y){y[(y %in% names0)]})()
@@ -313,6 +319,12 @@ aggregate_impacts <- function(
   driverCols    <- driverCols |> (function(y){y[(y %in% names0)]})()
   rm(names0)
 
+  # ### Add state_pop column
+  # if(aggNational){
+  #   groupByCols <- groupByCols |> c(regPopCols) |> unique()
+  #   summaryCols <- summaryCols |> (function(y){y[!(y %in% groupByCols)]})()
+  # } ### End if(aggNational)
+
   ### List of standardized columns
   standardCols  <- c(groupByCols, baseCols, regPopCols, natPopCols) |> unique()
   standardCols  <- standardCols |> c(driverCols, summaryCols) |> unique()
@@ -353,7 +365,7 @@ aggregate_impacts <- function(
   rm(data)
   # df_agg |> nrow() |> print(); df_agg |> head() |> glimpse()
 
-  ###### Impact Years ######
+  ###### ** Impact Years ######
   ### Separate into years after 2090 and before 2090
   if(aggImpYear){
     if(msgUser){msg0 (1) |> paste0("Interpolating between impact year estimates...") |> message()}
@@ -454,7 +466,8 @@ aggregate_impacts <- function(
   # paste0("Finished impact year interpolation: ", nrow(df_agg)) |> print(); df_agg |> head() |> glimpse()
   # "got here1" |> print()
 
-  ###### Model Averages ######
+  ###### ** Model Averages ######
+  # groupByCols |> print(); df_agg |> glimpse()
   ### Average values across models
   if(aveModels){
     modelAveMsg   <- "Calculating model averages..."
@@ -465,6 +478,7 @@ aggregate_impacts <- function(
     ### Group by columns
     group0        <- groupByCols |> (function(y){y[!(y %in% c("model", yearCol0))]})()
     group0        <- group0 |> c("year")
+    # group0 |> print()
     ### Separate model types
     df_gcm        <- df_agg |> filter(model_type |> tolower() == "gcm")
     df_agg        <- df_agg |> filter(model_type |> tolower() != "gcm")
@@ -500,23 +514,27 @@ aggregate_impacts <- function(
     df_agg <- df_gcm |> rbind(df_agg)
     rm(df_gcm, group0)
   } ### End if "model" %in% aggLevels
-  # paste0("Finished model aggregation: ", nrow(df_agg)) |> print(); df_agg |> head() |> glimpse()
+  # paste0("Finished model aggregation: ", nrow(df_agg)) |> print();
+  # df_agg |> glimpse()
   # "got here2" |> print()
 
-  ###### National Totals ######
+  ###### ** National Totals ######
   if(aggNational){
     if(msgUser){msg0 (1) |> paste0("Calculating national totals...") |> message()}
     ### Ungroup first
     df_agg      <- df_agg |> ungroup()
     ### Grouping columns
-    group0      <- groupByCols |> (function(y){y[!(y %in% c("region", stateCols0,popCol0, yearCol0))]})()
+    # group0      <- groupByCols |> (function(y){y[!(y %in% c("region", stateCols0, popCol0, yearCol0))]})()
+    group0      <- groupByCols |> (function(y){y[!(y %in% c("region", stateCols0, popCol0, yearCol0))]})()
     group0      <- group0 |> c("year")
     ### Calculate number of non missing values
     df_national <- df_agg |> (function(w){
       w |> mutate(not_isNA = 1 * (!(w[[summaryCol1]] |> is.na())))
     })()
     ### Group data, sum data, calculate averages, and drop NA column
-    sum0        <- summaryCols |> c(popCol0,"not_isNA")
+    # sum0        <- summaryCols |> c(popCol0) |> c("not_isNA")
+    sum0        <- summaryCols |> c("not_isNA")
+    # df_national |> glimpse()
     df_national <- df_national |>
       group_by_at(c(group0)) |>
       summarize_at(vars(sum0), sum, na.rm=T) |> ungroup()
@@ -531,15 +549,15 @@ aggregate_impacts <- function(
     df_national <- df_national |> select(-c("not_isNA"))
     df_national <- df_national |> mutate(region="National Total")
     ### Join with National Pop
-    #join0       <- natPopCols |> (function(y){y[!(y %in% c("national_pop"))]})()
-    #df_national <- df_national |> left_join(nationalPop, by = c(join0))
+    # join0       <- natPopCols |> (function(y){y[!(y %in% c("national_pop", popCol0))]})()
+    # df_national <- df_national |> left_join(nationalPop, by = c(join0))
     if(byState){
       df_national   <- df_national |> mutate(state ="All")
       df_national   <- df_national |> mutate(postal="US")
     } ### End if(byState)
 
     ### Add back into regional values and bind national population to impact types
-    #df_agg |> glimpse(); df_national |> glimpse()
+    # df_agg |> glimpse(); df_national |> glimpse()
     df_agg      <- df_agg |> rbind(df_national);
 
     ### Add national to total populations
@@ -551,7 +569,7 @@ aggregate_impacts <- function(
   # paste0("Finished national totals: ", nrow(df_agg)) |> print; df_agg |> head |> glimpse
   # "got here3" |> print()
 
-  ###### Impact Types ######
+  ###### ** Impact Types ######
   ### Summarize by Impact Type
   if(aggImpTypes){
     if(msgUser){msg0 (1) |> paste0("Summing across impact types...") |> message()}
@@ -613,12 +631,20 @@ aggregate_impacts <- function(
   df_agg        <- df_agg |> left_join(df_base      , by=c(join1))
   rm(names0, names1, join0, join1)
 
+  ###### Format Columns ######
   # ###### Reformat sectorprimary and includeaggregate, which were converted to character
   # mutate0         <- c("sectorprimary", "includeaggregate")
-  # mutate0         <- mutate0[mutate0 %in% names(df_return)]
+  # mutate0         <- mutate0[mutate0 %in% names(df_agg)]
   # doMutate        <- mutate0 |> length() > 0
-  # if(doMutate){df_return <- df_return |> mutate_at(c(mutate0), as.numeric)}
-  # if(doMutate){df_return <- df_return |> mutate_at(c(mutate0), as.numeric)}
+  # if(doMutate){df_agg <- df_agg |> mutate_at(c(mutate0), as.numeric)}
+  # if(doMutate){df_agg <- df_agg |> mutate_at(c(mutate0), as.numeric)}
+  mutate0       <- baseCols |> c(popCol0) |> c("driverValue")
+  mutate0       <- mutate0  |> c(summaryCols)
+  mutate0       <- mutate0  |> c("sectorprimary", "includeaggregate")
+  mutate0       <- mutate0  |> unique()
+  mutate0       <- mutate0  |> (function(y){y[y %in% (df_agg |> names())]})()
+  doMutate      <- (mutate0 |> length()) > 0
+  if(doMutate){df_agg <- df_agg |> mutate_at(c(mutate0), as.numeric)}
 
   ###### Order Columns ######
   ### Order the data frame and ungroup
diff --git a/FrEDI/R/run_fredi.R b/FrEDI/R/run_fredi.R
index 3a1a386a..d727a151 100644
--- a/FrEDI/R/run_fredi.R
+++ b/FrEDI/R/run_fredi.R
@@ -740,10 +740,10 @@ run_fredi <- function(
   ### For regular use (i.e., not impactYears), simplify the data: groupCols0
   if(requiresAgg){
     # df_results <- df_results |> aggregate_impacts(aggLevels=aggLevels, groupByCols=groupCols0)
-    # group0     <- groupCols0
+    group0     <- groupCols0
     # group0     <- select0 |> (function(x){x[!(x %in% driverCols0)]})()
     # select0 |> print(); df_results |> names() |> print()
-    group0     <- select0
+    # group0     <- select0
     df_results <- df_results |> aggregate_impacts(
       aggLevels   = aggLevels,
       groupByCols = group0,
diff --git a/FrEDI/R/utils_create_report_figures.R b/FrEDI/R/utils_create_report_figures.R
index eaffbcc0..f94462e2 100644
--- a/FrEDI/R/utils_create_report_figures.R
+++ b/FrEDI/R/utils_create_report_figures.R
@@ -1,46 +1,3 @@
-# ### Add names to list object
-# addListNames <- function(
-#     list0, ### List object
-#     names0 ### Names to give to list or data frame
-# ){
-#   names(list0) <- names0
-#   return(list0)
-# } ### End addListNames
-
-# ### This function makes it easier to get data objects from the sysdata.rda file
-# get_ciraDataObj <- function(
-#     x,       ### Object name
-#     listall  = FALSE,
-#     listName = "rDataList",
-#     pkg      = "FrEDI",
-#     lib.loc  = .libPaths()[1] ### Library path to look for packages
-#     ){
-#   ### Messaging
-#   msg0    <- "\t"
-#   ### Check if list name exists
-#   exists0 <- listName |> exists()
-#   ### If the listname exists in the name space, parse it
-#   ### Otherwise, grab it from a package name space
-#   if(exists0){new_x <- parse(text=listName) |> eval()}
-#   else       {
-#     ### Check if package & list name
-#     pkgList0    <- lib.loc |> installed.packages()
-#     pkgExists0  <- pkg %in% pkgList0
-#     if(!pkgExists0){
-#       msg0 |> paste0("Package doesn't exist...") |> message()
-#       msg0 |> paste0("Exiting...") |> message()
-#       return()
-#     } ### End if(!pkgExists0)
-#     else           {new_x <- getFromNamespace(listName, ns=pkg)}
-#   } ### End else(exists0)
-#
-#   ### Whether to list all items in data object or not
-#   if(listall) {return_x <- new_x |> names()}
-#   else        {return_x <- new_x[[x]]}
-#   ### Return
-#   return(return_x)
-# } ### End get_ciraDataObj
-
 ### Get column values from a tibble
 get_column_values <- function(
     df0,    ### Tibble
@@ -127,18 +84,25 @@ filter_years <- function(
   df0 <- df0 |> filter(year %in% years)
   ### Return
   return(df0)
-} ### End filter_years
+} ### End filter_years()
 
 ### Filter values
 ### Format values to specific number of decimal places
 format_values <- function(
     df0, ### data
-    cols0  = c("driverValue", "gdp_usd", "national_pop", "gdp_percap", "reg_pop", "annual_impacts"), ### Columns to format
+    byState = TRUE,
+    # cols0  = c("driverValue", "gdp_usd", "national_pop", "gdp_percap", "reg_pop", "annual_impacts"), ### Columns to format
     digits = 16
 ){
-  df0 <- df0 |> mutate_at(.vars=c(cols0), function(x){format(x, digits=digits)})
+  ### Pop columns
+  if(byState){popCols <- c("state", "postal")} else{c()}
+  popCol <- byState |> ifelse("state_pop", "reg_pop")
+  ### Columns
+  cols0 <- c("driverValue", "gdp_usd", "national_pop", "gdp_percap", popCol, "annual_impacts")
+  ### Mutate
+  df0 <- df0 |> mutate_at(vars(cols0), function(x){format(x, digits=digits)})
   return(df0)
-}
+} ### End format_values()
 
 ### Run CONUS scenarios
 create_constant_temp_scenario <- function(
@@ -189,7 +153,8 @@ create_constant_temp_scenario <- function(
 #### Get scenario inputs
 #### Get inputs list for a single scenario
 get_scenario_inputsList <- function(
-    df0   ### Data
+    df0,   ### Data
+    byState = TRUE
 ){
   ### df0 names
   names0  <- df0 |> names()
@@ -199,12 +164,15 @@ get_scenario_inputsList <- function(
   slr0    <- NULL
   gdp0    <- NULL
   pop0    <- NULL
+  ### Pop columns
+  if(byState){popCols <- c("state", "postal")} else{c()}
+  popCol <- byState |> ifelse("state_pop", "reg_pop")
   ### Columns for scenarios
   cTemp0  <- c("year", "temp_C_conus")
   cTemp1  <- c("year", "temp_C")
   cSlr    <- c("year", "slr_cm")
   cGdp    <- c("year", "gdp_usd")
-  cPop    <- c("year", "region", "reg_pop")
+  cPop    <- c("year", "region") |> c(popCols, popCol)
   ### Whether to create scenarios
   doTemp0 <- (cTemp0 %in% names0) |> all()
   doTemp1 <- (cTemp1 %in% names0) |> all()
@@ -217,28 +185,28 @@ get_scenario_inputsList <- function(
   if(doTemp){
     if(doTemp0){cTemp <- cTemp0}
     else       {cTemp <- cTemp1}
-    temp0   <- df0 |> select(c(all_of(cTemp)))
+    temp0   <- df0 |> select(all_of(cTemp))
     if(doTemp0){
-      temp0 <- temp0 |> rename_at(.vars=c("temp_C_conus"), ~c("temp_C"))
-    }
+      temp0 <- temp0 |> rename_at(vars("temp_C_conus"), ~c("temp_C"))
+    } ### End if(doTemp0)
     list0[["tempInput"]] <- temp0
     rm("temp0")
   } ### End if(doTemp)
 
   if(doSlr){
-    slr0   <- df0 |> select(c(all_of(cSlr)))
+    slr0   <- df0 |> select(all_of(cSlr))
     list0[["slrInput"]] <- slr0
     rm("slr0")
   } ### End if(doSlr)
 
   if(doGdp){
-    gdp0   <- df0 |> select(c(all_of(cGdp)))
+    gdp0   <- df0 |> select(all_of(cGdp))
     list0[["gdpInput"]] <- gdp0
     rm("gdp0")
   } ### End if(doGdp)
 
   if(doPop){
-    pop0   <- df0 |> select(c(all_of(cPop)))
+    pop0   <- df0 |> select(all_of(cPop))
     list0[["popInput"]] <- pop0
     rm("pop0")
   } ### End if(doPop)
@@ -275,11 +243,20 @@ agg_fredi_scenario <- function(
     joinCols  = c("year"),
     aggLevels = c("modelaverage", "national")
 ){
+  ### Pop cols
+  byState <- "state" %in% (df0 |> names())
+  if(byState){stateCols <- c("state", "postal")} else{stateCols <- c()}
+  popCol  <- byState |> ifelse("state_pop", "reg_pop")
   ### Filter to grouping columns
   drop0   <- scenCols[!(scenCols %in% joinCols)]
   ### Run FrEDI
-  group0  <- c("sector", "variant", "impactYear", "impactType", "model_type", "model", "region") |> c(drop0)
+  group0  <- c("sector", "variant", "impactType", "impactYear")
+  group0  <- group0 |> c("region", stateCols)
+  group0  <- group0 |> c("model_type", "model")
+  group0  <- group0 |> c("sectorprimary", "includeaggregate")
+  group0  <- group0 |> c(drop0)
   df0     <- df0 |> FrEDI::aggregate_impacts(aggLevels = aggLevels, groupByCols = group0)
+  # df0     <- df0 |> FrEDI::aggregate_impacts(aggLevels = aggLevels)
   ### Return
   return(df0)
 } ### End agg_fredi_scenario
@@ -301,22 +278,30 @@ run_scenario <- function(
   rm("df0")
   ### Run FrEDI
   if(fredi){
-    df_x0     <- df_x0 |> run_fredi_scenario(
+    df_x0 <- df_x0 |> run_fredi_scenario(
       sectors   = sectors,
       scenCols  = scenCols,
       joinCols  = joinCols
     ) ### End run_fredi_scenario
   } ### End if(fredi)
+  # "got here1" |> print(); df_x0 |> glimpse()
+
   ### Aggregate FrEDI
   agg0      <- !("none" %in% aggLevels)
   # agg0 |> print()
   if(agg0){
-    df_x0     <- df_x0 |> agg_fredi_scenario(
+    # "got here1" |> print()
+    df_x0 <- df_x0 |> agg_fredi_scenario(
       scenCols  = scenCols,
       joinCols  = joinCols,
       aggLevels = aggLevels
     ) ### End run_fredi_scenario
   } ### End if(agg0)
+  # "got here2" |> print(); df_x0 |> glimpse()
+
+  ### Format other values
+  mutate0   <- c("temp_C_conus", "temp_C_global", "slr_cm")
+  df_x0     <- df_x0 |> mutate_at(vars(mutate0), as.numeric)
 
   ### Return
   return(df_x0)
@@ -352,7 +337,8 @@ run_scenarios <- function(
     return(df_x)
   }) ### End function(.x), walk
   ### Bind values into a list
-  df0    <- list0 %>% (function(x){do.call(rbind, x)})
+  # df0    <- list0 %>% (function(x){do.call(rbind, x)})
+  df0    <- list0 |> bind_rows()
 
   ### Return
   return(df0)
@@ -407,7 +393,8 @@ sum_impacts_byDoW <- function(
     return(df_z)
   })
   ### Bind together
-  df0        <- list0 %>% (function(x){do.call(rbind, x)})
+  # df0        <- list0 %>% (function(x){do.call(rbind, x)})
+  df0        <- list0 |> bind_rows()
   rm(list0)
   ### Adjust values
   df0[[adjCol]] <- df0[["annual_impacts"]] * adjVal
@@ -468,7 +455,8 @@ sum_impacts_byDoW_years <- function(
     return(df_x)
   }) ### End walk
   ### Bind values together
-  df0 <- list0 %>% (function(x){do.call(rbind, x)})
+  # df0 <- list0 %>% (function(x){do.call(rbind, x)})
+  df0 <- list0 |> bind_rows()
   rm(list0)
   ### Convert to tibble
   df0 <- df0 |> as_tibble()
@@ -488,19 +476,19 @@ get_fig7_slrDataObj <- function(
   ### Sector Info
   ### Variant Info
   ### Model Info
-  dfSectors <- "co_sectors"  |> get_ciraDataObj()
-  dfVariant <- "co_variants" |> get_ciraDataObj()
-  slrRef    <- "co_models"   |> get_ciraDataObj()
+  dfSectors <- "co_sectors"  |> get_frediDataObj()
+  dfVariant <- "co_variants" |> get_frediDataObj()
+  slrRef    <- "co_models"   |> get_frediDataObj()
 
   ### SLR Driver values
   ### SLR Scaled impct values
-  if(drivers){slrCm  <- "slr_cm"     |> get_ciraDataObj()}
-  if(impacts){slrImp <- "slrImpacts" |> get_ciraDataObj()}
+  if(drivers){slrCm  <- "slr_cm"     |> get_frediDataObj()}
+  if(impacts){slrImp <- "slrImpacts" |> get_frediDataObj(listSub="stateData")}
 
   ###### SLR Models ######
   ### Format SLR Models
   slrRef    <- slrRef |> filter(modelType=="slr")
-  slrRef    <- slrRef |> rename_at(.vars=c("model_label"), ~c("model"))
+  slrRef    <- slrRef |> rename_at(vars("model_label"), ~c("model"))
 
   ###### Levels & Labels ######
   ### Initial levels & labels
@@ -513,7 +501,7 @@ get_fig7_slrDataObj <- function(
   # slrLevels |> print(); slrLabels |> print()
   ### Vector of model labels and number of models
   c_slrs    <- slrLabels
-  n_slrs    <- c_slrs %>% length
+  n_slrs    <- c_slrs |> length()
 
   ###### Sectors Data ######
   ### Format Sectors data
@@ -845,17 +833,19 @@ plot_DoW <- function(
   df_types   <- tibble()
   if(do_gcm){
     df_gcm   <- "GCM" %>%
-      map(function(.x){tibble(type=.x, year=years0, label=.x |> paste0("_", years0))}) %>%
-      (function(y){do.call(rbind, y)})
+      # map(function(.x){tibble(type=.x, year=years0, label=.x |> paste0("_", years0))}) %>%
+      # (function(y){do.call(rbind, y)})
+      map(function(.x){tibble(type=.x, year=years0, label=.x |> paste0("_", years0))}) |>
+      bind_rows()
     df_types <- df_types |> rbind(df_gcm)
     rm(df_gcm)
-  }
+  } ### if(do_gcm)
   ### SLR data
   if(do_slr){
     df_slr   <- tibble(type="SLR", year="all", label="SLR" |> paste0("_", "all"))
     df_types <- df_types |> rbind(df_slr)
     rm(df_slr)
-  }
+  } ### if(do_slr)
   # "got here" |> print()
   # df_types |> glimpse()
 
@@ -864,7 +854,7 @@ plot_DoW <- function(
   ### Initialize list
   list0      <- pList0 %>% pmap(function(x1, x2){
     x1 |> paste0("_", x2) |> print()
-    plot_y <- plot_DoW_by_modelYear(
+    plot_y   <- plot_DoW_by_modelYear(
       df0        = df0, ### Data (e.g., output from sum_impactsByDegree)
       type0      = x1,  ### Model type: GCM or SLR
       year0      = x2,
@@ -884,7 +874,7 @@ plot_DoW <- function(
   ### Add list names
   # list0 |> print()
   labels0    <- df_types[["label"]]
-  list0      <- list0 |> addListNames(labels0)
+  list0      <- list0 |> set_names(labels0)
 
   ### Return
   return(list0)
@@ -922,10 +912,12 @@ plot_DoW_by_sector <- function(
         df1      <- df0 |> filter(model_type=="GCM")
         sectors0 <- df1[["sector"]] |> unique()
         df_x     <- sectors0 |> map(function(.y){tibble(type=.x, sector=.y, year=years, label=.y |> paste0("_", years))})
-        df_x     <- df_x %>% (function(y){do.call(rbind, y)})
+        # df_x     <- df_x %>% (function(y){do.call(rbind, y)})
+        df_x     <- df_x |> bind_rows()
         return(df_x)
       })
-    df_gcm   <- df_gcm %>% (function(y){do.call(rbind, y)})
+    # df_gcm   <- df_gcm %>% (function(y){do.call(rbind, y)})
+    df_gcm   <- df_gcm |> bind_rows()
     df_types <- df_types |> rbind(df_gcm)
     rm(df_gcm)
   } ### End if(do_gcm)
@@ -972,12 +964,12 @@ plot_DoW_by_sector <- function(
     })
     ### Add names
     labels_x <- types_x[["label"]]
-    list_x   <- list_x |> addListNames(labels_x)
+    list_x   <- list_x |> set_names(labels_x)
     ### Return
     return(list_x)
   })
   ### Add names
-  list0   <- list0 |> addListNames(models)
+  list0   <- list0 |> set_names(models)
   ### Return
   return(list0)
 } ### End plot_DoW_by_sector
@@ -1010,10 +1002,10 @@ plot_slr_scenarios <- function(
     scale_x_continuous("Year") +
     scale_y_continuous("GMSL (cm)")
 
-  plot0      <- plot0 +
-    theme(panel.background = element_rect(fill="white")) +
-    theme(panel.grid       = element_line(color="lightgrey")) +
-    theme(axis.line        = element_line(color="darkgrey"))
+  # plot0      <- plot0 +
+  #   theme(panel.background = element_rect(fill="white")) +
+  #   theme(panel.grid       = element_line(color="lightgrey")) +
+  #   theme(axis.line        = element_line(color="darkgrey"))
 
   plot0      <- plot0 +
     ggtitle(title0, subTitle0) +
diff --git a/FrEDI/R/utils_plot_DOW_byImpactType.R b/FrEDI/R/utils_plot_DOW_byImpactType.R
index c73af167..82ec1982 100644
--- a/FrEDI/R/utils_plot_DOW_byImpactType.R
+++ b/FrEDI/R/utils_plot_DOW_byImpactType.R
@@ -157,39 +157,46 @@ plot_DOW_byImpactType <- function(
   # subtitle0 <- variant0
 
   ###### Create the plot ######
-  plot0  <- df0 |> ggplot(aes(x=.data[[xCol]], y=.data[[yCol]]))
+  plot0     <- df0 |> ggplot(aes(x=.data[[xCol]], y=.data[[yCol]]))
 
   ### Add Geoms
-  plot0  <- plot0 + geom_line (aes(color = model))
+  plot0     <- plot0 + geom_line (aes(color = model))
   if(do_slr){df_points0 <- df0 |> filter(year %in% x_breaks)}
   else      {df_points0 <- df0}
-  plot0  <- plot0 + geom_point(data=df_points0, aes(color = model, shape=model))
+  plot0     <- plot0 + geom_point(data=df_points0, aes(color = model, shape=model))
 
   ### Add Scales
-  plot0  <- plot0 + scale_color_discrete(lgdLbl)
-  plot0  <- plot0 + scale_shape_discrete(lgdLbl)
+  # plot0     <- plot0 + scale_shape_discrete(lgdTitle0)
+  shapeLvls <- df0[["model"]] |> unique() |> sort()
+  numShapes <- shapeLvls |> length()
+  shapeVals <- c(1:numShapes)
+  # shapeLvls |> print()
+  # plot0     <- plot0 + scale_shape_discrete(lgdLbl)
+  plot0     <- plot0 + scale_shape_manual(lgdLbl, breaks=shapeLvls, values=shapeVals)
+  plot0     <- plot0 + scale_color_discrete(lgdLbl)
+  # plot0     <- plot0 + scale_shape_discrete(lgdLbl)
 
   ###### Adjust legend title ######
   if(hasLgdPos){plot0  <- plot0  + guides(color = guide_legend(title.position = lgdPos))}
 
   ###### Add themes and title ######
-  plot0  <- plot0 + ggtitle(title0, subtitle0)
+  plot0     <- plot0 + ggtitle(title0, subtitle0)
 
   ###### Add scales ######
-  plot0  <- plot0 + scale_x_continuous(xTitle, breaks = x_breaks, limits = x_limits)
-  # plot0    <- plot0 + scale_y_continuous(yTitle)
-  plot0  <- plot0 + scale_y_continuous(y_label)
+  plot0     <- plot0 + scale_x_continuous(xTitle, breaks = x_breaks, limits = x_limits)
+  # plot0     <- plot0 + scale_y_continuous(yTitle)
+  plot0     <- plot0 + scale_y_continuous(y_label)
 
   ###### Adjust Appearance ######
-  # plot0  <- plot0 + theme(panel.background = element_rect(fill="white"))
-  # plot0  <- plot0 + theme(panel.grid = element_line(color="lightgrey"))
-  # plot0  <- plot0 + theme(plot.title    = element_text(hjust = 0.5, size=11))
-  plot0  <- plot0 + theme(plot.title    = element_text(hjust = 0.5, size=11))
-  plot0  <- plot0 + theme(plot.subtitle = element_text(hjust = 0.5, size=10))
-  # plot0  <- plot0 + theme(axis.title.x  = element_text(hjust = 0.5, size=9, color="white"))
-  plot0  <- plot0 + theme(axis.title.x  = element_text(hjust = 0.5, size=9))
-  plot0  <- plot0 + theme(axis.title.y  = element_text(hjust = 0.5, size=9))
-  plot0  <- plot0 + theme(legend.position = "bottom")
+  # plot0     <- plot0 + theme(panel.background = element_rect(fill="white"))
+  # plot0     <- plot0 + theme(panel.grid = element_line(color="lightgrey"))
+  # plot0     <- plot0 + theme(plot.title    = element_text(hjust = 0.5, size=11))
+  plot0     <- plot0 + theme(plot.title    = element_text(hjust = 0.5, size=11))
+  plot0     <- plot0 + theme(plot.subtitle = element_text(hjust = 0.5, size=10))
+  # plot0     <- plot0 + theme(axis.title.x  = element_text(hjust = 0.5, size=9, color="white"))
+  plot0     <- plot0 + theme(axis.title.x  = element_text(hjust = 0.5, size=9))
+  plot0     <- plot0 + theme(axis.title.y  = element_text(hjust = 0.5, size=9))
+  plot0     <- plot0 + theme(legend.position = "bottom")
 
   ###### Plot Index #####
   ###### If plotIndex, remove some plot elements
diff --git a/FrEDI/R/utils_plot_DOW_byImpactTypes.R b/FrEDI/R/utils_plot_DOW_byImpactTypes.R
index 6b35034e..ca384930 100644
--- a/FrEDI/R/utils_plot_DOW_byImpactTypes.R
+++ b/FrEDI/R/utils_plot_DOW_byImpactTypes.R
@@ -220,7 +220,7 @@ plot_DOW_byImpactTypes <- function(
       })
 
       ### Name the plots
-      listVars_j  <- listVars_j |> addListNames(c_variants)
+      listVars_j  <- listVars_j |> set_names(c_variants)
       # return(listVars_j)
       # "got here1..." |> print()
 
@@ -245,7 +245,7 @@ plot_DOW_byImpactTypes <- function(
 
     ### Name the plots
     # listTypes_i |> length() |> print(); c_impTypes |> print()
-    listTypes_i <- listTypes_i |> addListNames(c_impTypes)
+    listTypes_i <- listTypes_i |> set_names(c_impTypes)
     # "got here3..." |> print()
     # return(listTypes_i)
 
@@ -284,7 +284,7 @@ plot_DOW_byImpactTypes <- function(
     return(plotGrid_i)
   })
   ### Name the plots
-  listYears0 <- listYears0 |> addListNames(c_impYears)
+  listYears0 <- listYears0 |> set_names(c_impYears)
 
   ###### Return ######
   ### Return the plot
diff --git a/FrEDI/R/utils_plot_DOW_byModelType.R b/FrEDI/R/utils_plot_DOW_byModelType.R
index 014e674d..16bfa608 100644
--- a/FrEDI/R/utils_plot_DOW_byModelType.R
+++ b/FrEDI/R/utils_plot_DOW_byModelType.R
@@ -143,7 +143,7 @@ plot_DOW_byModelType <- function(
 
   ### Add list names
   # # x |> length() |> print()
-  plotList_x <- plotList_x |> addListNames(c_sectors)
+  plotList_x <- plotList_x |> set_names(c_sectors)
 
   ###### Get Reference Plot ######
   refPlot_x  <- c_sectors[1] |> plot_DOW_bySector(
diff --git a/FrEDI/R/utils_plot_DOW_bySector.R b/FrEDI/R/utils_plot_DOW_bySector.R
index 3f3ff6ae..7b22e64f 100644
--- a/FrEDI/R/utils_plot_DOW_bySector.R
+++ b/FrEDI/R/utils_plot_DOW_bySector.R
@@ -84,17 +84,25 @@ plot_DOW_bySector <- function(
   if(!hasMUnits){mUnit0 <- "cm"}
 
   ###### Create the plot ######
-  # df0 %>% names() %>% print()
-  # df0 %>% glimpse()
+  # df0 |> names() |> print()
+  # df0 |> glimpse()
   plot0      <- df0 |> ggplot(aes(x=.data[[xCol]], y=.data[[yCol]]))
+  # plot0      <- df0 |> ggplot(aes(x=.data[[xCol]], y=.data[[yCol]], group=interaction(sector, model)))
 
   ### Add Geoms
   plot0      <- plot0 + geom_line (aes(color = model))
+  # plot0      <- plot0 + geom_point(aes(color = model))
   plot0      <- plot0 + geom_point(aes(color = model, shape=model))
 
   ### Add Scales
   plot0      <- plot0 + scale_color_discrete(lgdTitle0)
-  plot0      <- plot0 + scale_shape_discrete(lgdTitle0)
+  # plot0      <- plot0 + scale_shape_discrete(lgdTitle0)
+  shapeLvls  <- df0[["model"]] |> unique() |> sort()
+  numShapes  <- shapeLvls |> length()
+  shapeVals  <- c(1:numShapes)
+  # shapeLvls |> print()
+  # plot0      <- plot0 + scale_shape_discrete(lgdTitle0)
+  plot0      <- plot0 + scale_shape_manual(lgdTitle0, breaks=shapeLvls, values=shapeVals)
   plot0      <- plot0 + scale_x_continuous(xTitle0, limits = x_limits, breaks = x_breaks)
   plot0      <- plot0 + scale_y_continuous(yTitle0, limits = y_limits, breaks = y_breaks)
 
diff --git a/FrEDI/scripts/create_DoW_results.R b/FrEDI/scripts/create_DoW_results.R
index 8fe5f851..656902b2 100644
--- a/FrEDI/scripts/create_DoW_results.R
+++ b/FrEDI/scripts/create_DoW_results.R
@@ -9,18 +9,32 @@ require(ggpubr)
 ###### create_DoW_results ######
 create_DoW_results <- function(
     sectors  = FrEDI::get_sectorInfo(), ### Which sectors
-    gcmYears = c(2090), ### Which years to report on for GCM sectors
+    gcmYears = c(2090),   ### Which years to report on for GCM sectors
     slrYears = c(2050, 2090), ### Which years to report on for SLR sectors
-    silent   = TRUE,  ### Degree of messaging
-    testing  = FALSE, ### Whether to print out extra diagnostic values
-    aggOnly  = TRUE, ### Whether to only include sectors for which "includeaggregate==1" in Fig 7 plots
+    byState  = TRUE,       ### Whether values are by state or just by region
+    silent   = TRUE,      ### Degree of messaging
+    testing  = FALSE,     ### Whether to print out extra diagnostic values
+    aggOnly  = TRUE,      ### Whether to only include sectors for which "includeaggregate==1" in Fig 7 plots
     loadCode = "project", ### Whether to load code as source or devtools
-    fpath    = "." ,  ### Path to main FrEDI directory to load code from if loadCode == "project" or loadCode == "package"
-    saveFile = FALSE, ### Save file
+    fpath    = "." ,      ### Path to main FrEDI directory to load code from if loadCode == "project" or loadCode == "package"
+    saveFile = FALSE,     ### Save file
     outPath  = "." |> file.path("report_figures"),  ### Path to save results if saveFile == TRUE
-    img_dev  = "pdf", ### Image device if saveFile == TRUE
-    return   = TRUE   ### Whether to return list object
+    img_dev  = "pdf",     ### Image device if saveFile == TRUE
+    return   = TRUE       ### Whether to return list object
 ){
+  # sectors  = FrEDI::get_sectorInfo() ### Which sectors
+  # gcmYears = c(2090) ### Which years to report on for GCM sectors
+  # slrYears = c(2050, 2090) ### Which years to report on for SLR sectors
+  # silent   = TRUE    ### Degree of messaging
+  # testing  = TRUE    ### Whether to print out extra diagnostic values
+  # byState  = TRUE    ### Whether values are by state or just by region
+  # aggOnly  = TRUE    ### Whether to only include sectors for which "includeaggregate==1" in Fig 7 plots
+  # loadCode = "project" ### Whether to load code as source or devtools
+  # fpath    = "."     ### Path to main FrEDI directory to load code from if loadCode == "project" or loadCode == "package"
+  # saveFile = FALSE   ### Save file
+  # outPath  = "." |> file.path("report_figures")  ### Path to save results if saveFile == TRUE
+  # img_dev  = "pdf"   ### Image device if saveFile == TRUE
+  # return   = TRUE    ### Whether to return list object
   ###### Initial values ######
   ### Messaging
   do_msg          <- !silent
@@ -100,12 +114,15 @@ create_DoW_results <- function(
   ###### ** Constants ######
   ### Numeric columns: Specify so that we can print out the associated data
   ### Number of digits to format
-  c_numVars      <- c("driverValue", "gdp_usd", "national_pop", "gdp_percap", "reg_pop", "annual_impacts")
+  c_popCol       <- byState |> ifelse("state_pop", "reg_pop")
+  # c_numVars      <- c("driverValue", "gdp_usd", "national_pop", "gdp_percap", "reg_pop", "annual_impacts")
+  c_numVars      <- c("driverValue", "gdp_usd", "national_pop", "gdp_percap") |> c(c_popCol) |> c("annual_impacts")
   ### Integer temperatures: data frame of inputs
   conusPrefix0   <- "Other_Integer"
   globalPrefix0  <- "preI_global"
   ### Temperatures
-  c_conusTemps   <- 0:7
+  # c_conusTemps   <- 0:7
+  c_conusTemps   <- 0:10
   c_globalTemps  <- c(1.487, 2.198)
   ### Numbers of scenarios
   n_conusTemps   <- c_conusTemps |> length()
@@ -136,7 +153,9 @@ create_DoW_results <- function(
     x = df_scenarios[["temp_C"  ]],
     y = df_scenarios[["tempType"]],
     z = df_scenarios[["prefix"  ]]
-    ) |> pmap(function(x, y, z){
+  )
+  ### Create constant temp scenarios
+  inputs_df_int  <- inputs_df_int |> pmap(function(x, y, z){
     create_constant_temp_scenario(
       temp0   = x,
       type0   = y,
@@ -162,6 +181,7 @@ create_DoW_results <- function(
     scenCols  = c("scenario", "year", "temp_C_conus", "temp_C_global", "slr_cm"),
     joinCols  = c("year")
   )
+
   ### Glimpse results
   if(return0) resultsList[["df_int_byType"]] <- df_int_byType
   if(testing) df_int_byType |> glimpse()
@@ -182,7 +202,7 @@ create_DoW_results <- function(
   ###### ** Result Totals ######
   if(testing|do_msg) "Aggregating integer scenario results..." |> message()
   #### Aggregate Impact Types, Impact Years
-  df_int_totals  <- df_int_byType %>% run_scenarios(
+  df_int_totals  <- df_int_byType |> run_scenarios(
     col0      = "scenario",
     fredi     = FALSE,
     aggLevels = c("impactyear", "impacttype"),
@@ -241,15 +261,15 @@ create_DoW_results <- function(
   # codePath  |> loadCustomFunctions()
   if(testing|do_msg) "Plotting GCM results by sector, degree of warming (DOW)..." |> message()
   plots_dow_gcm  <- sum_gcm_totals |> plot_DoW(
-      types0  = c("GCM"), ### Model type: GCM or SLR
-      years   = gcmYears,
-      xCol    = "driverValue",
-      yCol    = "annual_impacts",
-      thresh0 = breakChars
-    )
+    types0  = c("GCM"), ### Model type: GCM or SLR
+    years   = gcmYears,
+    xCol    = "driverValue",
+    yCol    = "annual_impacts",
+    thresh0 = breakChars
+  )
   ### Glimpse
+  if(testing) plots_dow_gcm[["GCM_2090"]] |> print()
   if(return0) resultsList[["plots_dow_gcm"]] <- plots_dow_gcm
-  if(testing) plots_dow_gcm[["GCM_2010"]] |> print()
   ### Save
   # codePath  |> loadCustomFunctions()
   if(saveFile){
@@ -280,8 +300,8 @@ create_DoW_results <- function(
     silent      = TRUE
   )
   ### Glimpse
-  if(return0) resultsList[["sum_gcm_byType"]] <- sum_gcm_byType
   if(testing) sum_gcm_byType |> glimpse()
+  if(return0) resultsList[["sum_gcm_byType"]] <- sum_gcm_byType
   ### Save summary table
   if(saveFile){
     if(do_msg) paste0("Saving summary of GCM results by sector, impact type, degree of warming...") |> message()
@@ -294,13 +314,14 @@ create_DoW_results <- function(
   if(testing|do_msg) "Plotting GCM results by sector, impact type, degree of warming (DOW)..." |> message()
   plots_gcm_byType <- sum_gcm_byType |>
     # filter(sector %in% c_sectorNames[c(10)]) |>
+    filter(!(sector %in% c("Roads"))) |>
     plot_DoW_by_sector(
       models  = c("GCM"),
       yCol    = "annual_impacts"
     )
   ### Glimpse
-  if(return0) resultsList[["plots_gcm_byType"]] <- plots_gcm_byType
   if(testing) plots_gcm_byType$GCM$`Extreme Temperature_2010`[["2010"]] |> print()
+  if(return0) resultsList[["plots_gcm_byType"]] <- plots_gcm_byType
   ### Save
   if(saveFile){
     if(do_msg) paste0("Saving plots of GCM results by  sector, impact type, degree of warming...") |> message()