07_temporal-sampling.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

# Checking Temporal Sampling Frequency

How often are checklists recorded in each grid cell?

## Load libraries

```{r}
# load libraries
library(tidyverse)
library(sf)

# for plotting
library(ggplot2)
library(colorspace)
library(scico)
library(ggthemes)
library(ggspatial)
library(patchwork)
```

## Load checklist data

Here we load filtered checklist data and convert to UTM 43N coordinates.

```{r}
# load checklist data
load("results/02_data_prelim_processing.rdata")

# get checklists
data <- distinct(
  dataGrouped, sampling_event_identifier, observation_date,
  longitude, latitude
)

# remove old data
rm(dataGrouped)

# transform to UTM 43N
data <- st_as_sf(data, coords = c("longitude", "latitude"), crs = 4326)
data <- st_transform(data, crs = 32643)

# get coordinates and bind to data
data <- cbind(
  st_drop_geometry(data),
  st_coordinates(data)
)

# bin to 1000m
data <- mutate(data,
  X = plyr::round_any(X, 2500),
  Y = plyr::round_any(Y, 2500)
)
```

## Get time differences per grid cell

```{r}
# get time differences in days
data <- mutate(data, observation_date = as.POSIXct(observation_date))
data <- nest(data, data = c("sampling_event_identifier", "observation_date"))

# map over data
data <- mutate(data,
  lag_metrics = lapply(data, function(df) {
    df <- arrange(df, observation_date)

    lag <- as.numeric(diff(df$observation_date, na.rm = TRUE) / (24 * 3600))

    data <- tibble(
      mean_lag = mean(lag, na.rm = TRUE),
      median_lag = median(lag, na.rm = TRUE),
      sd_lag = sd(lag, na.rm = TRUE),
      n_chk = nrow(df)
    )

    data
  })
)
```

```{r}
# unnest lag metrics
data_lag <- select(data, -data)
data_lag <- unnest(data_lag, cols = "lag_metrics")

# set the mean and median to infinity if nchk is 1
data_lag <- mutate(data_lag,
  mean_lag = ifelse(n_chk == 1, Inf, mean_lag),
  median_lag = ifelse(n_chk == 1, Inf, median_lag),
  sd_lag = ifelse(n_chk == 1, Inf, sd_lag)
)

# set all 0 to 1
data_lag <- mutate(data_lag,
  mean_lag = mean_lag + 1,
  median_lag = median_lag + 1
)
# melt data by tile
# data_lag = pivot_longer(data_lag, cols = c("mean_lag", "median_lag", "sd_lag"))
```

## Time Since Previous Checklist

### Get aux data

```{r}
# hills data
wg <- st_read("data/spatial/hillsShapefile/Nil_Ana_Pal.shp") %>%
  st_transform(32643)

roads <- st_read("data/spatial/roads_studysite_2019/roads_studysite_2019.shp") %>%
  st_transform(32643)

# add land
library(rnaturalearth)
land <- ne_countries(
  scale = 50, type = "countries", continent = "asia",
  country = "india",
  returnclass = c("sf")
) %>%
  st_transform(32643)

bbox <- st_bbox(wg)
```

### Histogram of lags

Figure code hidden in HTML and PDF versions.

```{r}
# get lags
data <- mutate(data,
  lag_hist = lapply(data, function(df) {
    df <- arrange(df, observation_date)

    lag <- as.numeric(diff(df$observation_date, na.rm = TRUE) / (24 * 3600))

    data <- tibble(
      lag = lag + 1,
      index = seq(lag)
    )

    data
  })
)

# unnest lags
data_hist <- select(data, X, Y, lag_hist) %>%
  unnest(cols = "lag_hist")
```

```{r echo=FALSE}
fig_hist_lag <-
  ggplot(data_hist) +
  geom_histogram(
    aes(x = lag),
    bins = 10, size = 0.2, fill = "steelblue"
  ) +
  scale_x_log10(
    breaks = c(1, 30, 365),
    labels = c("1 day", "1 mo.", "1 yr.")
  ) +
  scale_y_continuous(
    labels = scales::label_number(
      scale = 0.001, accuracy = 1,
      suffix = "K"
    ),
    limits = c(0, 10.5e3)
  ) +
  theme_few() +
  theme(
    plot.background = element_rect(fill = "white", colour = 1),
    panel.background = element_blank(),
    panel.border = element_blank(), axis.line = element_blank(),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    )
  ) +
  coord_cartesian(
    expand = F
  ) +
  labs(
    x = "Timelag prev. obs.",
    y = "# checklists"
  )
```

```{r echo=FALSE}
# make plot
fig_lag_temporal <-
  ggplot(data_lag) +
  geom_sf(data = land, fill = "grey90", col = NA) +
  geom_sf(
    data = wg,
    fill = NA,
    lty = 2
  ) +
  annotation_custom(
    grob = fig_hist_lag %>% ggplotGrob(),
    xmin = bbox["xmax"] - (bbox["xmax"] - bbox["xmin"]) / 2.5,
    xmax = bbox["xmax"],
    ymin = bbox["ymax"] - (bbox["ymax"] - bbox["ymin"]) / 3,
    ymax = bbox["ymax"]
  ) +
  geom_tile(
    aes(X, Y,
      fill = mean_lag
    ),
    col = "grey90"
  ) +
  geom_sf(
    data = roads, size = 0.2, col = "indianred"
  ) +
  scale_fill_scico(
    palette = "davos",
    trans = "log10",
    begin = 0.1,
    limits = c(1, 365),
    breaks = c(1, 7, 30, 180, 365),
    label = c("1 day", "1 wk.", "1 mo.", "6 mo.", "> 1 yr."),
    na.value = alpha("ivory")
  ) +
  annotation_north_arrow(
    location = "bl", which_north = "true",
    pad_x = unit(0.1, "in"), pad_y = unit(0.5, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    text_cex = 1,
    style = "ticks"
  ) +
  coord_sf(
    expand = FALSE,
    xlim = bbox[c("xmin", "xmax")],
    ylim = bbox[c("ymin", "ymax")]
  ) +
  theme_test() +
  theme(
    legend.position = c(0.85, 0.5),
    legend.background = element_blank(),
    legend.key = element_rect(fill = "grey90"),
    legend.key.width = unit(2, units = "mm"),
    legend.title = element_text(face = "bold"),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    ),
    axis.title = element_blank(),
    panel.background = element_rect(fill = "lightblue")
  ) +
  labs(
    fill = "Mean time\nb/w checklists"
  )
```

```{r echo=FALSE}
# save figure
ggsave(
  fig_lag_temporal,
  filename = "figs/fig_temporal_bias.png",
  device = png()
)
```

![Most sites are resurveyed at least once, but some are visited much more frequently than others. There does not appear to be a link between roads and visit frequency. eBird checklists are also strongly clustered in time, with some of the most sampled areas over the study period visited at intervals of > 1 week, and with some less intensively sampled areas visited frequently, at intervals of < 1 week. Overall, the majority of checklists are reported only a day after the previous checklist at that location (see inset).](figs/fig_temporal_bias.png)

## Main Text Figure 3

Combining figures for spatial and temporal clustering into main text figure 3. This overall figure is not shown here, see main text.

```{r}
# load spatial bias figure
load("data/fig_checklists_grid.Rds")
```

```{r echo=FALSE}
# make figure
fig_sampling_bias <-
  wrap_plots(
    fig_checklists_grid, fig_lag_temporal
  ) +
    plot_annotation(
      tag_levels = "a",
      tag_prefix = "(",
      tag_suffix = ")"
    ) &
    theme(
      plot.tag = element_text(
        size = 14,
        family = "Arial",
        face = "bold"
      )
    )

# save figure
ggsave(
  fig_sampling_bias,
  filename = "figs/fig_03.png",
  width = 300,
  height = 150,
  units = "mm"
)
```

![**Distribution of sampling effort in the form of eBird checklists in the Nilgiri and Anamalai Hills between 2013 and 2021.**
(a) Sampling effort across the Nilgiri and Anamalai Hills, in the form of eBird checklists reported by birdwatchers, mostly takes place along roads, with the majority of checklists located <1 km from a roadway (see distribution in inset), and therefore, only about 300m, on average, from the location of another checklist. (b)  eBird checklists are also strongly clustered in time, with some of the most sampled areas over the study period visited at intervals of > 1 week, and with some less intensively sampled areas visited frequently, at intervals of < 1 week. Overall, most checklists are reported only a day after the previous checklist at that location (see inset). Both spatial and temporal clustering make data thinning necessary. Both panels show counts or mean intervals in a 2.5km grid cell; the study area is bounded by a dashed line, and roads within it are shown as (a) blue or (b) red lines.](figs/fig_03.png) 

## Checklists per Month

We counted the checklists per month, pooled over years, to determine how sampling effort varies over the year.

```{r}
# get two week period by date
data <- select(data, X, Y, data)

# unnest
data <- unnest(data, cols = "data")

# get fortnight
library(lubridate)
data <- mutate(data,
  week = week(observation_date),
  week = plyr::round_any(week, 2),
  year = year(observation_date),
  month = month(observation_date)
)

# count checklists per fortnight
data_count <- count(data, month, year)
```

```{r echo=FALSE}
ggplot(data_count) +
  geom_boxplot(
    aes(
      x = factor(month),
      y = n
    ),
    fill = "steelblue"
  ) +
  scale_y_log10(
    limits = c(10, NA)
  ) +
  theme_classic() +
  theme(
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    )
  ) +
  labs(
    x = "Month",
    y = "# checklists"
  )

ggsave(filename = "figs/fig_chk_per_month.png")
```

![Observations peak in the early months of the year, and decline towards the rainy months, slowly increasing until the following winter.](figs/fig_chk_per_month.png)