knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = TRUE, error = TRUE)
knitr::opts_chunk$set(progress = TRUE, verbose = TRUE)
knitr::opts_chunk$set(
  echo = TRUE,
  fig.width = 9,
  fig.height = 9,
  fig.align = "center"
)
# Note that GhostScript needs to be installed at the system level for the PS files to be generated.
# MacOS users can use `brew install ghostscript`
# 
# Windows users can follow these directions:
# 1.    Go to the GhostScript website (https://www.ghostscript.com/download/gsdnld.html)
# 2.    Download the windows installer suitable for your machine
# 3.    Run the installer file which you downloaded and follow the prompts
# 4.    After running the installer click the windows "Start" button and type "Edit environment variables for your account" and open
# 5.    In the tab 'Advanced' click the button at the bottom 'Environment Variables...'
# 6.    Under 'System variables' find the variable 'Path', select 'Path' and click the 'Edit' button
# 7.    Select a new line and copy the Ghostscript 'bin' folder location into the field.
# 7.1   If you installed ghostscript to the default folder location; then the folder location will likely be "C:\Program Files\gs\gs9.52\bin", the version number (9.52) may differ.
# 8.    Save and exit the environmental variables window
# This chunk is then run only if knitting on new computer that the files have not been generated on
# this is necessary to embed fonts in .eps files 
library("extrafont")
if (.Platform$OS.type == "windows") {
   font_import(pattern = "arial", prompt = FALSE)
   loadfonts(device = "postscript", quiet = TRUE)
} else {
   font_import(pattern = "Arial", prompt = FALSE)
   loadfonts(device = "postscript", quiet = TRUE)
}

Load libraries

library("readxl")
library("cowplot")
library("ggplot2")
library("ggpubr")
library("grDevices")
library("dplyr")
library("lubridate")
library("clifro")
library("viridis")
library("showtext")
library("readr")
library("here")
library("patchwork")
library("classInt")
library("extrafont")

extrafont::loadfonts()

Import data

dat <-
   read_excel(
      system.file("extdata", "SpatioTemporalSpreadData_N.xlsx",
                  package = "spatiotemporaldynamics"),
      sheet = 1
   )

Examine data

str(dat)

Convert variables to their correct classes

cols_1 <-
   c(
      "location",
      "distance",
      "plot_number",
      "quadrat",
      "min_temp",
      "max_temp",
      "min_rh",
      "max_rh",
      "avg_wind_speed",
      "avg_rh",
      "avg_temp",
      "assessment_number",
      "total_rain"
   )
dat[cols_1] <- lapply(dat[cols_1], factor)

cols_2 <- c("infected_plants", "total_plants")
dat[cols_2] <- lapply(dat[cols_2], as.integer)
dat$assessment_date <- as.Date(dat$assessment_date)

Re-check class

sapply(dat, class)

Kernel denisty plots

Kernel density plot to visualise data distribution as by distance & location

ggplot(dat, aes(infected_plants)) +
   geom_density() +
   facet_grid(distance ~ location, labeller = label_both) +
   theme_pubclean(base_family = "Arial Unicode MS")

Kernel density plot to visualise overall data distribution

ggplot(dat, aes(x = infected_plants)) +
  geom_density(fill= "steelblue", alpha = 0.7) +
  geom_rug(alpha = 0.4) +
  xlab("Infected plants count")

Line plot with median and max/min values

fig_1 <- ggplot(data = dat,
                mapping = aes(x = assessment_number, y = infected_plants)) +
   geom_pointrange(
      stat = "summary",
      fun.min = min,
      fun.max = max,
      fun = median
   ) +
   stat_summary(fun = median,
                geom = "line",
                aes(group = location)) +
   facet_grid(distance ~ location) +
   xlab("Assessment number") +
   ylab("Infected plants") +
   theme_pubclean(base_family = "Arial Unicode MS")

fig_1
ggsave(here::here("man", "figures/Fig1.png"))
ggsave(here::here("man", "figures/Fig1.eps"), device = cairo_ps)

Wind rose

Import wind direction data and covert wind direction from text to degrees

# data munging
wind_direc_dat <-
   read_csv(
      system.file("extdata", "WindDirectionData.csv",
                  package = "spatiotemporaldynamics")
   ) %>%
   mutate(
      wind_direction_degrees = case_when(
         wind_direction == "N" ~ "0",
         wind_direction == "NbE" ~ "11.25",
         wind_direction == "NNE" ~ "22.5",
         wind_direction == "NEbN" ~ "33.75",
         wind_direction == "NE" ~ "45",
         wind_direction == "NEbE" ~ "56.25",
         wind_direction == "ENE" ~ "67.5",
         wind_direction == "EbN" ~ "73.5",
         wind_direction == "E" ~ "90",
         wind_direction == "EbS" ~ "101.2",
         wind_direction == "ESE" ~ "112.5",
         wind_direction == "SEbE" ~ "123.8",
         wind_direction == "SE" ~ "135.1",
         wind_direction == "SEbS" ~ "146.3",
         wind_direction == "SSE" ~ "157.6",
         wind_direction == "SbE" ~ "168.8",
         wind_direction == "S" ~ "180",
         wind_direction == "SbW" ~ "191.2",
         wind_direction == "SSW" ~ "202.5",
         wind_direction == "SWbS" ~ "213.8",
         wind_direction == "SW" ~ "225",
         wind_direction == "SWbW" ~ "236.2",
         wind_direction == "WSW" ~ "247.5",
         wind_direction == "WbS" ~ "258.8",
         wind_direction == "W" ~ "270",
         wind_direction == "WbN" ~ "281.2",
         wind_direction == "WNW" ~ "292.5",
         wind_direction == "NWbW" ~ "303.8",
         wind_direction == "NW" ~ "315",
         wind_direction == "NWbN" ~ "326.2",
         wind_direction == "NNW" ~ "337.5",
         wind_direction == "NbW" ~ "348.8",
         TRUE ~ wind_direction
      )
   ) %>%
   mutate(wind_direction_degrees = as.numeric(wind_direction_degrees)) %>%
   mutate(date = dmy(date))

Join wind speed and wind direction data

# Import wind speed data
wind_speed_dat <-
   read_csv(system.file("extdata", "WindSpeedData.csv",
                        package = "spatiotemporaldynamics")) %>%
   mutate(date = as_date(dmy_hm(date)))

### Join wind speed and wind direction data
wind_dat <- left_join(wind_speed_dat, wind_direc_dat)

Plot windrose over chickpea growing season

fig_2 <-
   with(
      wind_dat,
      windrose(
         wind_speed,
         wind_direction_degrees,
         location,
         n_col = 2,
         legend_title = "Wind speed (m/s)"
      )
   )

fig_2 <-
   fig_2 +
   scale_fill_viridis_d(name = "Wind Speed (m/s)", direction = -1) +
   xlab("") +
   theme_pubclean(base_family = "Arial Unicode MS")

fig_2
ggsave(here::here("man", "figures/Fig2.png"))
ggsave(here::here("man", "figures/Fig2.eps"), device = cairo_ps)

Data munging for ploting windrose by assessment number

# data munging
wind_dt <-
   wind_dat %>%
   mutate(
      assessment_number =
         case_when(
            date <= "2020-07-02" & location == "Billa Billa" ~ 1,
            date > "2020-07-02" &
               date <= "2020-07-16" & location == "Billa Billa" ~ 2,
            date > "2020-07-16" &
               date <= "2020-07-30" & location == "Billa Billa" ~ 3,
            date > "2020-07-30" &
               date <= "2020-08-13" & location == "Billa Billa" ~ 4,
            date > "2020-08-13" &
               date <= "2020-08-26" & location == "Billa Billa" ~ 5,
            date > "2020-08-26" &
               date <= "2020-09-10" & location == "Billa Billa" ~ 6,
            date > "2020-09-10" &
               date <= "2020-09-24" & location == "Billa Billa" ~ 7,
            date > "2020-09-24" &
               date <= "2020-10-08" & location == "Billa Billa" ~ 8,
            date > "2020-10-08" &
               date <= "2020-10-22" & location == "Billa Billa" ~ 9,
            date > "2020-10-22" &
               date <= "2020-10-27" &
               location == "Billa Billa" ~ 10,
            date <= "2020-07-30" & location == "Tosari" ~ 1,
            date > "2020-07-30" &
               date <= "2020-08-14" & location == "Tosari" ~ 2,
            date > "2020-08-14" &
               date <= "2020-08-27" & location == "Tosari" ~ 3,
            date > "2020-08-27" &
               date <= "2020-09-10" & location == "Tosari" ~ 4,
            date > "2020-09-10" &
               date <= "2020-09-25" & location == "Tosari" ~ 5,
            date > "2020-09-25" &
               date <= "2020-10-09" & location == "Tosari" ~ 6,
            date > "2020-10-09" &
               date <= "2020-10-23" & location == "Tosari" ~ 7,
            date > "2020-10-23" &
               date <= "2020-11-05" & location == "Tosari" ~ 8
         )
   ) %>%
   mutate(assessment_number  = as.factor(assessment_number))


Billa_Billa_wind_dt <- filter(wind_dt, location == "Billa Billa")
Tosari_wind_dt <- filter(wind_dt, location == "Tosari") %>%
   droplevels()

Plot windrose

# create breaks for the windroses
breaks <-
   classIntervals(wind_dat$wind_speed, n = 4, style = "jenks")

# Billa Billa windrose
fig_3.1 <-
   with(
      Billa_Billa_wind_dt,
      windrose(
         wind_speed,
         wind_direction_degrees,
         facet = assessment_number,
         n_col = 5,
         legend_title = "Wind speed (m/s)",
         speed_cuts = c(0, 1.4, 2.8, 4.2, 5.6, 7)
      )
   )

fig_3.1 <-
   fig_3.1 +
   scale_fill_viridis_d(name = "Wind Speed (m/s)",
                        direction = -1,
                        option = "cividis") +
   scale_y_continuous(name = "Proportion (%)",
                      labels = c(0, 10, 20, 30, 40, 50, 60)) +
   xlab("") +
   theme_pubclean(base_family = "Arial Unicode MS") +
   theme(
      axis.ticks.length = unit(0, "mm"),
      axis.line = element_blank(),
      panel.spacing.x = unit(0, "lines"),
      panel.spacing.y = unit(0, "lines"),
      plot.margin = margin(0, 0, 0, 0, "cm")
   )

# Tosari windrose
fig_3.2 <-
   with(
      Tosari_wind_dt,
      windrose(
         wind_speed,
         wind_direction_degrees,
         facet = assessment_number,
         n_col = 5,
         legend_title = "Wind speed (m/s)",
         speed_cuts = c(0, 1.4, 2.8, 4.2, 5.6, 7)
      )
   )

fig_3.2 <-
   fig_3.2 +
   scale_fill_viridis_d(name = "Wind Speed (m/s)",
                        direction = -1,
                        option = "cividis") +
   scale_y_continuous(name = "Proportion (%)",
                      labels = c(0, 10, 20, 30, 40, 50)) +
   xlab("") +
   theme_pubclean(base_family = "Arial Unicode MS") +
   theme(
      legend.position = "none",
      axis.ticks.length = unit(0, "mm"),
      axis.line = element_blank(),
      panel.spacing.x = unit(0, "lines"),
      panel.spacing.y = unit(0, "lines"),
      plot.margin = margin(0, 0, 0, 0, "cm")
   )

fig_3 <-
   fig_3.1 / fig_3.2

fig_3 <-
   fig_3 +
   plot_annotation(tag_levels = "A") &
   theme(plot.tag = element_text())

fig_3
ggsave(here::here("man", "figures/Fig3.png"))
ggsave(here::here("man", "figures/Fig3.eps"), device = cairo_ps)

Colophon

sessionInfo()


IhsanKhaliq/spatiotemporaldynamics documentation built on July 21, 2021, 3:13 a.m.