Maize_Spread_Processes-copy.Rmd

---
title: "Spread of Maize Agriculture to the Southwestern United States"
author: "Andrew Gillreath-Brown"
date: "`r format(Sys.time(), '%d %B %Y')`"
html_document:
  toc: yes
toc_depth: 2
toc_float: yes
number_sections: no
code_folding: hide
output:
  html_document:
  df_print: paged
mainfont: Calibri
editor_options:
  chunk_output_type: console
keep_md: yes
---

```{r setup, include = FALSE, results = 'hide', echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide'}

# devtools::install()
# devtools::load_all()
library(cropSpreadR)
require(ggh4x)
library(cropDiffusionR)
library(raster)
library(tidyverse)

```

```{r load_data, include = FALSE, results = 'hide', echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide'} 

# Load the temperature anomaly reconstruction from paleoMAT package.
temp_anom <-
  readr::read_csv(
    "https://raw.githubusercontent.com/Archaeo-Programmer/paleoMAT/master/vignettes/tables/supp_table_s2.csv",
    show_col_types = FALSE
  ) %>%
  dplyr::rename(
    year = 1,
    anom = 2,
    ci = 3,
    se = 4
  )

# Load location data for each sample. Site locations are not saved in this repository, as we want to ensure that archaeological site locations are protected.
sf::sf_use_s2(FALSE)

MDB <- cropDiffusionR::maizeDB %>%
  dplyr::left_join(.,
                   readr::read_csv(
                     stringr::str_remove(
                       here::here("Paper 2/Maize_Database_Locations_DONT_DISTRIBUTE.csv"),
                       "cropSpreadR/"
                     ),
                     col_types = readr::cols()
                   ),
                   by = "LabID") %>%
  dplyr::mutate(Long1 = Long, Lat1 = Lat) %>%
  # Convert to spatial object
  sf::st_as_sf(coords = c("Long1", "Lat1"),
               crs = 4326) %>%
  sf::st_make_valid() %>%
  # Join to counties data to get county and state information.
  sf::st_join(cropSpreadR::SWUS_counties) %>%
  dplyr::mutate(county_st = ifelse(!is.na(REGION) &
                                     !is.na(NAME), paste(NAME, REGION, sep = ", "), NA)) %>%
  sf::st_drop_geometry() %>% 
  suppressMessages()

# Get modern annual accumulated growing degree days at each sample location. 
# For Mexico, we use WorldClim 30 sec. (~1 km) for GDD data.
MDB_gdd_wc <- MDB %>%
  sf::st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  dplyr::select(LabID) %>%
  cropDiffusionR::extract_GDD(gdd_raster = cropDiffusionR::WorldClim_annual_gdd,
                              sites = .)

# For USA, we use PRISM (~800 meters) for GDD data.
MDB_gdd_prism <- MDB %>%
  sf::st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  dplyr::select(LabID) %>%
  cropDiffusionR::extract_GDD(gdd_raster = cropSpreadR::PRISM_annual_gdd,  sites = .)

# Join the GDD together.
MDB_gdd <-
  dplyr::left_join(MDB_gdd_prism, MDB_gdd_wc, by = "LabID") %>%
  dplyr::mutate(GDD = dplyr::coalesce(GDD.x, GDD.y)) %>%
  dplyr::select(-c(GDD.x, GDD.y))

# Add GDD and elevation to maize dataset.
MDB_env <- MDB %>%
  dplyr::select(-Elevation_masl) %>%
  dplyr::left_join(., MDB_gdd, by = "LabID") %>%
  # Add elevation to dataset.
  dplyr::left_join(., cropSpreadR::MDB_elev, by = "LabID")

# Get sites that are on the northern maize frontier through time from cropDiffusionR package for the SWUS.
temp <- tempfile()
download.file(
  "https://raw.githubusercontent.com/Archaeo-Programmer/cropDiffusionR/master/vignettes/output/maize_frontier_swus.rds",
  temp
)
maize_frontier_swus <-
  readr::read_rds(temp, "maize_frontier_swus.rds")
unlink(temp)

maize_frontier_swus <- MDB %>%
  dplyr::filter(LabID %in% maize_frontier_swus$LabID) %>%
  dplyr::mutate(Long1 = Long, Lat1 = Lat) %>%
  # Convert to spatial object
  sf::st_as_sf(coords = c("Long1", "Lat1"),
               crs = 4326) %>%
  sf::st_make_valid()

# Get sites that are on the northern maize frontier through time from cropDiffusionR package for Mexico and SWUS.
temp <- tempfile()
download.file(
  "https://raw.githubusercontent.com/Archaeo-Programmer/cropDiffusionR/master/vignettes/output/maize_frontier_mexico.rds",
  temp
)
maize_frontier_mexico <-
  readr::read_rds(temp, "maize_frontier_mexico.rds")
unlink(temp)

maize_frontier_mexico <- MDB_env %>%
  dplyr::filter(LabID %in% maize_frontier_mexico$LabID) %>%
  dplyr::mutate(Long1 = Long, Lat1 = Lat) %>%
  # Convert to spatial object
  sf::st_as_sf(coords = c("Long1", "Lat1"),
               crs = 4326) %>%
  sf::st_make_valid() %>% 
  dplyr::mutate(centroid_long = ifelse(is.na(centroid_long), Long, centroid_long),
                centroid_lat = ifelse(is.na(centroid_lat), Lat, centroid_lat))

```

```{r Figure_1, fig.cap="Figure 1. Frequency distribution of directly dated maize macrofossils in our maize database for sites that are north of 30° latitude and older than AD 1500 (n = 985), distinguishing between the earliest maize sample at each site (n = 275), and all other maize samples fitting the criteria above (n = 710). Two of the earliest maize sites, Las Capas and Old Corn, are individually labeled.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

maize_earliest_distribution <- MDB_env %>%
  dplyr::filter(Lat > 30 & Date < 1500) %>%
  dplyr::group_by(SiteIDName) %>%
  dplyr::slice_min(Date, with_ties = FALSE) %>%
  dplyr::mutate(sample = "Earliest Maize Subset")

full_distribution <- MDB_env %>%
  dplyr::filter(Lat > 30 & Date < 1500) %>%
  dplyr::mutate(sample = "Full Maize Database") %>%
  dplyr::filter(!LabID %in% maize_earliest_distribution$LabID)

swus_date_distribution_stacked <- maize_earliest_distribution %>%
  dplyr::bind_rows(., full_distribution) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(sample = factor(sample, levels = c(
    "Full Maize Database", "Earliest Maize Subset"
  ))) %>%
  ggplot(aes(x = Date, fill = sample)) +
  geom_histogram(
    binwidth = 100,
    color = "#000000",
    position = "stack",
    center = -3650
  ) +
  scale_fill_manual(values = c("cornflowerblue", "#DDCC77")) +
  xlab("Year BC/AD") +
  ylab("Count") +
  scale_x_continuous(breaks = seq(-4000, 1500, 500),
                     limits = c(-4000, 1500)) +
  scale_y_continuous(breaks = seq(0, 90, 10), limits = c(0, 90)) +
  theme_classic() +
  annotate(
    geom = "text",
    x = -3550,
    y = 15,
    label = "Las Capas",
    hjust = 0,
    angle = 45,
    size = 6.75
  ) +
  annotate(
    "segment",
    x = -3650,
    xend = -3550,
    y = 1,
    yend = 13,
    colour = "black"
  ) +
  annotate(
    "segment",
    x = -3450,
    xend = -3550,
    y = 1,
    yend = 13,
    colour = "black"
  ) +
  # Annotation for Old Corn.
  annotate(
    geom = "text",
    x = -2150,
    y = 15,
    label = "Old Corn",
    hjust = 0,
    angle = 45,
    size = 6.75
  ) +
  annotate(
    "segment",
    x = -2250,
    xend = -2150,
    y = 1,
    yend = 13,
    colour = "black"
  ) +
  theme(
    legend.title = element_blank(),
    legend.text = element_text(size = 14,
                               family = "Helvetica"),
    panel.spacing = unit(1, "lines"),
    strip.text.x = element_text(
      size = 16,
      family = "Helvetica",
      face = "bold"
    ),
    axis.ticks.length = unit(0.125, "inch"),
    #ggh4x.axis.ticks.length.minor = rel(0.5),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_text(
      size = 20,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 5,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 5,
        r = 10,
        b = 10,
        l = 10
      )
    ),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"),
    plot.margin = margin(
      t = 5,
      r = 25,
      b = 5,
      l = 5
    )
  )

swus_date_distribution_stacked

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_1.png"),
#   plot = swus_date_distribution_stacked,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

```

```{r Figure_2, fig.cap="Figure 2. Spread rates for maize beginning from Guilà Naquitz to sites in the US Southwest on maize’s northern frontier as defined in the text. A) Only change in northing contributes to the distance categories differentiated in the legend. B) As in A, except rates are expressed as Euclidean distances.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

nodes <-
  maize_frontier_mexico %>%
  dplyr::mutate(
    SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName),
    SiteName = ifelse(SiteName == "Guila Naquitz", "Guilá Naquitz", SiteName)
  ) %>%
  dplyr::arrange(Date) %>%
  dplyr::mutate(id = dplyr::row_number()) %>%
  dplyr::mutate(
    lead_row = geometry[dplyr::row_number() + 1],
    dist_deg = as.numeric(dplyr::lag(
      sf::st_distance(geometry, lead_row, by_element = TRUE, which = "Euclidean")
    )),
    dist_km = dist_deg * 110.574,
    delta_deg = dist_deg / (Date - dplyr::lag(Date)),
    delta_km = dist_km / (Date - dplyr::lag(Date)),
    delta_lat_deg = ((Lat - dplyr::lag(Lat))) / (Date - dplyr::lag(Date)),
    delta_lat_deg = as.numeric(sprintf("%.4f", delta_lat_deg)),
    delta_lat_km = delta_lat_deg * 110.574
  ) %>%
  dplyr::select(-lead_row) %>%
  sf::st_drop_geometry()

maize_spread_table <- nodes %>%
  dplyr::arrange(Date) %>%
  dplyr::select(
    LabID,
    SiteName,
    NAME,
    REGION,
    Country,
    Date,
    Age,
    Error,
    MedianBP,
    GDD,
    elevation,
    dist_deg:delta_lat_km,
    dplyr::starts_with("TwoSigma")
  ) %>%
  dplyr::mutate(
    across(c(NAME, REGION), tools::toTitleCase),
    REGION = ifelse(SiteName == "Guilá Naquitz", "Oaxaca", REGION)
  )

sigma_range <- 
  maize_spread_table %>% 
  dplyr::select(LabID, dplyr::starts_with("TwoSigma")) %>% 
  tidyr::pivot_longer(-LabID) %>% 
    dplyr::filter(!is.na(value)) %>% 
    dplyr::group_by(LabID) %>% 
    dplyr::filter(dplyr::row_number() %in% c(1, dplyr::n())) %>% 
    tidyr::separate(value, into = c("early", "late"), sep = " to ") %>% 
    tidyr::pivot_longer(c(early, late), names_to = "date") %>% 
    dplyr::mutate(range = paste0(max(value), " to ", min(value))) %>% 
    dplyr::select(-c(name, date, value)) %>% 
    dplyr::ungroup() %>% 
    dplyr::distinct()


nodes1 <- nodes %>%
  dplyr::select(
    c(
      "id",
      "Long",
      "Lat",
      "SiteName",
      "Date",
      "dist_deg",
      "dist_km",
      "delta_deg",
      "delta_km",
      "delta_lat_deg",
      "delta_lat_km"
    )
  )

names(nodes1) <-
  c(
    "id",
    "lon",
    "lat",
    "name",
    "date",
    "dist_deg",
    "dist_km",
    "delta_deg",
    "delta_km",
    "delta_lat_deg",
    "delta_lat_km"
  )

country_shapes1 <-
  ggplot2::geom_polygon(
    aes(x = long, y = lat, group = group),
    data = map_data('world'),
    fill = "white",
    color = "black",
    size = 0.15
  )

country_shapes2 <-
  ggplot2::geom_polygon(
    aes(x = long, y = lat, group = group),
    data = map_data('state'),
    fill = "white",
    color = "black",
    size = 0.15
  )

mapcoords1 <-
  ggplot2::coord_fixed(xlim = c(-114.814521,-92),
                       ylim = c(14.72716, 42.001167))

legend.breaks = 3


# Figure 2a. Change in N latitude (°) per year.
edges2a <- data.frame(from = 1:13,
                     to = 2:14,
                     color = nodes1$delta_lat_km[!is.na(nodes1$delta_lat_km)])

edges_for_plot2a <- edges2a %>%
  dplyr::inner_join(nodes1 %>% dplyr::select(id, lon, lat), by = c('from' = 'id')) %>%
  dplyr::rename(x = lon, y = lat) %>%
  dplyr::inner_join(nodes1 %>% dplyr::select(id, lon, lat), by = c('to' = 'id')) %>%
  dplyr::rename(xend = lon, yend = lat)

assertthat::assert_that(nrow(edges_for_plot2a) == nrow(edges2a))

maize_frontier_mexico_delta_Lat <- ggplot(nodes1) +
  country_shapes1 +
  country_shapes2 +
  geom_curve(
    aes(
      x = x,
      y = y,
      xend = xend,
      yend = yend,
      color = color
    ),
    arrow = arrow(length = unit(0.02, "npc"), type = "closed"),
    data = edges_for_plot2a,
    curvature = 0.3,
    alpha = 1.0,
    size = 1.0
  ) +
  scale_color_distiller(palette = "RdYlGn",
                        direction = 1,
                        name = "\u0394 km N/year") +
  geom_point(
    aes(x = lon, y = lat),
    # draw nodes
    shape = 21,
    fill = 'white',
    color = 'black',
    stroke = 0.5,
    size = 3.5
  ) +
  ggrepel::geom_label_repel(
    aes(x = lon, y = lat, label = name),
    # draw text labels
    size = 2,
    color = "black",
    min.segment.length = 0.15,
    fill = alpha(c("white"), 0.7),
    force = 34,
    family = "Helvetica",
    max.overlaps = 15
  ) +
  mapcoords1 +
  theme_void() +
  theme(
    plot.margin = grid::unit(c(2, 2, 2, 2), "mm"),
    legend.box.spacing = unit(20.0, "pt"),
    legend.key.size = unit(0.55, "cm"),
    legend.text = element_text(size = 14, family = "Helvetica"),
    legend.title = element_text(
      size = 18,
      face = "bold",
      family = "Helvetica"
    )
  )

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_2a.png"),
#   plot = maize_frontier_mexico_delta_Lat,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )


# Figure 2b. Change in euclidean distance (km) per year.
edges2b <- data.frame(from = 1:13,
                     to = 2:14,
                     color = nodes1$delta_km[!is.na(nodes1$delta_km)])

edges_for_plot2b <- edges2b %>%
  dplyr::inner_join(nodes1 %>% dplyr::select(id, lon, lat), by = c('from' = 'id')) %>%
  dplyr::rename(x = lon, y = lat) %>%
  dplyr::inner_join(nodes1 %>% dplyr::select(id, lon, lat), by = c('to' = 'id')) %>%
  dplyr::rename(xend = lon, yend = lat)

assertthat::assert_that(nrow(edges_for_plot2b) == nrow(edges2b))

maize_frontier_mexico_euclidean_dist <- ggplot(nodes1) +
  country_shapes1 +
  country_shapes2 +
  geom_curve(
    aes(
      x = x,
      y = y,
      xend = xend,
      yend = yend,
      color = color
    ),
    arrow = arrow(length = unit(0.02, "npc"), type = "closed"),
    data = edges_for_plot2b,
    curvature = 0.3,
    alpha = 1.0,
    size = 1.0
  ) +
  scale_color_distiller(palette = "RdYlGn",
                        direction = 1,
                        name = "\u0394 km ED/year") +
  geom_point(
    aes(x = lon, y = lat),
    # draw nodes
    shape = 21,
    fill = 'white',
    color = 'black',
    stroke = 0.5,
    size = 3.5
  ) +
  ggrepel::geom_label_repel(
    aes(x = lon, y = lat, label = name),
    # draw text labels
    size = 2,
    color = "black",
    min.segment.length = 0.15,
    fill = alpha(c("white"), 0.7),
    force = 34,
    family = "Helvetica",
    max.overlaps = 15
  ) +
  mapcoords1 +
  theme_void() +
  theme(
    plot.margin = grid::unit(c(2, 2, 2, 2), "mm"),
    legend.box.spacing = unit(20.0, "pt"),
    legend.key.size = unit(0.55, "cm"),
    legend.text = element_text(size = 14, family = "Helvetica"),
    legend.title = element_text(
      size = 18,
      face = "bold",
      family = "Helvetica"
    )
  )

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_2b.png"),
#   plot = maize_frontier_mexico_euclidean_dist,
#   # width = 15.51,
#   # height = 10.32,
#   dpi = 600
# )

Figure_2 <-
  cowplot::plot_grid(
    maize_frontier_mexico_delta_Lat,
    maize_frontier_mexico_euclidean_dist,
    labels = "AUTO",
    label_x = 0.25,
    label_y = 0.99,
    align = "v",
    nrow = 2,
    ncol = 1
  )

Figure_2

# ggplot2::ggsave(
#   Figure_2,
#   filename = here::here("vignettes/figures/Figure_2.png"),
#   width = 13.7,
#   height = 10.2,
#   dpi = 900
# )

```

```{r Table1, echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'asis'}

table1 <- maize_spread_table %>% 
  dplyr::select(LabID:elevation) %>%
  dplyr::left_join(., sigma_range, by = "LabID") %>% 
  dplyr::select(LabID:Country, Age, Error, range, MedianBP, Date, everything()) %>% 
  dplyr::mutate(References = "")

names(table1) <-
  c(
    "Lab ID",
    "Site Name",
    "County",
    "State",
    "Country",
    "14C Age",
    "14C Error",
    "95% Confidence Interval BP",
    "Median cal BP",
    "Median cal BC/AD",
    "GDD (°C)",
    "Elevation (m)",
    "References"
  )

# readr::write_excel_csv(table1,
#                        here::here("vignettes/tables/Table_1.csv"))

knitr::kable(table1,
             caption = "Table 1. Sites on the Maize Northern Frontier (MNF), beginning with Guilá Naquitz. Here, we include information on the uncalibrated and calibrated dates, as well as other information on elevation and GDD.")

```

```{r Table2, echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'asis'}

table2 <- maize_spread_table %>% 
  dplyr::select(LabID:SiteName, Date, dist_deg:delta_lat_km) %>% 
  dplyr::mutate(dplyr::across(dist_deg:delta_lat_km, ~ sprintf("%.3f", .x)))

names(table2) <-
  c(
    "Lab ID",
    "Site Name",
    "Median cal BC/AD",
    "Euclidean Distance (Degrees)",
    "Euclidean Distance (km)",
    "\u0394 Degrees Euclidean Distance/Year",
    "\u0394 km Euclidean Distance/Year",
    "\u0394 Degrees North/Year",
    "\u0394 km North/Year"
  )

# readr::write_excel_csv(table2,
#                        here::here("vignettes/tables/Table_2.csv"))

knitr::kable(table2,
             caption = "Table 2. Distances and spread rates for sites on the Maize Northern Frontier (MNF), beginning with Guilá Naquitz. Distances and maize spread rates are between a previous row and a current row. For example, Guilá Naquitz is blank because the distance between Guilá Naquitz and Las Capas is provided in the Las Capas row.")
             
```

```{r Figure_3, fig.cap="Figure 3. Spread rates for maize beginning from Las Capas to sites in the US Southwest on the MNF as defined in the text. A) Only change in northing contributes to the distance categories differentiated in the legend. B) As in A, except rates are expressed as Euclidean distances.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

nodes2 <- nodes1 %>%
  dplyr::filter(name != "Guilá Naquitz") %>%
  dplyr::mutate(id = dplyr::row_number())

country_shapes1 <-
  ggplot2::geom_polygon(
    aes(x = long, y = lat, group = group),
    data = map_data('state'),
    fill = "white",
    color = "black",
    size = 0.15
  )

mapcoords1 <-
  ggplot2::coord_fixed(xlim = c(-114.814521, -105),
                       ylim = c(31.333254, 42))

legend.breaks = 3

# Figure 3a. Change in km North latitude per year.
edges3a <- data.frame(from = 1:12,
                      to = 2:13,
                      color = nodes2$delta_lat_km[-1])

edges_for_plot3a <- edges3a %>%
  dplyr::inner_join(nodes2 %>% dplyr::select(id, lon, lat), by = c('from' = 'id')) %>%
  dplyr::rename(x = lon, y = lat) %>%
  dplyr::inner_join(nodes2 %>% dplyr::select(id, lon, lat), by = c('to' = 'id')) %>%
  dplyr::rename(xend = lon, yend = lat)

assertthat::assert_that(nrow(edges_for_plot3a) == nrow(edges3a))

maize_frontier_SW_delta_Lat <- ggplot(nodes2) +
  country_shapes1 +
  geom_curve(
    aes(
      x = x,
      y = y,
      xend = xend,
      yend = yend,
      color = color
    ),
    arrow = arrow(length = unit(0.02, "npc"), type = "closed"),
    data = edges_for_plot3a,
    curvature = 0.3,
    alpha = 1.0,
    size = 1.0
  ) +
  scale_color_distiller(palette = "RdYlGn",
                        direction = 1,
                        name = "\u0394 km N/year") +
  geom_point(
    aes(x = lon, y = lat),
    # draw nodes
    shape = 21,
    fill = 'white',
    color = 'black',
    stroke = 0.5,
    size = 3.5
  ) +
  ggrepel::geom_label_repel(
    aes(x = lon, y = lat, label = name),
    # draw text labels
    size = 2.75,
    color = "black",
    min.segment.length = 0.15,
    fill = alpha(c("white"), 0.7),
    force = 34,
    family = "Helvetica"
  ) +
  mapcoords1 +
  theme_void() +
  theme(
    plot.margin = grid::unit(c(2, 2, 2, 2), "mm"),
    legend.box.spacing = unit(20.0, "pt"),
    legend.key.size = unit(0.6, "cm"),
    legend.text = element_text(size = 14, family = "Helvetica"),
    legend.title = element_text(
      size = 18,
      face = "bold",
      family = "Helvetica"
    )
  )

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/spread_rate.png"),
#   plot = maize_frontier_SW_delta_Lat,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

# Figure 3b. Change in euclidean distance (km) per year.
edges3b <- data.frame(from = 1:12,
                      to = 2:13,
                      color = nodes2$delta_km[-1])

edges_for_plot3b <- edges3b %>%
  dplyr::inner_join(nodes2 %>% dplyr::select(id, lon, lat), by = c('from' = 'id')) %>%
  dplyr::rename(x = lon, y = lat) %>%
  dplyr::inner_join(nodes2 %>% dplyr::select(id, lon, lat), by = c('to' = 'id')) %>%
  dplyr::rename(xend = lon, yend = lat)

assertthat::assert_that(nrow(edges_for_plot3b) == nrow(edges3b))

maize_frontier_SW_delta_euclidean <- ggplot(nodes2) +
  country_shapes1 +
  geom_curve(
    aes(
      x = x,
      y = y,
      xend = xend,
      yend = yend,
      color = color
    ),
    arrow = arrow(length = unit(0.02, "npc"), type = "closed"),
    data = edges_for_plot3b,
    curvature = 0.3,
    alpha = 1.0,
    size = 1.0
  ) +
  scale_color_distiller(palette = "RdYlGn",
                        direction = 1,
                        name = "\u0394 km ED/year") +
  geom_point(
    aes(x = lon, y = lat),
    # draw nodes
    shape = 21,
    fill = 'white',
    color = 'black',
    stroke = 0.5,
    size = 3.5
  ) +
  ggrepel::geom_label_repel(
    aes(x = lon, y = lat, label = name),
    # draw text labels
    size = 2.75,
    color = "black",
    min.segment.length = 0.15,
    fill = alpha(c("white"), 0.7),
    force = 34,
    nudge_x = -0.51,
    nudge_y = 0.15,
    family = "Helvetica"
  ) +
  mapcoords1 +
  theme_void() +
  theme(
    legend.box.spacing = unit(20.0, "pt"),
    legend.key.size = unit(0.6, "cm"),
    legend.text = element_text(size = 14, family = "Helvetica"),
    legend.title = element_text(
      size = 18,
      face = "bold",
      family = "Helvetica"
    )
  )

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/spread_rate_euclidean.png"),
#   plot = maize_frontier_SW_delta_euclidean,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

Figure_3 <-
  cowplot::plot_grid(
    maize_frontier_SW_delta_Lat,
    maize_frontier_SW_delta_euclidean,
    labels = "AUTO",
    label_x = 0.233,
    label_y = 0.99,
    align = "v",
    nrow = 2,
    ncol = 1
  )

Figure_3

# ggplot2::ggsave(
#   Figure_3,
#   filename = here::here("vignettes/figures/Figure_3.png"),
#   width = 13.7,
#   height = 10.2,
#   dpi = 900
# )

```

```{r Figure_4, fig.cap="Figure 4. Location of sites on the maize northern frontier and other sites mentioned in the text relative to A) elevation and B) annual accumulated GDDs (°C).", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

# GDD Map

# Get raster for the North American Southwest and convert to a dataframe.
gdd_cropped <-
  raster::crop(cropSpreadR::PRISM_annual_gdd, cropDiffusionR::swus_states)
gdd_cropped <-
  raster::mask(gdd_cropped, cropDiffusionR::swus_states)

NASW_gdd_data <-
  raster::as.data.frame(gdd_cropped, xy = TRUE)
# Assign names to the columns.
names(NASW_gdd_data) <- c("long", "lat", "GDD")
#Transform the data for plotting in ggplot.
NASW_gdd_data_transformed <-
  usmap::usmap_transform(
    NASW_gdd_data,
    input_names = c("long", "lat"),
    output_names = c("x", "y")
  )

maize_frontier_swus_tp <- maize_frontier_swus %>%
  dplyr::mutate(`Site Type` = "Maize Frontier") %>%
  dplyr::bind_rows(
    .,
    MDB %>%
      dplyr::filter(SiteName == "Turkey Pen Ruin") %>%
      dplyr::slice_min(Date) %>%
      dplyr::mutate(Long1 = Long, Lat1 = Lat) %>%
      # Convert to spatial object
      sf::st_as_sf(coords = c("Long1", "Lat1"),
                   crs = 4326) %>%
      sf::st_make_valid() %>%
      dplyr::mutate(`Site Type` = "Other")
  ) %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName))

GDD_output <- ggplot2::ggplot(data = NASW_gdd_data_transformed) +
  geom_raster(aes(x = long, y = lat, fill = GDD),
              alpha = 0.8,
              na.rm = TRUE) +
  # For GDD raster in greyscale.
  # scale_fill_gradientn(
  #   colors = grey.colors(
  #     10,
  #     start = 0.15,
  #     end = 0.9,
  #     rev = TRUE
  #   ),
  #   na.value = "transparent",
  #   name = "Annual GDDs "
  # ) +
#For GDD raster in color.
scale_fill_gradientn(colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(255),
                     na.value = "transparent",
                     name = "Annual GDD (°C)") +
  coord_equal() +
  theme_bw() +
  geom_sf(
    data = cropDiffusionR::swus_states,
    color = "#2b2b2b",
    fill = "transparent",
    size = 0.4,
    inherit.aes = FALSE
  ) +
  geom_sf(
    data = maize_frontier_swus_tp,
    aes(geometry = geometry,
        shape = `Site Type`),
    inherit.aes = FALSE,
    size = 4,
    show.legend = "point"
  ) +
  scale_shape_manual(values = c(16, 17)) +
  ggrepel::geom_label_repel(
    data = maize_frontier_swus_tp,
    aes(label = SiteName, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0.4,
    force = 47,
    label.size = 0,
    size = 4.5,
    inherit.aes = FALSE,
    max.overlaps = 20,
    fill = alpha(c("white"), 0.7)
  ) +
  theme(
    #legend.position = c(.15, .2),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_blank(),
    axis.text.x = element_blank(),
    axis.title = element_blank(),
    axis.line = element_blank(),
    axis.ticks = element_blank(),
    legend.text = element_text(size = 15, family = "Helvetica"),
    legend.title = element_text(
      size = 16,
      family = "Helvetica",
      face = "bold"
    )
  )


# Elevation Map

NASW_ned_data <-
  raster::as.data.frame(swus_NED, xy = TRUE)
# Assign names to the columns.
names(NASW_ned_data) <- c("long", "lat", "elev")
# Transform the data for plotting in ggplot.
swus_NED_transformed <-
  usmap::usmap_transform(
    NASW_ned_data,
    input_names = c("long", "lat"),
    output_names = c("x", "y")
  )

NED_output <- ggplot2::ggplot(data = swus_NED_transformed) +
  geom_raster(aes(x = long, y = lat, fill = elev),
              alpha = 0.8,
              na.rm = TRUE) +
  # For GDD raster in greyscale.
  # scale_fill_gradientn(
  #   colors = grey.colors(
  #     10,
  #     start = 0.15,
  #     end = 0.9,
  #     rev = TRUE
  #   ),
  #   na.value = "transparent",
  #   name = "Annual GDDs "
  # ) +
#For GDD raster in color.
scale_fill_gradientn(
  colors = grDevices::colorRampPalette(MetBrewer::met.brewer("VanGogh3"))(256),
  na.value = "transparent",
  name = "Elevation (m)"
) +
  coord_equal() +
  theme_bw() +
  geom_sf(
    data = cropDiffusionR::swus_states,
    color = "#2b2b2b",
    fill = "transparent",
    size = 0.4,
    inherit.aes = FALSE
  ) +
  geom_sf(
    data = maize_frontier_swus_tp,
    aes(geometry = geometry,
        shape = `Site Type`),
    inherit.aes = FALSE,
    size = 4,
    show.legend = "point"
  ) +
  scale_shape_manual(values = c(16, 17)) +
  ggrepel::geom_label_repel(
    data = maize_frontier_swus_tp,
    aes(label = SiteName, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0.4,
    force = 47,
    label.size = 0,
    size = 4.5,
    inherit.aes = FALSE,
    max.overlaps = 20,
    fill = alpha(c("white"), 0.7)
  ) +
  theme(
    #legend.position = c(.15, .2),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_blank(),
    axis.text.x = element_blank(),
    axis.title = element_blank(),
    axis.line = element_blank(),
    axis.ticks = element_blank(),
    legend.text = element_text(size = 15, family = "Helvetica"),
    legend.title = element_text(
      size = 16,
      family = "Helvetica",
      face = "bold"
    )
  )

fig_envA <-
  cowplot::plot_grid(
    NED_output,
    GDD_output,
    labels = "AUTO",
    label_x = 0.26,
    align = "v",
    nrow = 2,
    ncol = 1
  )

fig_envA

# ggplot2::ggsave(
#   fig_envA,
#   filename = here::here("vignettes/figures/Figure_4.png"),
#   dpi = 600
# )

```

```{r Figure_5, fig.cap="Figure 5. Box plots showing elevation, GDD, and latitude through time for the earliest maize sample from each of 275 sites in our database (A. Gillreath-Brown and Kohler, n.d.), excluding 19 sites post-dating AD 1500 or south of 30° Latitude. Samples are binned by 500-year periods and plotted on the midpoint of each bin.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

MDB_env2 <- MDB_env %>%
  dplyr::filter(Lat >= 30 & Date < 1500) %>%
  dplyr::group_by(SiteIDName) %>%
  dplyr::slice_min(Date, with_ties = FALSE) %>%
  dplyr::group_by(gr = cut(
    Date,
    breaks = seq(-4000, 1500, 500),
    labels = seq(-3750, 1250, 500)
  )) %>%
  dplyr::mutate(date2 = as.numeric(as.character(gr))) %>%
  dplyr::ungroup() %>%
  dplyr::select(-gr) %>%
  dplyr::arrange(date2) %>%
  dplyr::rename(Latitude = Lat, Elevation = elevation)

env_summaries <- MDB_env2 %>%
  dplyr::summarise(dplyr::across(c(Latitude, GDD, Elevation), ~ mean(.x, na.rm = TRUE)))

env_summaries_bp <-
  MDB_env2 %>%
  dplyr::select(date2, Latitude, GDD, Elevation) %>%
  dplyr::rename(
    `Latitude (°)` = Latitude,
    `Annual GDD (°C)` = GDD,
    `Elevation (m)` = Elevation
  ) %>%
  tidyr::pivot_longer(-date2) %>%
  ggplot(aes(
    x = factor(date2),
    y = value,
    fill = factor(name)
  )) +
  geom_boxplot() +
  facet_grid(vars(name), scales = "free") +
  labs(fill = "") +
  xlab("Years BC/AD") +
  ylab("Value") +
  theme_bw() +
  theme(
    strip.text.y = element_text(size = 12),
    legend.key.size = unit(1.5, "cm"),
    legend.text = element_text(size = 12),
    axis.ticks.length = unit(0.125, "inch"),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
    axis.text = element_text(
      size = 26,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_blank(),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 5,
        r = 10,
        b = 10,
        l = 10
      )
    ),
    axis.line = element_line(colour = "black"),
    plot.margin = margin(
      t = 5,
      r = 25,
      b = 5,
      l = 5
    )
  )

env_summaries_bp

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_5.png"),
#   plot = env_summaries_bp,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

```

```{r multiple_linear_regression, echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'asis'}

mlr <- MDB_env %>%
  dplyr::group_by(SiteIDName) %>%
  dplyr::slice_min(Date) %>%
  dplyr::filter(Lat >= 30 & Date < 1500) %>%
  dplyr::mutate(bin = cut(Date,
                          breaks = seq(-3700, 1500, 100))) %>%
  dplyr::group_by(bin) %>%
  dplyr::mutate(group = dplyr::cur_group_id()) %>%
  dplyr::filter(LabID %in% maize_frontier_swus$LabID) %>%
  dplyr::ungroup()

mlr2 <- mlr %>%
  dplyr::select(Date, Lat, elevation, GDD) %>% 
  dplyr::mutate(dplyr::across(everything(), scales::rescale))

gdd_lat <- lm(GDD ~ Lat, data = mlr2)
gdd_elev <- lm(GDD ~ elevation, data = mlr2)
date_lat <- lm(Date ~ Lat, data = mlr2)
date_elev <- lm(Date ~ elevation, data = mlr2)
date_gdd <- lm(Date ~ GDD, data = mlr2)

stats_table <-
  purrr::map_dfr(list(gdd_lat, gdd_elev, date_lat, date_elev, date_gdd), function(x) {
    x %>%
      summary() %>%
      {
        data.frame(
          formula = as.character(.$call)[[2]],
          slope = formatC(.$coefficients[2, 1], digits = 3, format = "f"),
          r2 = formatC(.$r.squared, digits = 2, format = "f"),
          pvalue = formatC(.$coefficients[2, 4],
                           digits = 3,
                           format = "f"),
          Fstatistic = formatC(.$fstatistic[[1]],
                               digits = 2,
                               format = "f")
        )
      } })

knitr::kable(stats_table,
             caption = "Statistics Table. Summary of statistics for the linear models shown in Figure 6.")

```

```{r Figure_6, fig.cap="Figure 6. A directed acyclic graph of the relationships between the ages of maize on the MFN and the latitude, elevation, and GDD for the sites with this maize. Beta weights on each arrow indicate the size of the effect of the independent variable on the dependent. Solid lines indicate a significant effect (at p < 0.01); the dashed line indicates an effect that is not significant.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

dag_coords <-
  tibble::tibble(
    name = c("Latitude", "Elevation", "GDD", "Date"),
    x    = c(1, 2, 1.5, 1.5),
    y    = c(3, 3, 2, 1)
  )

dag_out <- ggdag::dagify(GDD ~ Latitude + Elevation,
                         Date ~ GDD,
                         Date ~ Latitude + Elevation,
                         coords = dag_coords) %>%
  ggdag::tidy_dagitty() %>%
  cropSpreadR::shorten_dag_arrows(proportion = 0.27)

dag_out$data <-
  dplyr::left_join(
    x = dag_out$data,
    stats_table %>%
      tidyr::separate(formula, into = c("to", "name"), sep = " ~ ") %>%
      dplyr::mutate(across(
        c(to, name),
        ~ dplyr::case_when(.x == "Lat" ~ "Latitude",
                           .x == "elevation" ~ "Elevation",
                           TRUE ~ .x)
      )) %>%
      dplyr::mutate(sig = pvalue <= 0.05) %>%
      dplyr::select(name, to, slope, sig),
    by = c("name", "to")
  ) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(
    slope = sprintf("%.2f", as.numeric(slope)),
    x_label = ((xstart - xend) / 2) + xend,
    y_label = ((ystart - yend) / 2) + yend,
    edge_linetype = ifelse(sig == TRUE, "solid", "dashed")
  ) %>%
  dplyr::ungroup()

dag_figure <- dag_out %>%
  ggplot(aes(
    x = x,
    y = y,
    xend = xend,
    yend = yend
  )) +
  ggdag::geom_dag_point(color = "cornflowerblue",
                        alpha = 1 / 4,
                        size = 45) +
  ggdag::geom_dag_text(color = "cornflowerblue",
                       size = 7,
                       fontface = "bold") +
  ggdag::geom_dag_edges(
    aes(
      x = xstart,
      y = ystart,
      edge_linetype = edge_linetype
    ),
    edge_color = "cornflowerblue",
    arrow_directed = grid::arrow(length = grid::unit(12, "pt"), type = "closed"),
    edge_width = 1.35
  ) +
  ggdag::geom_dag_label_repel(
    aes(x = x_label, y = y_label, label = slope),
    segment.size = 0,
    force = 0,
    size = 6,
    label.size = NA
  ) +
  scale_x_continuous(NULL, breaks = NULL, expand = c(0.1, 0.1)) +
  scale_y_continuous(NULL, breaks = NULL, expand = c(0.2, 0.2)) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        panel.border = element_blank())

dag_figure

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_6.png"),
#   plot = dag_figure,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

```

```{r Figure_7, fig.cap="Figure 7. ", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

nodes7 <-
  maize_frontier_swus %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName)) %>%
  dplyr::arrange(Date) %>%
  dplyr::mutate(id = dplyr::row_number()) %>%
  dplyr::select(id, Long, Lat, SiteName, Date) %>%
  sf::st_drop_geometry() %>%
  dplyr::mutate(delta = ((Lat - dplyr::lag(Lat)) * 110.574) / (Date - dplyr::lag(Date)),
                delta = as.numeric(sprintf("%.2f", delta)))

anom_point <- temp_anom %>%
  dplyr::select(year, anom) %>%
  dplyr::mutate(
    anom_25 = dplyr::lag(anom, 25),
    anom_50 = dplyr::lag(anom, 50),
    anom_75 = dplyr::lag(anom, 75),
    anom_100 = dplyr::lag(anom, 100),
    anom_200 = dplyr::lag(anom, 200)
  )

lags_point <- nodes7 %>%
  dplyr::filter(SiteName != "Las Capas") %>%
  dplyr::left_join(., anom_point, by = c("Date" = "year"))

# Get fit, and make a variable with labels for points outside conf. interval
fit <- lm(delta ~ anom_25, data = lags_point)
dat <- predict(fit, interval = "confidence")

lags_point$inside <-
  ifelse(
    lags_point$delta < dat[, "upr"] &
      lags_point$delta > dat[, "lwr"],
    "",
    as.character(lags_point$SiteName)
  )

temp_spreadrate <- lags_point %>%
  ggplot2::ggplot(aes(x = anom_25, y = delta)) +
  geom_point(size = 4) +
  stat_smooth(method = "lm", size = 2) +
  scale_x_continuous(breaks = seq(-1.5, 0, 0.25)) +
  scale_y_continuous(breaks = seq(-0.4, 1.0, 0.2)) +
  ggrepel::geom_text_repel(
    aes(label = inside),
    cex = 6.5,
    show.legend  = F,
    min.segment.length = unit(0.1, 'lines'),
    box.padding = 1
  ) +
  theme_bw() +
  xlab("Temperature Anomaly (°C)") +
  ylab("\u0394 km N/year") +
  theme(
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_text(
      size = 26,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 5,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 5,
        r = 10,
        b = 10,
        l = 10
      )
    ),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"),
    plot.margin = margin(
      t = 5,
      r = 25,
      b = 5,
      l = 5
    )
  )

temp_spreadrate

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_7.png"),
#   plot = temp_spreadrate,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

```

```{r Figure_8, fig.cap="Figure 8. The 12 counties in the four-state US Southwest with a site on maize’s northern frontier (2 sites fall within Daggett County, Utah). The northern frontier of maize expansion was located almost entirely within the western Southwest.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

sw_rc <- cropSpreadR::sw_rc

#frontier_counties <- unique(maize_frontier_swus$county_st)[!grepl("delta, colorado", unique(maize_frontier_swus$county_st))]
frontier_counties <-
  c(
    "sevier, utah",
    "delta, colorado",
    "daggett, utah",
    "catron, new mexico",
    "navajo, arizona",
    "salt lake, utah",
    "pima, arizona",
    "kane, utah",
    "montrose, colorado",
    "uintah, utah",
    "box elder, utah",
    "san juan, utah"
  )

sw_rc_subset <- sw_rc %>%
  dplyr::filter(county_st %in% frontier_counties) %>% 
  dplyr::group_by(county_st) %>% 
  dplyr::summarise(no_samples = dplyr::n(),
                   no_sites = length(unique(SiteIDName)))

# Get sites that are on the northern maize frontier through time from cropDiffusionR package.
temp <- tempfile()
download.file(
  "https://raw.githubusercontent.com/Archaeo-Programmer/cropDiffusionR/master/vignettes/output/maize_frontier_swus.rds",
  temp
)
maize_frontier_swus <-
  readr::read_rds(temp, "maize_frontier_swus.rds") %>%
  dplyr::left_join(., MDB %>% dplyr::select(LabID, Long, Lat), by = "LabID")
unlink(temp)

# Create new column for fill, which will be used for plotting. Grey if the county occurs on the maize northern frontier, and blank if it does not.
SWUS_counties_fill <- cropSpreadR::SWUS_counties %>%
  dplyr::left_join(., sw_rc_subset, by = "county_st") %>% 
  dplyr::mutate(
    #fill = ifelse(county_st %in% maize_frontier_swus$county_st, "grey", NA),
    label = ifelse(
      county_st %in% maize_frontier_swus$county_st,
      tools::toTitleCase(NAME),
      NA
    )
  )

# Plot.
spd_counties_samples <- ggplot(data = SWUS_counties_fill) +
  geom_sf(color = "black",
          aes(fill = no_samples),
          size = 0.25,
          na.rm = TRUE
          ) +
  scale_fill_continuous(low="#FFC20A", high="#0C7BDC",
                       guide="colorbar",na.value="transparent",
                       name = "no. radiocarbon dates") +
  ggrepel::geom_label_repel(
    aes(label = label, geometry = geom),
    stat = "sf_coordinates",
    min.segment.length = 0.1,
    size = 2.5
  ) +
coord_sf() +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.key.width=unit(0.05,"npc"),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_blank(),
    axis.text.x = element_blank(),
    axis.title = element_blank(),
    axis.line = element_blank(),
    axis.ticks = element_blank(),
    legend.text = element_text(size = 12, family = "Helvetica"),
    legend.title = element_text(
      size = 14,
      family = "Helvetica",
      face = "bold"
    )
  )

spd_counties_sites <- ggplot(data = SWUS_counties_fill) +
  geom_sf(color = "black",
          aes(fill = no_sites),
          size = 0.25,
          na.rm = TRUE
          ) +
  scale_fill_continuous(low="#FFC20A", high="#0C7BDC",
                       guide="colorbar",na.value="transparent",
                       name = "no. sites with \nradiocarbon dates") +
  ggrepel::geom_label_repel(
    aes(label = label, geometry = geom),
    stat = "sf_coordinates",
    min.segment.length = 0.1,
    size = 2.5
  ) +
coord_sf() +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.key.width=unit(0.05,"npc"),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_blank(),
    axis.text.x = element_blank(),
    axis.title = element_blank(),
    axis.line = element_blank(),
    axis.ticks = element_blank(),
    legend.text = element_text(size = 12, family = "Helvetica"),
    legend.title = element_text(
      size = 14,
      family = "Helvetica",
      face = "bold"
    )
  )

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_8.png"),
#   plot = spd_counties,
#   width = 15.51,
#   height = 10.32,
#   dpi = 600
# )

Figure_8 <-
  cowplot::plot_grid(
    spd_counties_samples,
    spd_counties_sites,
    labels = "AUTO",
    label_x = 0.36,
    label_y = 0.99,
    align = "v",
    nrow = 2,
    ncol = 1
  )

Figure_8

# ggplot2::ggsave(
#   Figure_8,
#   filename = here::here("vignettes/figures/Figure_8.png"),
#   width = 13.7,
#   height = 10.2,
#   dpi = 900
# )

```

```{r Figure_9, fig.cap="Figure 9. ", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.align = "left"}

sw_rc <- cropSpreadR::sw_rc

#frontier_counties <- unique(maize_frontier_swus$county_st)[!grepl("delta, colorado", unique(maize_frontier_swus$county_st))]
frontier_counties <-
  c(
    "sevier, utah",
    "daggett, utah",
    "catron, new mexico",
    "navajo, arizona",
    "salt lake, utah",
    "pima, arizona",
    "kane, utah",
    "montrose, colorado",
    "uintah, utah",
    "box elder, utah",
    "san juan, utah"
  )

sw_rc_subset <- sw_rc %>%
  dplyr::filter(county_st %in% frontier_counties)

nsim = 100

sw.caldates <- sw_rc_subset %>%
  dplyr::select(county_st, Age, Error, calib) %>%
  dplyr::group_by(county_st) %>%
  tidyr::nest() %>%
  dplyr::mutate(monte_carlo = purrr::map(data, function(x) {
    sw.caldates <-
      rcarbon::calibrate(x = x$Age,
                         errors = x$Error,
                         calCurves = 'intcal20')
    rcarbon::modelTest(
      sw.caldates,
      errors = x$Error,
      nsim = nsim,
      timeRange = c(7000, 0),
      model = "exponential",
      runm = 100,
      method = c("calsample")
    )
  }))

monte_carlo_sim <- sw.caldates %>%
  dplyr::mutate(
    pvalue = purrr::map_dbl(monte_carlo, ~ .x$pval),
    n = purrr::map_dbl(monte_carlo, ~ .x$n),
    pvalue = ifelse(pvalue < 0.01, paste0("p < 0.01"), paste0("p = ", sprintf("%.2f", pvalue))),
    monte_carlo = purrr::map(
      monte_carlo,
      ~ .x$result %>% dplyr::mutate(Date = rcarbon::BPtoBCAD(calBP))
    )
  ) %>%
  dplyr::select(-data) %>%
  tidyr::unnest(monte_carlo) %>%
  dplyr::mutate(
    boom_bust = dplyr::case_when(
      PrDens < hi & PrDens > lo ~ "medium",
      PrDens > hi ~ "indianred",
      PrDens < lo ~ "royalblue",
      TRUE ~ "medium"
    )
  )

maize_oldest <- maize_frontier_swus %>%
  sf::st_drop_geometry() %>%
  dplyr::mutate(SiteNameID = dplyr::coalesce(SiteName, SiteID)) %>%
  dplyr::select(county_st, SiteNameID, Date) %>%
  dplyr::mutate(date_maize = Date) %>%
  dplyr::group_by(county_st) %>%
  dplyr::slice_min(Date) %>%
  dplyr::ungroup() %>%
  dplyr::arrange(Date) %>%
  dplyr::filter(county_st %in% frontier_counties)

ckde_test3 <- monte_carlo_sim %>%
  dplyr::left_join(., maize_oldest, by = c("county_st", "Date")) %>%
  dplyr::mutate(
    county_st2 = county_st,
    county_st = tools::toTitleCase(county_st),
    county_st = paste0(county_st, " (n = ", n, "; ",  pvalue, ")")
  ) %>%
  dplyr::mutate(SiteNameID = ifelse(SiteNameID == "Elsinore Burial", "42SV2111", SiteNameID))

county_order <- cropSpreadR::SWUS_counties %>%
  sf::st_drop_geometry() %>%
  dplyr::select(county_st2 = county_st, centroid_lat) %>%
  dplyr::arrange(desc(centroid_lat)) %>%
  dplyr::filter(county_st2 %in% frontier_counties) %>%
  dplyr::left_join(
    .,
    ckde_test3 %>% dplyr::ungroup() %>% dplyr::select(county_st, county_st2) %>% dplyr::distinct(),
    by = "county_st2"
  ) %>%
  dplyr::pull(county_st)

condensed <- monte_carlo_sim %>%
  dplyr::group_by(county_st) %>%
  dplyr::mutate(change = cumsum(boom_bust != lag(boom_bust, def = first(boom_bust)))) %>%
  dplyr::group_by(change, .add = TRUE) %>%
  dplyr::filter(dplyr::row_number() %in% c(1, dplyr::n()) &
                  boom_bust != "medium") %>%
  dplyr::mutate(time = dplyr::row_number()) %>% 
  dplyr::ungroup() %>% 
  tidyr::complete(county_st, change, time) %>%
  tidyr::fill(dplyr::everything()) %>% 
  dplyr::select(-time) %>% 
  dplyr::group_by(grp = stringr::str_c('Column', rep(1:2, length.out = dplyr::n()))) %>%
  dplyr::mutate(rn = dplyr::row_number()) %>%
  dplyr::ungroup() %>%
  dplyr::select(county_st, Date, boom_bust, grp, rn) %>%
  tidyr::pivot_wider(names_from = grp, values_from = Date) %>%
  dplyr::select(-rn,
                xmin = Column1,
                xmax = Column2,
                fill = boom_bust) %>%
  dplyr::mutate(ymin = -Inf, ymax = Inf) %>%
  dplyr::left_join(
    .,
    ckde_test3 %>% dplyr::ungroup() %>% dplyr::select(county_st, county_st2) %>% dplyr::distinct(),
    by = c("county_st" = "county_st2")
  ) %>%
  dplyr::select(-county_st, county_st = county_st.y) %>%
  dplyr::mutate(county_st = factor(county_st, levels = county_order))

monte_carlo_plot <-
  ggplot(data = ckde_test3, aes(y = PrDens, x = Date)) +
  geom_ribbon(aes(ymin = lo, ymax = hi),
              alpha = 0.4,
              fill = "grey70") +
  geom_rect(
    aes(
      xmin = xmin,
      xmax = xmax,
      ymin = ymin,
      ymax = ymax,
      fill = I(fill)
    ),
    alpha = 0.4,
    inherit.aes = FALSE,
    data = condensed
  ) +
  geom_line() +
  geom_point(
    aes(y = PrDens, x = date_maize),
    size = 2.75,
    color = "red",
    inherit.aes = FALSE
  ) +
  ggrepel::geom_label_repel(
    aes(y = PrDens, x = date_maize, label = SiteNameID),
    color = "red",
    cex = 5,
    show.legend  = F,
    min.segment.length = unit(0.1, 'lines'),
    box.padding = 1,
    nudge_x = -550,
    nudge_y = 0.35,
    fill = alpha(c("white"), 0.6),
    label.size = NA
  ) +
  theme_bw() +
  scale_x_continuous(breaks = seq(-4000, 2000, 500),
                     limits = c(-4100, 2000)) +
  scale_y_continuous(
    breaks = function(x)
      seq(0, max(x), length.out = 3),
    labels = function(y)
      sprintf("%.2f", y)
  ) +
  facet_wrap(~ factor(county_st, levels = county_order),
             scales = "free_y",
             ncol = 1) +
  xlab("Years BC/AD") +
  ylab("Summed Probability") +
  theme(
    legend.position = "none",
    panel.spacing = unit(1, "lines"),
    strip.text.x = element_text(
      size = 16,
      family = "Helvetica",
      face = "bold"
    ),
    axis.ticks.length = unit(0.125, "inch"),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    axis.text = element_text(
      size = 20,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 5,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 5,
        r = 10,
        b = 10,
        l = 10
      )
    ),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"),
    plot.margin = margin(
      t = 5,
      r = 25,
      b = 5,
      l = 5
    )
  )

monte_carlo_plot

# ggplot2::ggsave(
#   filename = here::here("vignettes/figures/Figure_9.png"),
#   plot = monte_carlo_plot,
#   width = 16,
#   height = 18,
#   dpi = 600
# )

```

```{r estimating_maize_departure, echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'asis'}

maize_spread_rates <- readr::read_csv(here::here("vignettes/tables/Table_2.csv"), show_col_types = FALSE) %>% 
  dplyr::select(`Lab ID`, `Site Name`, delta_kmN = `Δ km North/Year`)

max_SW_spread <- maize_spread_rates %>% 
  dplyr::filter(!`Site Name` %in% c("Guilá Naquitz", "Las Capas")) %>% 
  dplyr::slice_max(delta_kmN)

max_spread <- max_SW_spread %>% 
  dplyr::slice(2) %>% 
  dplyr::pull(delta_kmN)

LC_ST_rate <- maize_frontier_swus %>% 
  dplyr::filter(SiteName %in% c("Las Capas", "Six Toe Shelter")) %>% 
  dplyr::arrange(Date) %>% 
  dplyr::mutate(delta_lat_deg = ((Lat - dplyr::lag(Lat))) / (Date - dplyr::lag(Date)), 
                delta_lat_deg = as.numeric(sprintf("%.4f", delta_lat_deg)), 
                delta_lat_km = delta_lat_deg * 110.574) %>% 
  dplyr::filter(!is.na(delta_lat_km)) %>% 
  dplyr::pull(delta_lat_km)

LC_WM_rate <- maize_frontier_swus %>% 
  dplyr::filter(SiteName %in% c("Las Capas", "Willard Mounds")) %>% 
  dplyr::arrange(Date) %>% 
  dplyr::mutate(delta_lat_deg = ((Lat - dplyr::lag(Lat))) / (Date - dplyr::lag(Date)), 
                delta_lat_deg = as.numeric(sprintf("%.4f", delta_lat_deg)), 
                delta_lat_km = delta_lat_deg * 110.574) %>% 
  dplyr::filter(!is.na(delta_lat_km)) %>% 
  dplyr::pull(delta_lat_km)

balsas_dissolve <- cropSpreadR::Balsas_Basin %>% 
  dplyr::summarize()

balsas_centroid <- balsas_dissolve %>% 
  sf::st_centroid() %>% 
  sf::st_coordinates()

dist <- (maize_frontier_swus %>% 
  dplyr::filter(SiteName == "Las Capas") %>% 
  dplyr::pull(Lat) - balsas_centroid[1,2]) * 110.574
  
out <- data.frame(type = c("Max SWUS from MNF", "Las Capas to Six Toe Shelter", "Las Capas to Willard Mounds"),
                  rate = c(max_spread, LC_ST_rate, LC_WM_rate)) %>% 
  dplyr::mutate(time = dist / rate)

names(out) <- c("Spread Rate", "Rate (km N/Year)", "Time Elapsed (Years)")

knitr::kable(out,
             caption = "Table S2. Summary for the spread rates of maize. Here, we show the maximum rate between any two sites on the MNF, then we show the rate between Las Capas and two other sites, Six Toe Shelter and Willard Mounds. Time elapsed is the amount of time that it would take to reach the Balsas River Basin at the given spread rate. We use the distance between Las Capas to the centroid of the Balses River Basin to calculate the time elapsed.")

```