06_spatial-sampling-bias-distRoads.Rmd

---
editor_options: 
  chunk_output_type: console
---

# Examining Spatial Sampling Bias

The goal of this section is to show how far each checklist location is from the nearest road, and how far each site is from its nearest neighbour. This follows finding the pairwise distance between a large number of unique checklist locations to a vast number of roads, as well as to each other.

## Prepare libraries

```{r setup_libs}
# load libraries
# for data
library(sf)
library(rnaturalearth)
library(dplyr)
library(readr)
library(purrr)

# for plotting
library(scales)
library(ggplot2)
library(ggspatial)
library(colorspace)

# round any function
round_any <- function(x, accuracy = 20000) {
  round(x / accuracy) * accuracy
}
# ci function
ci <- function(x) {
  qnorm(0.975) * sd(x, na.rm = TRUE) / sqrt(length(x))
}
```

## Read checklist data

Read in checklist data with distance to nearest neighbouring site, and the distance to the nearest road.

```{r read_chk_dists}
# read from local file
chkCovars <- read_csv("data/03_data-covars-perChklist.csv")
```

### Spatially explicit filter on checklists

We filter the checklists by the boundary of the study area. This is *not* the extent.

```{r spatial_filter_nnb_dist}
chkCovars <- st_as_sf(chkCovars, coords = c("longitude", "latitude")) %>%
  `st_crs<-`(4326) %>%
  st_transform(32643)

# read wg
wg <- st_read("data/spatial/hillsShapefile/Nil_Ana_Pal.shp") %>%
  st_transform(32643)
# get bounding box
bbox <- st_bbox(wg)

# spatial subset
chkCovars <- chkCovars %>%
  mutate(id = 1:nrow(.)) %>%
  filter(id %in% unlist(st_contains(wg, chkCovars)))
```

### Get background land for plotting

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

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

## Prepare Main Text Figure 3

### Prepare histogram of distance to roads

Figure code is hidden in versions rendered as HTML or PDF.

```{r plot_histogram, echo=FALSE}
# make histogram
hist_roads <-
  ggplot(chkCovars) +
  geom_histogram(
    aes(
      dist_road / 1e3
    ),
    bins = 20, size = 0.2,
    fill = "steelblue"
  ) +
  scale_x_log10(
    # label = label_number(accuracy = 0.1),
    breaks = c(0.1, 1, 10),
    labels = c("0.1", "1", "10")
  ) +
  scale_y_continuous(
    breaks = c(0, 2500, 5000),
    label = label_number(
      scale = 0.001,
      accuracy = 1, suffix = "K"
    )
  ) +
  coord_cartesian(
    expand = F
  ) +
  theme_test() +
  theme(
    plot.background = element_rect(fill = "white", colour = 1),
    panel.background = element_blank(),
    panel.border = element_blank(), axis.line = element_blank()
  ) +
  labs(
    x = "Distance to roads (km)",
    y = "# Locations"
  )
```

### Table: Distance to roads

```{r save_global_mean_dist_roads}
# write the mean and ci95 to file
chkCovars %>%
  st_drop_geometry() %>%
  dplyr::select(dist_road, nnb) %>%
  tidyr::pivot_longer(
    cols = c("dist_road", "nnb"),
    names_to = "variable"
  ) %>%
  group_by(variable) %>%
  summarise_at(
    vars(value),
    list(~ mean(.), ~ sd(.), ~ min(.), ~ max(.))
  ) %>%
  write_csv("data/results/distance_roads_sites.csv")
```

```{r show_dist_roads_nnb, , message=FALSE, echo=FALSE}
# read in and show
library(magrittr)
readr::read_csv("data/results/distance_roads_sites.csv") %>%
  kableExtra::kbl(
    booktabs = TRUE
  ) %>%
  kableExtra::kable_styling(latex_options = c("hold_position"))
```

### Distance to nearest neighbouring site

```{r get_nn_site}
# get unique locations from checklists
locs_unique <- cbind(
  st_drop_geometry(chkCovars),
  st_coordinates(chkCovars)
) %>%
  as_tibble()

locs_unique <- distinct(locs_unique, X, Y, .keep_all = T)
```

Figure code is hidden in versions rendered as HTML and PDF.

```{r echo=FALSE, eval=FALSE}
# make histogram of nearest neighbours
hist_sites <-
  ggplot(locs_unique) +
  geom_histogram(aes(nnb / 1e3),
    bins = 100, size = 0.2, fill = "steelblue"
  ) +
  labs(x = "dist. nearest site (km)", y = "# sites") +
  # scale_x_log10(label=label_number(accuracy = 0.1),
  #               breaks = c(0.1, 1, 10))+
  coord_cartesian(xlim = c(0, 10)) +
  scale_y_continuous(label = label_number(
    scale = 0.001, accuracy = 1,
    suffix = "K"
  )) +
  theme_test() +
  theme(
    plot.background = element_rect(fill = NA, colour = 1),
    panel.background = element_blank(),
    panel.border = element_blank(), axis.line = element_blank()
  )
```

### Spatial distribution of distances to neighbours

Figure code is hidden in HTML and PDF versions, consult the Rmarkdown file.

```{r plot_map_nnb, echo=FALSE, eval=FALSE}
# transform points to utm
locs_unique <- locs_unique %>%
  st_as_sf(coords = c("X", "Y"), crs = 32643)

# add nnb to locations
fig_site_nnb <-
  ggplot() +
  geom_sf(data = land, fill = "grey90", col = NA) +
  geom_sf(data = wg, fill = NA, col = 1) +
  annotation_custom(
    grob = hist_sites %>% 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_sf(data = roads, size = 0.2, col = "steelblue") +
  geom_sf(
    data = locs_unique, aes(col = nnb / 1000),
    alpha = 0.5
  ) +
  scico::scale_colour_scico(
    palette = "lajolla",
    values = c(0, 1), direction = 1, limits = c(0, 5),
    na.value = "dodgerblue"
  ) +
  annotation_north_arrow(
    location = "br", which_north = "true",
    pad_x = unit(0.1, "in"), pad_y = unit(0.5, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  annotation_scale(location = "br", width_hint = 0.4, text_cex = 1) +
  theme_test() +
  theme(
    legend.position = c(0.9, 0.53),
    legend.background = element_blank(),
    legend.key = element_rect(fill = "grey90"),
    legend.key.width = unit(2, units = "mm"),
    legend.title = element_text(face = "bold"),
    axis.title = element_blank(),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    ),
    panel.background = element_rect(fill = "lightblue")
  ) +
  coord_sf(
    expand = FALSE, xlim = bbox[c("xmin", "xmax")],
    ylim = bbox[c("ymin", "ymax")]
  ) +
  labs(colour = "Dist NN site")

ggsave(
  fig_site_nnb,
  filename = "figs/fig_site_nnb.png",
  width = 6, height = 6
)
```

![Most observation sites are within 300m of another site.](figs/fig_site_nnb.png)

## Figure: Spatial sampling bias

```{r load_data_in_r}
# get locations
points <- chkCovars %>%
  bind_cols(as_tibble(st_coordinates(.))) %>%
  st_drop_geometry() %>%
  mutate(X = round_any(X, 2500), Y = round_any(Y, 2500))

# count points
points <- count(points, X, Y)
```

Figure code is hidden in versions rendered as HTML and PDF.

```{r plot_map_dist_roads, , echo=FALSE}
# plot on maps
fig_checklists_grid <-
  ggplot() +
  geom_sf(
    data = land,
    fill = "grey90",
    col = NA
  ) +
  geom_sf(
    data = wg,
    fill = NA,
    col = 1,
    lty = 2
  ) +
  annotation_custom(
    grob = hist_roads %>% 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(data = points, aes(X, Y, fill = n), col = "grey90") +
  geom_sf(data = roads, size = 0.2, col = "steelblue") +

  # scale_colour_manual(values = "steelblue", labels = "roads")+

  scale_fill_continuous_sequential(
    palette = "Lajolla",
    trans = "log10",
    rev = F
  ) +
  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"
  ) +
  theme_test() +
  theme(
    legend.position = c(0.9, 0.53),
    legend.background = element_blank(),
    legend.key = element_rect(fill = "grey90"),
    legend.key.width = unit(2, units = "mm"),
    legend.title = element_text(face = "bold"),
    axis.title = element_blank(),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    ),
    panel.background = element_rect(fill = "lightblue")
  ) +
  coord_sf(
    expand = FALSE,
    xlim = bbox[c("xmin", "xmax")],
    ylim = bbox[c("ymin", "ymax")]
  ) +
  labs(fill = "# Checklists", colour = NULL)
```

```{r}
# save as png
ggsave(
  fig_checklists_grid,
  filename = "figs/fig_spatial_bias.png"
)

# save figure as Robject for next plot
save(fig_checklists_grid, file = "data/fig_checklists_grid.Rds")
```

![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. Each cell here is 2.5km x 2.5km.](figs/fig_spatial_bias.png)