Maize_Analysis.Rmd

---
title: "cropDiffusionR: How Maize Farmers in the US Southwest Grew and Prospered under El Niño but Suffered under La Niña"
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'}
knitr::opts_chunk$set(echo = TRUE)

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

```

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

# Load the temperature anomaly reconstruction from paleoMAT.
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)

# Get county shapefiles for SWUS to add as an additional column to MDB.
SWUS_counties <-
  sf::st_as_sf(
    maps::map(
      "county",
      fill = TRUE,
      plot = FALSE,
      regions = "colorado|utah|arizona|new mexico"
    )
  ) %>%
  sf::st_transform(crs = 4326) %>%
  tidyr::separate(ID, into = c("REGION", "NAME"), sep = ",") %>%
  sf::st_make_valid() %>%
  dplyr::mutate(county_st = paste(NAME, REGION, sep = ", "))

# Create the Supplementary Table 1, which does not have site locations.
# write.csv(cropDiffusionR::maizeDB,
#           here::here("vignettes/tables/tableS1.csv"),
#           row.names = FALSE)

# Load location data for each sample. Site location 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(
                         "maize_database/Maize_Database_Locations_DONT_DISTRIBUTE.csv"
                       ),
                       "cropDiffusionR/"
                     ),
                     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(SWUS_counties) %>%
  dplyr::mutate(county_st = ifelse(!is.na(REGION) &
                                     !is.na(NAME), paste(NAME, REGION, sep = ", "), NA)) %>%
  sf::st_drop_geometry()

# If you do not have your own location information, then you can run the code below using the coarse scale location data provided by the p3k14c package. If so, then you can uncomment this code and run instead of lines 52-60 above. However, some of the samples may not be available in the p3k14c package, and locations could vary by at least 10 km from the actual site. p3k14c_data has a column that details the location accuracy. Please refer to Table 2 in https://doi.org/10.1038/s41597-022-01118-7 for more information on the location accuracy codes.
# MDB <- cropDiffusionR::maizeDB %>%
#   dplyr::left_join(
#     .,
#     p3k14c::p3k14c_data %>% dplyr::select(LabID, Long, Lat),
#     by = "LabID"
#   ) %>%
#   # Remove any rows that do not have a location.
#   dplyr::filter(!is.na(Lat) & !is.na(Long)) %>%
#   # Currently, there is an error with coordinates for AA-3308, so correcting them here.
#   dplyr::mutate(Long = dplyr::case_when(SiteName == "Coxcatlan" ~ -97,
#                                  TRUE ~ Long),
#                 Lat = dplyr::case_when(SiteName == "Coxcatlan" ~ 18,
#                                  TRUE ~ Lat))

# Get modern annual accumulated growing degree days at each sample location.
MDB_gdd <- MDB %>%
  sf::st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  dplyr::select(LabID) %>%
  cropDiffusionR::extract_GDD(gdd_raster = cropDiffusionR::WorldClim_annual_gdd,
                              sites = .)

```

# Figures

```{r Figure_1, fig.cap="Figure 1. Low-frequency July temperature anomaly reconstruction for the southwestern United States as defined in inset (based on Gillreath-Brown et al. 2024:Figure 5). The extent of the SWUS inset is approximately: xmin = -112.82, xmax = -105.52, ymin = 32.66, ymax = 38.90. Blue line fitted by loess (local polynomial regression fitting) smoothing (α = 0.15). Negative values on the y-axis indicate temperatures cooler than the 1961–1990 average; the horizontal line at 0 marks the average July temperature from 1961¬–1990. Negative values on x-axis indicate years BC.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 15.51, fig.height = 10.32, fig.align = "left"}

# Create a bounding box for limits of the Gillreath-Brown et al. 2024 study.
bb <-
  c(
    "xmin" = -112.81446,
    "xmax" = -105.52416,
    "ymin" = 32.65757,
    "ymax" = 38.089
  ) %>%
  sf::st_bbox() %>%
  sf::st_as_sfc() %>%
  sf::st_as_sf(crs = 4326) %>%
  sf::st_transform(crs = 4326)

# Get state boundaries for the four corner states.
usa <-
  sf::st_as_sf(maps::map("state", fill = TRUE, plot = FALSE)) %>%
  sf::st_transform(., crs = 4326) %>%
  dplyr::filter(ID %in% c("colorado", "utah", "arizona", "new mexico"))

# Plot study area for the low-frequency temperature reconstruction.
paleomat_extent <-
  ggplot2::ggplot() +
  theme_bw() +
  scale_x_continuous(breaks = seq(-115, -103, 2.0),
                     limits = c(-115, -103.2)) +
  scale_y_continuous(breaks = seq(31, 39, 2.0),
                     limits = c(31.25, 39)) +
  geom_sf(
    data = usa,
    color = "#2b2b2b",
    fill = "transparent",
    size = 0.25
  ) +
  geom_sf(
    data = bb,
    linewidth = 1,
    fill = NA,
    aes(colour = "slategrey",
        linetype = "dashed"),
    show.legend = "line",
    inherit.aes = FALSE
  ) +
  scale_color_manual(
    values = c("slategrey" = "slategrey"),
    labels = c("Gillreath-Brown et al. 2024 \n Study Area"),
    name = "Data Extent"
  ) +
  scale_linetype_manual(
    values = c("dashed" = "dashed"),
    labels = c("Gillreath-Brown et al. 2024 \n Study Area"),
    name = "Data Extent"
  ) +
  labs(colour = "") +
  annotate(
    "text",
    x = -109.2,
    y = 37.5,
    label = "UTAH",
    size = 3.5,
    fontface = 2,
    hjust = 1.0
  ) +
  annotate(
    "text",
    x = -108.90,
    y = 37.5,
    label = "COLORADO",
    size = 3.5,
    fontface = 2,
    hjust = 0
  ) +
  annotate(
    "text",
    x = -109.2,
    y = 36.3,
    label = "ARIZONA",
    size = 3.5,
    fontface = 2,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -108.90,
    y = 36.3,
    label = "NEW MEXICO",
    size = 3.5,
    fontface = 2,
    hjust = 0
  ) +
  theme(
    legend.position = "right",
    legend.title = element_blank(),
    legend.text = element_text(size = 12),
    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(),
    panel.background = element_rect(fill = "transparent",
                                    colour = NA_character_),
    plot.background = element_rect(fill = "transparent",
                                   colour = NA_character_),
    # Necessary to avoid drawing plot outline.
    legend.background = element_rect(fill = "transparent", color = NA),
    legend.box.background = element_rect(fill = "transparent", color = NA),
    legend.key = element_rect(fill = "transparent", color = NA),
    legend.box.margin = margin(-10,-10,-10,-20)
  )

# Plot low-frequency temperature anomaly reconstruction.
anom_plot <- temp_anom %>%
  ggplot2::ggplot(aes(x = year, y = anom)) +
  geom_ribbon(aes(ymax = anom + ci, ymin = anom - ci), alpha = 0.2) +
  geom_line(size = 2, colour = "cornflowerblue") +
  geom_hline(
    yintercept = 0,
    linetype = 'dashed',
    color = c('darkslategray')
  ) +
  scale_x_continuous(
    limits = c(-3000, 2100),
    breaks = seq(-3000, 2000, 500),
    minor_breaks = seq(-2900, 2000, by = 100),
    guide = "axis_minor",
    labels = c(seq(-3000, -500, 500), 1, seq(500, 2000, 500))
  ) +
  scale_y_continuous(
    limits = c(-2.0, 3.0),
    breaks = seq(-2.0, 3.0, 0.5),
    minor_breaks = seq(-2.0, 3.0, by = 0.25),
    guide = "axis_minor"
  ) +
  theme_bw() +
  xlab("Years BC/AD") +
  ylab("Temperature Anomaly (°C)") +
  theme(
    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 = 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
    )
  )

# Plot the two plots together.
plot.with.inset <-
  cowplot::ggdraw() +
  cowplot::draw_plot(anom_plot) +
  cowplot::draw_plot(
    paleomat_extent,
    x = 0.067,
    y = 0.68,
    width = .50,
    height = .30
  )

plot.with.inset

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

```

```{r Figure_2, fig.cap="Figure 2. Low-frequency July temperature anomaly for the SWUS (Gillreath-Brown et al. 2024) and the Northern Hemisphere low-frequency component (Moberg et al. 2005) from AD 133 to 1900. Both series have been standardized to a mean of 0 and an s of 1 for the period in the graph.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 17.25, fig.height = 11.5, fig.align = "left"}

temp_anom_standard <-
  readr::read_csv(
    "https://raw.githubusercontent.com/Archaeo-Programmer/paleoMAT/master/vignettes/tables/supp_table_s2.csv",
    show_col_types = FALSE
  ) %>%
  dplyr::mutate(Year = as.integer(Year)) %>%
  dplyr::filter(Year >= 133 & Year <= 1900) %>%
  dplyr::mutate(Series = "Gillreath-Brown et al. 2024") %>%
  dplyr::select(year = Year, anom = 2, Series)

Moberg_standard <- paleomat::moberg_2005 %>%
  dplyr::mutate(Series = "Moberg et al. 2005") %>%
  dplyr::select(year = date, anom = LF, Series) %>%
  dplyr::filter(year >= 133 & year <= 1900)

Moberg_paleomat <-
  dplyr::bind_rows(temp_anom_standard, Moberg_standard)

# Moberg LF and paleomat
Moberg_paleomat_fig <- Moberg_paleomat %>%
  dplyr::group_by(Series) %>%
  dplyr::mutate(update = (anom - mean(anom)) / sd(anom)) %>%
  ggplot2::ggplot(aes(
    x = year,
    y = update,
    colour = Series,
    linetype = Series
  )) +
  geom_hline(
    yintercept = 0,
    color = "darkgrey",
    lty = 3,
    lwd = 1.25
  ) +
  geom_line(size = 2.0) +
  scale_color_manual(
    values = c("#A60F2D", "cornflowerblue"),
    labels = c("Gillreath-Brown et al. 2024", "Moberg et al. 2005")
  ) +
  scale_linetype_manual(
    values = c("solid", "twodash"),
    labels = c("Gillreath-Brown et al. 2024", "Moberg et al. 2005")
  ) +
  xlab("Years AD") +
  ylab("Standardized Temperature Values (Z-Scores)") +
  scale_x_continuous(
    limits = c(0, 2000),
    breaks = seq(0, 2000, 200),
    minor_breaks = seq(0, 2000, by = 100),
    guide = "axis_minor",
    labels = c(1, seq(200, 2000, 200))
  ) +
  scale_y_continuous(breaks = seq(-2, 3, 0.5)) +
  theme_bw() +
  guides(color = guide_legend(override.aes = list(fill = NA), direction = "horizontal")) +
  theme(
    legend.key.width = unit(1.25, 'cm'),
    legend.text = element_text(size = 19, family = "Helvetica"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.margin = margin(t = -14),
    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 = 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
    )
  )

Moberg_paleomat_fig

# ggplot2::ggsave(
#   Moberg_paleomat_fig,
#   filename = here::here("vignettes/figures/Figure_2.png"),
#   width = 17.25,
#   height = 11.5,
#   dpi = 600
# )

```

```{r Figure_3, fig.cap="Figure 3. The ENSO index is calculated as the difference between the temperature reconstruction by Gillreath-Brown et al. (2024) and the low-frequency component of the Moberg et al. (2005) temperature reconstruction for the Northern Hemisphere from AD 133 to 1900. This index is negative when the Gillreath-Brown et al. reconstruction is cooler than that of Moberg et al. (2005), suggesting El Nino-like conditions prevailed in the SWUS. This index is positive when the Gillreath-Brown et al. (2024) reconstruction is warmer than that of Moberg et al. (2005), suggesting La Niña-like conditions. The ENSO Index is computed from the z-scores for the Gillreath-Brown et al. (2024) and the Moberg et al. (2005) reconstructions as displayed in Figure 2.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 17.25, fig.height = 11.5, fig.align = "left"}

# In this vignette, the temperature anomaly is drawn from the 100 year intervals.
temp_anom_standard2 <-
  temp_anom_standard %>%
  dplyr::mutate(update = (anom - mean(anom)) / sd(anom)) %>%
  dplyr::select(-Series)

# For the figure in the publication, we draw on the annual version of the Gillreath-Brown et al. 2024 temperature reconstruction to fit with the Moberg et al. 2005 reconstruction, which also has annual.
# temp_anom_standard2 <-
#   readr::read_csv(
#     "../Paper 1/supp_table_s2_full.csv",
#     show_col_types = FALSE
#   ) %>%
#   dplyr::mutate(Year = as.integer(Year)) %>%
#   dplyr::filter(Year >= 133 & Year <= 1900) %>%
#   dplyr::select(year = Year, anom = 2) %>%
#   dplyr::mutate(update = (anom - mean(anom)) / sd(anom))


Moberg_standard2 <-
  Moberg_standard %>%
  dplyr::mutate(update = (anom - mean(anom)) / sd(anom)) %>%
  dplyr::select(-Series)

Moberg_paleomat_diff <-
  temp_anom_standard2 %>%
  dplyr::select(-anom) %>%
  dplyr::left_join(., Moberg_standard2 %>% dplyr::select(-anom), by = "year") %>%
  dplyr::mutate(anom = update.x - update.y) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(anom3 = dplyr::case_when(
    update.y > update.x ~ "Lower",
    update.y < update.x ~ "Upper",
    TRUE ~ NA
  ))

# Moberg LF and paleomat
Moberg_paleomat_diff_fig <- Moberg_paleomat_diff %>%
  ggplot2::ggplot(aes(x = year, y = anom)) +
  geom_hline(
    yintercept = 0,
    color = "darkgrey",
    lty = 3,
    lwd = 1.25
  ) +
  geom_line(size = 2.0) +
  xlab("Years AD") +
  ylab("ENSO Index") +
  scale_x_continuous(
    limits = c(0, 2000),
    breaks = seq(0, 2000, 200),
    minor_breaks = seq(0, 2000, by = 100),
    guide = "axis_minor",
    labels = c(1, seq(200, 2000, 200))
  ) +
  scale_y_continuous(limits = c(-3, 4),
                     breaks = seq(-3, 4, 1.0)) +
  annotate(
    "text",
    x = 400,
    y = 2.0,
    label = "La Niña-like",
    size = 12
  ) +
  annotate(
    "text",
    x = 400,
    y = -2.0,
    label = "El Niño-like",
    size = 12
  ) +
  theme_bw() +
  guides(color = guide_legend(override.aes = list(fill = NA), direction = "horizontal")) +
  theme(
    legend.key.width = unit(1.25, 'cm'),
    legend.text = element_text(size = 19, family = "Helvetica"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.margin = margin(t = -14),
    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 = 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
    )
  )

Moberg_paleomat_diff_fig

ggplot2::ggsave(
  Moberg_paleomat_diff_fig,
  filename = here::here("vignettes/figures/Figure_3.png"),
  width = 17.25,
  height = 11.5,
  dpi = 600
)

```

```{r Figure_5, fig.cap="Figure 5. Earliest site with maize in each of the labeled political subdivisions, using data from the Ancient Maize Map (Blake et al. 2017) with additions listed in Supplemental Table S1. The linear model shows the regression of latitude on date for the labeled sites. The slope of the regression line implies an average spread rate of 0.001 degrees/year for this sample (r2 = 0.01; p > F = 0.72).", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 17.25, fig.height = 11.5, fig.align = "left"}

# Subset the dataframe to only have the 2 sigma data and the unique lab code.
MDB_min_fig4 <- MDB %>%
  cropDiffusionR::extract_earliest_crop() %>%
  # Remove any sites that are to the south of Guila Naquitz. You can also replace this with another site name if you have a different starting point for a crop. If the site is not present (e.g., Guila Naquitz is not in the p3k14C data), then this will not filter out any sites.
  dplyr::filter(
    if ("Guila Naquitz" %in% SiteName)
      Lat >= MDB %>% dplyr::filter(SiteName == "Guila Naquitz") %>% dplyr::pull(Lat) %>% dplyr::first()
    else
      TRUE
  )

# Write output to file, but without location data.
# tmp4 <- MDB_min_fig4 %>%
#   dplyr::select(-c(Long, Lat))
# readr::write_rds(tmp4, here::here("vignettes/output/earliest_maize_province.rds"))

# Plot using a linear model.
figure5 <-
  MDB_min_fig4 %>%
  dplyr::mutate(Province = as.factor(Province)) %>%
  dplyr::arrange(Date) %>%
  dplyr::mutate(labels = paste0(
    stringr::str_replace(SiteName, " \\s*\\([^\\)]+\\)", ""),
    ", \n ",
    Province
  )) %>%
  ggplot(aes(x = Date, y = Lat)) +
  geom_point(size = 4, aes(shape = forcats::fct_rev(Country))) +
  scale_shape_manual(values = c(19, 15), name = "Country") +
  geom_smooth(
    aes(group = 1),
    method = 'lm',
    se = FALSE,
    linetype = 1,
    color = "black"
  ) +
  ggrepel::geom_label_repel(
    aes(label = labels),
    max.iter = 4500,
    cex = 6.7,
    show.legend  = F,
    min.segment.length = unit(0, 'lines'),
    force = 25,
    box.padding = 1
  ) +
  scale_x_continuous(
    breaks = seq(-4500, 1500, 500),
    minor_breaks = seq(-4500, 1500, 100),
    limits = c(-4500, 1500),
    guide = "axis_minor",
    labels = c(seq(-4500,-500, 500), 1, seq(500, 1500, 500))
  ) +
  scale_y_continuous(
    breaks = seq(10, 45, 5),
    minor_breaks = seq(10, 45, 1),
    limits = c(10, 45),
    guide = "axis_minor"
  ) +
  xlab("Years BC/AD") +
  ylab("Latitude") +
  # ggpubr::color_palette(palette = "jco") +
  theme_bw()  +
  guides(color = guide_legend(override.aes = list(fill = NA), direction = "horizontal")) +
  theme(
    legend.title = element_text(
      size = 19,
      family = "Helvetica",
      face = "bold"
    ),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.text = element_text(size = 18, colour = "black"),
    legend.margin = margin(t = -14),
    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 = 26,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 7,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 7,
        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
    )
  )

figure5

# ggplot2::ggsave(
#   figure5,
#   filename = here::here("vignettes/figures/Figure_5.png"),
#   width = 17.25,
#   height = 11.5,
#   dpi = 600
# )

```

```{r Figure_6, fig.cap="Figure 6. Location of maize’s northern frontier through time. Shown is the earliest site with maize in each century, with the constraint that only sites further north than earlier sites appear in the plot. Data from the Ancient Maize Map (Blake et al. 2017) with additions listed in Supplemental Table S1. The linear model shows the regression of latitude on date for the labeled sites. The slope of the regression line implies an average spread rate of 0.003 degrees/year for this sample (r2 = 0.76; p > F < 0.001). Between just Guilà Naquitz and Las Capas however the implied rate of spread is ~8.7 times more rapid (0.026 degrees/year).", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 17.25, fig.height = 11.5, fig.align = "left"}

# Subset the dataframe to only have the 2 sigma data and the unique lab code.
MDB_min_fig5 <- MDB %>%
  cropDiffusionR::extract_diffusion(direction = "north") %>%
  dplyr::filter(Date <= 1500) %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName))

# Write output to file, but without location data.
# tmp5 <- MDB_min_fig5 %>%
#   dplyr::select(-c(Long, Lat))
# readr::write_rds(tmp5, here::here("vignettes/output/maize_frontier_mexico.rds"))

# Plot using a linear model.
figure6 <-
  MDB_min_fig5 %>%
  dplyr::mutate(Province = forcats::fct_reorder(Province, Date, .fun = 'min')) %>%
  dplyr::arrange(Date) %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName)) %>%
  dplyr::mutate(labels = stringr::str_replace(SiteName, " \\s*\\([^\\)]+\\)", "")) %>%
  ggplot(aes(x = Date, y = Lat)) +
  geom_point(size = 4, aes(color = Province, shape = Province)) +
  geom_smooth(
    aes(group = 1),
    method = 'lm',
    se = FALSE,
    linetype = 1,
    color = "black"
  ) +
  # stat_poly_eq(aes(label = paste(after_stat(eq.label),
  #                              after_stat(rr.label), sep = "*\", \"*"))) +
  ggrepel::geom_label_repel(
    aes(colour = Province, label = labels),
    max.iter = 4500,
    cex = 6.7,
    show.legend  = F,
    min.segment.length = unit(0, 'lines'),
    force = 25,
    box.padding = 1,
    max.overlaps = 16
  ) +
  scale_x_continuous(
    breaks = seq(-4500, 1500, 500),
    minor_breaks = seq(-4500, 1500, 100),
    limits = c(-4500, 1500),
    guide = "axis_minor",
    labels = c(seq(-4500, -500, 500), 1, seq(500, 1500, 500))
  ) +
  scale_y_continuous(
    breaks = seq(15, 45, 5),
    minor_breaks = seq(15, 45, 1),
    limits = c(15, 45),
    guide = "axis_minor"
  ) +
  xlab("Years BC/AD") +
  ylab("Latitude") +
  scale_shape_manual(
    values = c(3, 19, 15, 17, 9),
    name = "States",
    labels = c("Oxaca", "Arizona", "New Mexico", "Utah", "Colorado")
  ) +
  scale_color_manual(
    values = ggsci::pal_jco("default")(10)[c(1, 3, 2, 4, 6)],
    name = "States",
    labels = c("Oxaca", "Arizona", "New Mexico", "Utah", "Colorado")
  ) +
  theme_bw()  +
  theme(
    legend.position = "bottom",
    legend.title = element_text(
      size = 19,
      family = "Helvetica",
      face = "bold"
    ),
    legend.direction = "horizontal",
    legend.text = element_text(size = 18, colour = "black"),
    legend.margin = margin(t = -14),
    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 = 26,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 7,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 7,
        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
    )
  )

figure6

# ggplot2::ggsave(
#   figure6,
#   filename = here::here("vignettes/figures/Figure_6.png"),
#   width = 17.25,
#   height = 11.5,
#   dpi = 600
# )

```


```{r Figure_7, fig.cap="Figure 7. Location of maize’s northern frontier through time within or near the SWUS as defined here. The slope of the regression line implies an average spread rate of 0.002 degrees/year for this sample (r2 = 0.95; p > F < 0.001).", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 17.25, fig.height = 11.5, fig.align = "left"}

# Subset the dataframe to only have the 2 sigma data and the unique lab code.
MDB_min_fig6 <- MDB %>%
  dplyr::filter(Country == "USA") %>%
  cropDiffusionR::extract_diffusion(direction = "north") %>%
  dplyr::filter(Date <= 1500) %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName))

# Write output to file, but without location data.
# tmp6 <- MDB_min_fig6 %>%
#   dplyr::select(-c(Long, Lat))
# readr::write_rds(tmp6, here::here("vignettes/output/maize_frontier_swus.rds"))

# Plot using a linear model.
figure7 <-
  MDB_min_fig6 %>%
  dplyr::mutate(Province = forcats::fct_reorder(Province, Date, .fun = 'min')) %>%
  dplyr::arrange(Date) %>%
  dplyr::mutate(SiteName = ifelse(SiteName == "Elsinore Burial", SiteID, SiteName)) %>%
  dplyr::mutate(labels = stringr::str_replace(SiteName, " \\s*\\([^\\)]+\\)+$", "")) %>%
  ggplot(aes(x = Date, y = Lat)) +
  geom_point(size = 4, aes(color = Province, shape = Province)) +
  geom_smooth(
    aes(group = 1),
    method = 'lm',
    se = FALSE,
    linetype = 1,
    color = "black"
  ) +
  # stat_poly_eq(aes(label = paste(after_stat(eq.label),
  #                              after_stat(rr.label), sep = "*\", \"*"))) +
  ggrepel::geom_label_repel(
    aes(colour = Province, label = labels),
    max.iter = 4500,
    cex = 6.7,
    show.legend  = F,
    min.segment.length = unit(0, 'lines'),
    force = 25,
    box.padding = 1,
    max.overlaps = 15
  ) +
  scale_shape_manual(
    values = c(19, 15, 17, 9),
    name = "States",
    labels = c("Arizona", "New Mexico", "Utah", "Colorado")
  ) +
  scale_color_manual(
    values = ggsci::pal_jco("default")(10)[c(3, 2, 4, 6)],
    name = "States",
    labels = c("Arizona", "New Mexico", "Utah", "Colorado")
  ) +
  scale_x_continuous(
    breaks = seq(-4500, 1500, 500),
    minor_breaks = seq(-4500, 1500, 100),
    limits = c(-4500, 1500),
    guide = "axis_minor",
    labels = c(seq(-4500, -500, 500), 1, seq(500, 1500, 500))
  ) +
  scale_y_continuous(
    breaks = seq(30, 42, 2),
    minor_breaks = seq(30, 42, 1),
    limits = c(30, 42),
    guide = "axis_minor"
  ) +
  xlab("Years BC/AD") +
  ylab("Latitude") +
  theme_bw()  +
  theme(
    legend.position = "bottom",
    legend.title = element_text(
      size = 19,
      family = "Helvetica",
      face = "bold"
    ),
    legend.direction = "horizontal",
    legend.text = element_text(size = 18, colour = "black"),
    legend.margin = margin(t = -14),
    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 = 26,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 7,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 28,
      family = "Helvetica",
      margin = margin(
        t = 7,
        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
    )
  )

figure7

# ggplot2::ggsave(
#   figure7,
#   filename = here::here("vignettes/figures/Figure_7.png"),
#   width = 17.25,
#   height = 11.5,
#   dpi = 600
# )

```

```{r Figure_4, fig.cap="Figure 4. Annual accumulated growing degree days across the southwestern United States and Mexico, with locations of sites mentioned in the text or the following figures. Growing degree days were extracted and calculated from the 30 seconds (~1 km2) average temperature (°C) from WorldClim version 2.1 (https://www.worldclim.org/data/worldclim21.html). We use 10°C as the base and 30°C as the maximum for the GDD calculation.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.width = 15.43, fig.height = 13.15, fig.align = "left"}

# Placing Figure 3 after Figures 4-6 since we use the sites from those figures to plot on to the GDD map.
NASW_sites <-
  dplyr::bind_rows(MDB_min_fig4, MDB_min_fig5, MDB_min_fig6) %>%
  dplyr::ungroup() %>%
  dplyr::select(SiteName, SiteID, Long, Lat) %>%
  dplyr::mutate(SiteName = dplyr::coalesce(SiteName, SiteID)) %>%
  dplyr::distinct(SiteName, .keep_all = TRUE) %>%
  sf::st_as_sf(coords = c("Long", "Lat"), crs = 4326)

# Get raster for the North American Southwest and convert to a dataframe.
NASW_gdd_data <-
  raster::as.data.frame(cropDiffusionR::WorldClim_annual_gdd, 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")
  )

figure4 <-
  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 GDDs") +
  coord_equal() +
  theme_bw() +
  geom_sf(
    data = cropDiffusionR::swus_mexico_states,
    color = "#2b2b2b",
    fill = "transparent",
    size = 0.4,
    inherit.aes = FALSE
  ) +
  geom_sf(
    data = NASW_sites,
    aes(geometry = geometry),
    inherit.aes = FALSE,
    size = 4
  ) +
  ggrepel::geom_label_repel(
    data = NASW_sites,
    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
  ) +
  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"
    )
  )

figure4

# ggplot2::ggsave(
#   figure4,
#   filename = here::here("vignettes/figures/Figure_4.png"),
#   width = 15.43,
#   height = 13.15,
#   dpi = 900
# )

```

# Statistics for Figures 5–7

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

stats_table <-
  purrr::map2_dfr(
    list(.x = MDB_min_fig4, MDB_min_fig5, MDB_min_fig6),
    .y = list("Figure 5", "Figure 6", "Figure 7"),
    .f = function(x, y) {
      lm(Lat ~ Date, data = x) %>%
        summary() %>%
        {
          data.frame(
            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")
          )
        } %>%
        dplyr::mutate(Figure = as.character(y)) %>%
        dplyr::relocate(Figure, .before = "slope")
    }
  )

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

```



```{r SuppFig1, fig.cap="Supplemental Figure S1. Every sample in our maize database is plotted on this graph by date and latitude. Additionally, we include the modern annual accumulated growing degree days for the site location of the sample. Growing degree days were extracted and calculated from the 30 seconds (~1 km2) average temperature (°C) from WorldClim version 2.1 (https://www.worldclim.org/data/worldclim21.html). We use 10°C as the base and 30°C as the maximum for the GDD calculation.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.height = 7.73, fig.align = "left", eval = FALSE}

# Subset the dataframe to only have the 2 sigma data and the unique lab code.
MDB_min_rayshader <- MDB %>%
  dplyr::select(LabID, SiteName, Date, Lat, everything()) %>%
  dplyr::ungroup() %>%
  dplyr::filter(
    Lat >= MDB %>% dplyr::filter(SiteName == "Guila Naquitz") %>% dplyr::pull(Lat) %>% dplyr::first()
  ) %>%
  dplyr::left_join(., MDB_gdd, by = "LabID")

# Using rayshader package to show all maize samples and the modern GDD at each sample location.
mtplot <-
  MDB_min_rayshader %>%
  dplyr::mutate(Province = as.factor(Province)) %>%
  dplyr::arrange(Date) %>%
  ggplot(aes(x = Date, y = Lat, color = GDD)) +
  geom_point(size = 2) +
  scale_colour_gradientn(limits = c(0, 7000), colors = colorRampPalette(rev(RColorBrewer::brewer.pal(8, "RdBu")))(150)) +
  scale_x_continuous(
    breaks = c(seq(-4500,-500, 1000), 1, seq(500, 1500, 1000)),
    minor_breaks = c(seq(-4000,-1000, 1000), 1000),
    limits = c(-4500, 1500),
    guide = "axis_minor",
    labels = c(seq(-4500,-500, 1000), 1, seq(500, 1500, 1000))
  ) +
  scale_y_continuous(
    breaks = seq(10, 45, 5),
    minor_breaks = seq(10, 45, 1),
    limits = c(10, 45),
    guide = "axis_minor"
  ) +
  xlab("Years BC/AD") +
  ylab("Latitude") +
  labs(colour = "GDD") +
  theme_bw()  +
  theme(
    legend.title = element_text(
      size = 14,
      family = "Helvetica",
      face = "bold"
    ),
    legend.title.align = 1.0,
    legend.position = "right",
    legend.direction = "vertical",
    legend.text = element_text(size = 12, colour = "black"),
    legend.margin = margin(t = -14),
    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 = 12,
      colour = "black",
      family = "Helvetica"
    ),
    axis.title.y = element_text(
      size = 14,
      family = "Helvetica",
      margin = margin(
        t = 10,
        r = 7,
        b = 10,
        l = 10
      )
    ),
    axis.title.x = element_text(
      size = 14,
      family = "Helvetica",
      margin = margin(
        t = 7,
        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
    )
  )

rayshader::plot_gg(
  mtplot,
  width = 6.5,
  height = 4,
  multicore = TRUE,
  windowsize = c(1000, 1000),
  zoom = 0.77,
  phi = 50,
  theta = 35,
  sunangle = 225,
  soliddepth = -0.5,
  pointcontract = 0.7,
  scale = 200
)
rayshader::render_snapshot(
  clear = TRUE,
  filename = here::here("vignettes/figures/rayshader_gdd.png")
)

```


```{r SuppVidS1, fig.cap="Supplemental Video S1. Each maize site with a maize sample is plotted on a map of modern annual accumulated growing degree days. Samples were binned into 100 year intervals for plotting. Sites are accumulated across all centuries; therefore, the last century shows all sites with samples in our database. Many centuries have no recorded samples; thus, they are not plotted here. Growing degree days were extracted and calculated from the 30 seconds (~1 km2) average temperature (°C) from WorldClim version 2.1 (https://www.worldclim.org/data/worldclim21.html). We use 10°C as the base and 30°C as the maximum for the GDD calculation.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.height = 7.73, fig.align = "left", eval = FALSE}

# Here, we do not evaluate this code block, as it can take a long time to create due to the GDD raster. 
NASW_sites2 <- MDB %>%
  dplyr::select(LabID, SiteName, SiteID, Long, Lat, Date) %>%
  dplyr::mutate(SiteName = dplyr::coalesce(SiteName, SiteID)) %>%
  sf::st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  dplyr::group_by(period = cut(
    Date,
    breaks = seq(-4300, 1900, 100),
    labels = seq(-4300, 1800, 100)
  )) %>%
  # Only label sites that show up in Figures 3-6.
  dplyr::mutate(
    period = as.numeric(as.character(period)),
    label = ifelse(SiteName %in% NASW_sites$SiteName, SiteName, NA_character_)
  ) %>%
  dplyr::arrange(period) %>%
  dplyr::distinct(geometry, .keep_all = TRUE) %>%
  split(.$period)

# We create an accumulation of sites with samples across all centuries.
NASW_sites2_cumulative <- dplyr::bind_rows(NASW_sites2) %>%
  dplyr::filter(period < 1600) %>%
  dplyr::mutate(label = ifelse(LabID %in% MDB_min_fig5$LabID, label, NA_character_)) %>%
  tidyr::complete(period = seq(-4300, 1500, 100)) %>%
  tidyr::fill(., c(LabID:label), .direction = "down") %>%
  # Geometry column will not fill so need to group by lab id, then fill in with first geometry value.
  dplyr::group_by(LabID) %>%
  dplyr::mutate(geometry = geometry[1]) %>%
  dplyr::ungroup() %>%
  dplyr::arrange(period) %>%
  dplyr::group_by(period) %>%
  dplyr::distinct(SiteName, .keep_all = TRUE) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(Long = unlist(map(.$geometry, 1)),
                Lat = unlist(map(.$geometry, 2))) %>%
  split(.$period) %>%
  # Bind rows to the previous period.
  purrr::accumulate(., function(a, b) {
    bind_rows(a, b)
  }) %>%
  # Keep only unique sites in each period, since they are accumulated.
  purrr::map(., function(x)
    x %>%
      dplyr::distinct(SiteName, .keep_all = TRUE) %>%
      dplyr::mutate(label = ifelse(label == SiteName, label, NA_character_)) %>%
      dplyr::arrange(Date) %>%
      dplyr::mutate(
        counter = ifelse(!is.na(label), TRUE, NA_real_),
        counter = row_number(counter)
      ))

NASW_sites2_cumulative_names <-
  data.frame(century = as.numeric(names(NASW_sites2_cumulative))) %>%
  dplyr::mutate(
    century = ifelse(century == 0, 1, century),
    century = ifelse(
      century < 0,
      paste(century, "to", century + 100, "BC"),
      paste("AD", century, "to", century + 100)
    ),
    century = stringr::str_replace_all(century, "[-]", "")
  ) %>%
  dplyr::pull(century)

# Plot the maps for when there is only 1 site (i.e., that have no segments).
maize_through_time_gdd2_1 <-
  purrr::map2(NASW_sites2_cumulative[1:6], NASW_sites2_cumulative_names[1:6], function(x, y) {
    ggplot2::ggplot(data = NASW_gdd_data_transformed) +
      geom_raster(aes(x = long, y = lat, fill = GDD),
                  alpha = 0.8,
                  na.rm = TRUE) +
      scale_fill_gradientn(colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(255),
                           na.value = "transparent",
                           name = "Annual GDDs") +
      coord_equal() +
      theme_bw() +
      geom_sf(
        data = cropDiffusionR::swus_mexico_states,
        color = "#2b2b2b",
        fill = "transparent",
        size = 0.4,
        inherit.aes = FALSE
      ) +
      geom_sf(
        data = x,
        aes(geometry = geometry),
        size = 4,
        inherit.aes = FALSE
      ) +
      ggrepel::geom_label_repel(
        data = x,
        aes(label = label, geometry = geometry),
        stat = "sf_coordinates",
        min.segment.length = 0.4,
        force = 40,
        label.size = 0,
        size = 4.5,
        inherit.aes = FALSE,
        max.overlaps = 20
      ) +
      ggtitle(y) +
      theme(
        #legend.key.width = unit(1.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          size = 20,
          family = "Helvetica",
          face = "bold"
        ),
        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"
        )
      )
  })

# Plot sites that have segments (i.e., more than 2 sites).
maize_through_time_gdd2_2 <-
  purrr::map2(NASW_sites2_cumulative[7:59], NASW_sites2_cumulative_names[7:59], function(x, y) {
    route <- purrr::discard(x$counter, is.na)
    
    df_seg <- x %>%
      dplyr::filter(!is.na(label))
    
    df_seg <-
      data.frame(
        lat1 = df_seg$Lat[route[-length(route)]],
        long1 = df_seg$Long[route[-length(route)]],
        lat2 = df_seg$Lat[route[2:length(route)]],
        long2 = df_seg$Long[route[2:length(route)]]
      )
    
    ggplot2::ggplot(data = NASW_gdd_data_transformed) +
      geom_raster(aes(x = long, y = lat, fill = GDD),
                  alpha = 0.8,
                  na.rm = TRUE) +
      scale_fill_gradientn(colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(255),
                           na.value = "transparent",
                           name = "Annual GDDs") +
      coord_equal() +
      theme_bw() +
      geom_sf(
        data = cropDiffusionR::swus_mexico_states,
        color = "#2b2b2b",
        fill = "transparent",
        size = 0.4,
        inherit.aes = FALSE
      ) +
      geom_sf(
        data = x,
        aes(geometry = geometry),
        size = 4,
        inherit.aes = FALSE
      ) +
      geom_segment(
        data = df_seg,
        aes(
          x = long1,
          y = lat1,
          xend = long2,
          yend = lat2
        ),
        arrow = arrow(angle = 12, type = "closed"),
        linetype = 1,
        size = 1
      ) +
      ggrepel::geom_label_repel(
        data = x,
        aes(label = label, geometry = geometry),
        stat = "sf_coordinates",
        min.segment.length = 0.4,
        force = 40,
        label.size = 0,
        size = 4.5,
        inherit.aes = FALSE,
        max.overlaps = 20
      ) +
      ggtitle(y) +
      theme(
        #legend.key.width = unit(1.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          size = 20,
          family = "Helvetica",
          face = "bold"
        ),
        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"
        )
      )
  })

maize_through_time_gdd2 <-
  c(maize_through_time_gdd2_1, maize_through_time_gdd2_2)

invisible(
  mapply(
    ggsave,
    file = paste0(
      stringr::str_remove(here::here("Paper 2/animation2/"), "cropDiffusionR/"),
      purrr::map_chr(
        1:length(maize_through_time_gdd2),
        ~ stringr::str_pad(., 4, pad = "0")
      ),
      ".png"
    ),
    width = 15.43,
    height = 13.15,
    dpi = 300,
    plot = maize_through_time_gdd2
  )
)

#Gif animation
# list file names and read in
imgs2 <-
  list.files(stringr::str_remove(here::here("Paper 2/animation2/"), "cropDiffusionR/"),
             full.names = TRUE)
img_list2 <- lapply(imgs2, magick::image_read)

# join the images together
img_joined2 <- magick::image_join(img_list2)

# animate at 4 frames per second
img_animated2 <- magick::image_animate(img_joined2, fps = 1)

# save to disk
magick::image_write(image = img_animated2,
                    path = stringr::str_remove(
                      here::here(
                        "Paper 2/animation2/maize_anomaly_plots_arranged_cumulative_1fps.gif"
                      ),
                      "cropDiffusionR/"
                    ))

```


```{r SuppVidS2, fig.cap="Supplemental Video S2. Each maize site with a maize sample is plotted on a map of modern annual accumulated growing degree days. Samples were binned into 100 year intervals for plotting. Many centuries have no recorded samples; thus, they are not plotted here. Growing degree days were extracted and calculated from the 30 seconds (~1 km2) average temperature (°C) from WorldClim version 2.1 (https://www.worldclim.org/data/worldclim21.html). We use 10°C as the base and 30°C as the maximum for the GDD calculation.", echo = FALSE, warning = FALSE, message = FALSE, error = FALSE, results = 'hide',fig.keep = 'all', fig.height = 7.73, fig.align = "left", eval = FALSE}

NASW_sites2_names <-
  data.frame(century = as.numeric(names(NASW_sites2))) %>%
  dplyr::mutate(
    century = ifelse(century == 0, 1, century),
    century = ifelse(
      century < 0,
      paste(century, "to", century + 100, "BC"),
      paste("AD", century, "to", century + 100)
    ),
    century = stringr::str_replace_all(century, "[-]", "")
  ) %>%
  dplyr::pull(century)


maize_through_time_gdd <-
  purrr::map2(NASW_sites2, NASW_sites2_names, function(x, y)
    ggplot2::ggplot(data = NASW_gdd_data_transformed) +
      geom_raster(
        aes(x = long, y = lat, fill = GDD),
        alpha = 0.8,
        na.rm = TRUE
      ) +
      scale_fill_gradientn(
        colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(255),
        na.value = "transparent",
        name = "Annual GDDs"
      ) +
      coord_equal() +
      theme_bw() +
      geom_sf(
        data = cropDiffusionR::swus_mexico_states,
        color = "#2b2b2b",
        fill = "transparent",
        size = 0.4,
        inherit.aes = FALSE
      ) +
      geom_sf(
        data = x,
        aes(geometry = geometry),
        size = 4,
        inherit.aes = FALSE
      ) +
      ggrepel::geom_label_repel(
        data = x,
        aes(label = label, geometry = geometry),
        stat = "sf_coordinates",
        min.segment.length = 0.4,
        force = 47,
        label.size = 0,
        size = 4.5,
        inherit.aes = FALSE,
        max.overlaps = 20
      ) +
      ggtitle(y) +
      theme(
        #legend.key.width = unit(1.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          size = 20,
          family = "Helvetica",
          face = "bold"
        ),
        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"
        )
      ))

invisible(
  mapply(
    ggsave,
    file = paste0(
      stringr::str_remove(here::here("Paper 2/animation/"), "cropDiffusionR/"),
      purrr::map_chr(
        1:length(maize_through_time_gdd),
        ~ stringr::str_pad(., 4, pad = "0")
      ),
      ".png"
    ),
    width = 15.43,
    height = 13.15,
    dpi = 300,
    plot = maize_through_time_gdd
  )
)

# Plot empty centuries.
NASW_sites2_names_empty <-
  data.frame(century = setdiff(seq(-4300, 1800, 100), as.numeric(names(NASW_sites2)))) %>%
  dplyr::mutate(
    century = ifelse(century == 0, 1, century),
    century = ifelse(
      century < 0,
      paste(century, "to", century + 100, "BC"),
      paste("AD", century, "to", century + 100)
    ),
    century = stringr::str_replace_all(century, "[-]", "")
  ) %>%
  dplyr::pull(century)

maize_through_time_gdd2 <-
  purrr::map(NASW_sites2_names_empty, function(x)
    ggplot2::ggplot(data = NASW_gdd_data_transformed) +
      geom_raster(
        aes(x = long, y = lat, fill = GDD),
        alpha = 0.8,
        na.rm = TRUE
      ) +
      scale_fill_gradientn(
        colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(255),
        na.value = "transparent",
        name = "Annual GDDs"
      ) +
      coord_equal() +
      theme_bw() +
      geom_sf(
        data = cropDiffusionR::swus_mexico_states,
        color = "#2b2b2b",
        fill = "transparent",
        size = 0.4,
        inherit.aes = FALSE
      ) +
      ggtitle(x) +
      theme(
        #legend.key.width = unit(1.25, 'cm'),
        plot.title = element_text(
          hjust = 0.5,
          size = 20,
          family = "Helvetica",
          face = "bold"
        ),
        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"
        )
      ))

invisible(
  mapply(
    ggsave,
    file = paste0(
      stringr::str_remove(here::here("Paper 2/animation/empty/"), "cropDiffusionR/"),
      purrr::map_chr(
        1:length(maize_through_time_gdd2),
        ~ stringr::str_pad(., 4, pad = "0")
      ),
      ".png"
    ),
    width = 15.43,
    height = 13.15,
    dpi = 300,
    plot = maize_through_time_gdd2
  )
)

#Gif animation
# list file names and read in
imgs <-
  list.files(stringr::str_remove(here::here("Paper 2/animation"), "cropDiffusionR/"),
             full.names = TRUE)

#Here, we use gifski, as it is a lot quicker than magick.
gifski(
  imgs,
  gif_file = stringr::str_remove(
    here::here("Paper 2/animation/maize_anomaly_plots_arranged.gif"),
    "cropDiffusionR/"
  ),
  width = 1117,
  height = 952,
  delay = 1
)

# For using magick to create gif.
# img_list <- lapply(imgs, magick::image_read)
#
# # join the images together
# img_joined <- magick::image_join(img_list)
#
# # animate at 4 frames per second
# img_animated <- magick::image_animate(img_joined, fps = 1)
#
# # save to disk
# magick::image_write(image = img_animated,
#                     path = stringr::str_remove(
#                       here::here("Paper 2/animation/maize_anomaly_plots_arranged.gif"),
#                       "cropDiffusionR/"
#                     ))

```