Load your packages

### An alternative way of making sure your pacakges are installed and loaded
install.packages("pacman")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)
pacman::p_load(tidyverse, DT, sf, lubridate, vroom, janitor, visdat, mapview)
plot_desc <- read_csv("data/msh/MSH_PLOT_DESCRIPTORS.csv") %>% 
  janitor::clean_names()
## Parsed with column specification:
## cols(
##   PLOT_NAME = col_character(),
##   PLOT_CODE = col_character(),
##   FIRST_YEAR = col_double(),
##   LAST_YEAR = col_double(),
##   UTMGRID = col_character(),
##   UTMEAST = col_double(),
##   UTMNORTH = col_double(),
##   LONG = col_double(),
##   LAT = col_double(),
##   POT_RAD = col_double(),
##   HEAT_LOAD = col_double(),
##   ELEVATION_M = col_double(),
##   ASPECT = col_character(),
##   SLOPE_DEG = col_double(),
##   IMPACT_TYPE = col_character(),
##   SUCCESSION_TYPE = col_character()
## )
spp_desc <- read_csv("data/msh/MSH_SPECIES_DESCRIPTORS.csv") %>% 
  janitor::clean_names()
## Warning: Missing column names filled in: 'X17' [17], 'X18' [18], 'X19' [19]
## Parsed with column specification:
## cols(
##   Family = col_character(),
##   Species_authority = col_character(),
##   Common_name = col_character(),
##   `Synonymy _comment` = col_character(),
##   Raw_code = col_character(),
##   NLSPN_code = col_character(),
##   Succession_status = col_character(),
##   Higher_taxon = col_character(),
##   Duration = col_character(),
##   Growth_form = col_character(),
##   Life_form = col_character(),
##   Clonal = col_character(),
##   Dispersal = col_character(),
##   Origin = col_character(),
##   Wet = col_character(),
##   Location_disturbance = col_character(),
##   X17 = col_logical(),
##   X18 = col_logical(),
##   X19 = col_logical()
## )
spp_plot_yr <- read_csv("data/msh/MSH_SPECIES_PLOT_YEAR.csv") %>% 
  janitor::clean_names()
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   PLOT_ID_YEAR = col_character(),
##   PLOT_NAME = col_character()
## )
## See spec(...) for full column specifications.
## clean 
spp_desc <- spp_desc %>% 
  select(-starts_with("x"))

Plot description

Field names

plot_desc %>% 
  names() %>% 
  enframe() %>% 
  select(value) %>% 
  datatable()

plot types

plot_desc %>% 
  select(1:5) %>% 
  datatable(filter = "top")
plot_desc %>% 
  visdat::vis_dat() +
  labs(title = "plot description")

spp_desc %>% 
  visdat::vis_dat()+
  labs(title = "species description")

spp_plot_yr %>% 
  visdat::vis_dat() +
  labs(title = "species plot year") +
  coord_flip()

spp_plot_yr_tidy <- spp_plot_yr %>% 
  pivot_longer(cols = -(1:4),names_to = "spp", values_to = "cover")
spp_plot_yr_tidy %>% 
  janitor::tabyl(plot_name) %>% 
  datatable()
spp_plot_yr_tidy %>% 
  janitor::tabyl(year, plot_name) %>% 
  datatable()
## wrangling! need to add a plot code to join. Extractd from the first part of the plot_id_year
spp_plot_yr_tidy <- spp_plot_yr_tidy %>% 
  mutate(plot_code = str_sub(plot_id_year,start = 1L, end = 6)) 

## JOIN in the plot info to the species table
spp_plot_yr_tidy <- spp_plot_yr_tidy %>% 
  select(-plot_name) %>% 
  left_join(.,plot_desc,by = "plot_code")

Let’s make things spatial!

plot_desc_sf <- plot_desc %>% 
  filter(!is.na(utmeast)) %>% 
  st_as_sf(., coords = c("utmeast", "utmnorth"), crs = 26910)

Let’s make a map of plots…

## set some options
mapviewOptions(basemaps = c("OpenTopoMap","Esri.WorldShadedRelief","Esri.WorldImagery"),
               raster.palette = grey.colors,
               vector.palette = colorRampPalette(c("snow", "cornflowerblue", "grey10")),
               na.color = "magenta",
               layers.control.pos = "topright")

plot_desc_sf %>% 
  mapview(zcol='elevation_m')

Look how easy we can change the variable that is symbolized.

plot_desc_sf %>% 
  mapview(zcol='succession_type')
## create spatial
spp_plot_yr_tidy_sf <- spp_plot_yr_tidy %>% 
  filter(!is.na(utmeast)) %>% 
  st_as_sf(., coords = c("utmeast", "utmnorth"), crs = 26910)

## the missing data
# spp_plot_yr_tidy %>% 
#   filter(!is.na(utmeast))
# spp_plot_yr_tidy_sf %>% 
#   sample_n(1000) %>% 
#   mapview()
mean_cover_by_impact_type <- spp_plot_yr_tidy %>%
  filter(cover != 0) %>% 
  group_by(year, spp, impact_type) %>% 
  summarise(mean.cover = mean(cover, na.rm = TRUE)) %>% 
  ungroup()

mean_cover_by_succession_type <- spp_plot_yr_tidy %>%
  filter(cover != 0) %>% 
  group_by(year, spp, succession_type) %>% 
  summarise(mean.cover = mean(cover, na.rm = TRUE)) %>% 
  ungroup()

Here’s a fig showing mean cover by spp, year, and succession type

mean_cover_by_succession_type %>% 
  filter(!is.na(succession_type)) %>% 
  mutate(year = as.character(year)) %>% 
  ggplot(aes(year, spp)) +
  geom_tile(aes(fill = mean.cover), color = "white") +
  scale_fill_viridis_c() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45)) +
  facet_wrap(~succession_type, ncol=3)

Here’s a fig showing mean cover by species by year and impact type

Not superuseful, but a start… Just filter on two types to demo.

mean_cover_by_impact_type %>% 
  filter(impact_type == "Tephra" | impact_type == "Blast") %>% 
  filter(!is.na(impact_type)) %>% 
  mutate(year = as.character(year)) %>% 
  ggplot(aes(year, spp)) +
  geom_tile(aes(fill = mean.cover), color = "white") +
  scale_fill_viridis_c() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45)) +
  facet_wrap(~impact_type, ncol=2)