popovs/turkey_map.R

## turkey_map.R
# Interactive Wild Turkey (Meleagris gallopavo, B-WITU) map

# 00 SETUP ====

library(magrittr) # for %>%

# Set up reproducible renv environment so anyone using this code
# will achieve the same result, regardless of machine

# This means any libraries I download now will also be downloaded
# for any other user who uses this .Rproj file in future (no need
# to get them to run `install.packages("all_the_packages_please")`)
#install.packages("renv")
renv::activate()

# 01 DOWNLOAD DATA SOURCES ====

# There are three primary turkey datasets collated in this map:
# 1) iNaturalist sightings
# 2) Bird watcher surveys
# 3) Regional biologist reports

# In addition, this map incorporates BC regional spatial boundaries.

# 01-0 BC boundaries ----
# BC Boundaries can be downloaded from the BC Data Catalogue (https://catalogue.data.gov.bc.ca/dataset/natural-resource-nr-district)
# Rather than manually downloading the files and dumping them in the "data"
# folder (not reproducible), we can use R + WMS web service to download the data.

# XML: https://openmaps.gov.bc.ca/geo/pub/WHSE_ADMIN_BOUNDARIES.ADM_NR_DISTRICTS_SPG/ows?service=WMS&request=GetCapabilities
# Make URL object
wms_url <- "https://openmaps.gov.bc.ca/geo/pub/WHSE_ADMIN_BOUNDARIES.ADM_NR_DISTRICTS_SPG/ows"

# Open connection to the WMS url
wms_client <- ows4R::WFSClient$new(wms_url, serviceVersion = "1.3.0")

# Extract the name(s) of the data layers
layer <- wms_client$getFeatureTypes(pretty = TRUE)[[1]]

# And save the data as an sf object
BC <- wms_client$getCapabilities()$findFeatureTypeByName(layer)$getFeatures()

# No CRS by default - add it
sf::st_crs(BC) <- "EPSG:3005"

# The shapefile is very high-res and takes a looong time to load.
# Use rmapshaper to simplify features.
BC <- rmapshaper::ms_simplify(BC)

# Extract the bounds of BC in lat/long format
bbox <- sf::st_bbox(sf::st_transform(BC, 4326))

# Cleanup
rm(wms_client, wms_url, layer)

# 01-1 iNaturalist data ----

# Download data via `rinat` R API
inat <- rinat::get_inat_obs(query = "Meleagris gallopavo",
                            bounds = c(bbox$ymin, bbox$xmin, bbox$ymax, bbox$xmax), # S, W, N, E vector
                            maxresults = 1000,
                            geo = TRUE)

# Add source column
inat$source <- "iNaturalist"

# Cut out domesticated/captive turkeys, poor quality IDs
inat <- inat[inat$captive_cultivated == "false" & inat$quality_grade == "research" & inat$scientific_name != "Meleagris gallopavo domesticus",]

# Dropping irrelevant columns
# The goal here is to see turkey hotspots & geographic distribution
inat <- inat %>% dplyr::select(source,
                               scientific_name,
                               datetime,
                               place_guess,
                               latitude,
                               longitude,
                               url)

names(inat)[4] <- "location"

# Coerce datetime column
inat$datetime <- lubridate::as_datetime(inat$datetime, tz = "Canada/Pacific")

# Turn into sf object
inat <- sf::st_as_sf(inat,
                     coords = c("longitude", "latitude"),
                     crs = "EPSG:4326")

# Change CRS
inat <- sf::st_transform(inat, crs = "EPSG:3005")

# Only include sightings in BC
inat <- sf::st_intersection(inat, BC)

# 01-2 Bird watcher surveys ----

# Download birdwatcher .xls file from: http://a100.gov.bc.ca/pub/siwe/details.do?projectId=4212&surveyId=13509&pagerOffset=0

# URL for "Point data for the 1974-2001 surveys formatted for entry into SPI"
url <- "http://a100.gov.bc.ca/pub/siwe/download.do;jsessionid=9FC8DAC86955209E088C26E0D1F50915?docId=26426"

# Data needs to be downloaded to a temp file first, then unzipped,
# then finally imported into the R session
temp <- tempfile() # create temp download location
temp_dir <- tempfile() # create temp directory location
download.file(url, temp) # download the linked .zip to the temp location
unzip(zipfile = temp, exdir = temp_dir) # extract the zipped file to temp directory

# Read the .xls from the unzipped temp location (requires you to
# know the name of the downloaded file + sheet of interest)
# `sip` for 'species inventory program'
sip <- readxl::read_xls(file.path(temp_dir, "wsi_4212_dct.xls"),
                        sheet = "Location Behaviour and Sign",
                        na = c("", "NA", "N/A", "#N/A", "na"))

# Close temp connections
unlink(c(temp, temp_dir))
rm(url, temp, temp_dir)

# Clean up data
sip <- janitor::clean_names(sip)

# Filter down to only turkeys
sip <- sip[sip$species == "B-WITU", ]

# Remove empty records
sip <- janitor::remove_empty(sip, c("cols"))

# Adjust cols a bit to match inat cols
sip$location <- paste0(sip$study_area_name, " - ", sip$design_component_label)
sip$scientific_name <- "Meleagris gallopavo"
sip$source <- "BC Species Inventory Web Explorer 4212"

# Create datetime col
lubridate::date(sip$time) <- sip$date
attr(sip$time, "tzone") <- "Canada/Pacific" # Update timezone
names(sip)[4] <- "datetime"

# Going off the "Definitions of Fields and Codes"
# and "New Field Definitions" - dropping irrelevant columns
# The goal here is to see turkey hotspots & geographic distribution
sip <- sip %>% dplyr::select(source,
                             scientific_name,
                             datetime,
                             location,
                             count,
                             utm_easting,
                             utm_northing,
                             comments)

# Make sf object
# utm_zone col indicated that everything is 11N, i.e. EPSG:32611
sip <- sf::st_as_sf(sip, crs = "EPSG:32611", coords = c("utm_easting", "utm_northing"))

# Change CRS
sip <- sf::st_transform(sip, "EPSG:3005")

# Intersect with BC
sip <- sf::st_intersection(sip, BC)

# 01-3 Regional biologist reports ----

# Similar to the birdwatching data, we can download directly via R
url <- "https://www.env.gov.bc.ca/perl/soft/dl.pl/20221202132550-14-gp-4b15a990-2e54-45af-a2e2-fd2092ef?simple=y"
temp <- tempfile()
download.file(url, temp)
br <- readxl::read_excel(file.path(temp),
                         na = c("", "NA", "N/A", "#N/A", "na"),
                         trim_ws = TRUE)
unlink(temp)
rm(url)

# Clean up
br <- janitor::clean_names(br)
br %<>% dplyr::mutate_at(c("long", "lat"),
                         stringr::str_trim)
br %<>% dplyr::mutate_at(c("long", "lat"),
                         as.numeric)

# Adjust cols to match other two datasets
names(br)[4] <- "reporting_method" # change colname from 'source' to 'reporting_method'
br$source <- "Regional Biologist Report"
br$scientific_name <- "Meleagris gallopavo"

# Add datetime col
br$datetime <- lubridate::as_datetime(br$date)

# Drop irrelevant columns
br <- br %>% dplyr::select(source,
                           scientific_name,
                           datetime,
                           number_of_birds,
                           reporting_method,
                           long,
                           lat)

# Convert to sf object, transform, intersect
br <- br %>%
  sf::st_as_sf(crs = "EPSG:4326", coords = c("long", "lat")) %>%
  sf::st_transform(crs = "EPSG:3005") %>%
  sf::st_intersection(BC)

# 02 MERGE DATASETS ====

turkeys <- dplyr::bind_rows(inat, sip, br)
turkeys <- turkeys %>% dplyr::select(source:url, count:reporting_method, dplyr::everything())

# 03 MAKE MAP ====

library(ggplot2)
library(plotly)
library(ggsci)

bb <- sf::st_bbox(BC)

world <- rnaturalearth::ne_countries(scale = 50, returnclass = 'sf')
world <- sf::st_transform(world, crs = "EPSG:3005")
world <- sf::st_crop(world,
                     xmin = bb[[1]] - 5000,
                     ymin = bb[[2]] - 5000,
                     xmax = bb[[3]] + 5000,
                     ymax = bb[[4]] + 5000
                     )

bc_pal <- cartography::carto.pal(pal1 = "green.pal", n1 = 20, pal2 = "kaki.pal", n2 = 3)

p <- ggplot() +
  geom_sf(data = world, fill = grey(0.9)) +
  geom_sf(data = BC,
          aes(fill = DISTRICT_NAME),
          alpha = 0.3,
          show.legend = F) +
  scale_fill_manual(values = bc_pal) +
  geom_sf(data = turkeys, aes(color = source)) +
  scale_color_viridis_d() +
  #scale_color_lancet(name = "Data source") +
  ggtitle("Wild turkey (Meleagris gallopavo) sightings in BC",
          subtitle = "1985-2022") +
  theme_light() +
  guides(fill = "none")

ggplotly(p)
	# Interactive Wild Turkey (Meleagris gallopavo, B-WITU) map

	# 00 SETUP ====

	library(magrittr) # for %>%

	# Set up reproducible renv environment so anyone using this code
	# will achieve the same result, regardless of machine

	# This means any libraries I download now will also be downloaded
	# for any other user who uses this .Rproj file in future (no need
	# to get them to run `install.packages("all_the_packages_please")`)
	#install.packages("renv")
	renv::activate()

	# 01 DOWNLOAD DATA SOURCES ====

	# There are three primary turkey datasets collated in this map:
	# 1) iNaturalist sightings
	# 2) Bird watcher surveys
	# 3) Regional biologist reports

	# In addition, this map incorporates BC regional spatial boundaries.

	# 01-0 BC boundaries ----
	# BC Boundaries can be downloaded from the BC Data Catalogue (https://catalogue.data.gov.bc.ca/dataset/natural-resource-nr-district)
	# Rather than manually downloading the files and dumping them in the "data"
	# folder (not reproducible), we can use R + WMS web service to download the data.

	# XML: https://openmaps.gov.bc.ca/geo/pub/WHSE_ADMIN_BOUNDARIES.ADM_NR_DISTRICTS_SPG/ows?service=WMS&request=GetCapabilities
	# Make URL object
	wms_url <- "https://openmaps.gov.bc.ca/geo/pub/WHSE_ADMIN_BOUNDARIES.ADM_NR_DISTRICTS_SPG/ows"

	# Open connection to the WMS url
	wms_client <- ows4R::WFSClient$new(wms_url, serviceVersion = "1.3.0")

	# Extract the name(s) of the data layers
	layer <- wms_client$getFeatureTypes(pretty = TRUE)[[1]]

	# And save the data as an sf object
	BC <- wms_client$getCapabilities()$findFeatureTypeByName(layer)$getFeatures()

	# No CRS by default - add it
	sf::st_crs(BC) <- "EPSG:3005"

	# The shapefile is very high-res and takes a looong time to load.
	# Use rmapshaper to simplify features.
	BC <- rmapshaper::ms_simplify(BC)

	# Extract the bounds of BC in lat/long format
	bbox <- sf::st_bbox(sf::st_transform(BC, 4326))

	# Cleanup
	rm(wms_client, wms_url, layer)

	# 01-1 iNaturalist data ----

	# Download data via `rinat` R API
	inat <- rinat::get_inat_obs(query = "Meleagris gallopavo",
	bounds = c(bbox$ymin, bbox$xmin, bbox$ymax, bbox$xmax), # S, W, N, E vector
	maxresults = 1000,
	geo = TRUE)

	# Add source column
	inat$source <- "iNaturalist"

	# Cut out domesticated/captive turkeys, poor quality IDs
	inat <- inat[inat$captive_cultivated == "false" & inat$quality_grade == "research" & inat$scientific_name != "Meleagris gallopavo domesticus",]

	# Dropping irrelevant columns
	# The goal here is to see turkey hotspots & geographic distribution
	inat <- inat %>% dplyr::select(source,
	scientific_name,
	datetime,
	place_guess,
	latitude,
	longitude,
	url)

	names(inat)[4] <- "location"

	# Coerce datetime column
	inat$datetime <- lubridate::as_datetime(inat$datetime, tz = "Canada/Pacific")

	# Turn into sf object
	inat <- sf::st_as_sf(inat,
	coords = c("longitude", "latitude"),
	crs = "EPSG:4326")

	# Change CRS
	inat <- sf::st_transform(inat, crs = "EPSG:3005")

	# Only include sightings in BC
	inat <- sf::st_intersection(inat, BC)

	# 01-2 Bird watcher surveys ----

	# Download birdwatcher .xls file from: http://a100.gov.bc.ca/pub/siwe/details.do?projectId=4212&surveyId=13509&pagerOffset=0

	# URL for "Point data for the 1974-2001 surveys formatted for entry into SPI"
	url <- "http://a100.gov.bc.ca/pub/siwe/download.do;jsessionid=9FC8DAC86955209E088C26E0D1F50915?docId=26426"

	# Data needs to be downloaded to a temp file first, then unzipped,
	# then finally imported into the R session
	temp <- tempfile() # create temp download location
	temp_dir <- tempfile() # create temp directory location
	download.file(url, temp) # download the linked .zip to the temp location
	unzip(zipfile = temp, exdir = temp_dir) # extract the zipped file to temp directory

	# Read the .xls from the unzipped temp location (requires you to
	# know the name of the downloaded file + sheet of interest)
	# `sip` for 'species inventory program'
	sip <- readxl::read_xls(file.path(temp_dir, "wsi_4212_dct.xls"),
	sheet = "Location Behaviour and Sign",
	na = c("", "NA", "N/A", "#N/A", "na"))

	# Close temp connections
	unlink(c(temp, temp_dir))
	rm(url, temp, temp_dir)

	# Clean up data
	sip <- janitor::clean_names(sip)

	# Filter down to only turkeys
	sip <- sip[sip$species == "B-WITU", ]

	# Remove empty records
	sip <- janitor::remove_empty(sip, c("cols"))

	# Adjust cols a bit to match inat cols
	sip$location <- paste0(sip$study_area_name, " - ", sip$design_component_label)
	sip$scientific_name <- "Meleagris gallopavo"
	sip$source <- "BC Species Inventory Web Explorer 4212"

	# Create datetime col
	lubridate::date(sip$time) <- sip$date
	attr(sip$time, "tzone") <- "Canada/Pacific" # Update timezone
	names(sip)[4] <- "datetime"

	# Going off the "Definitions of Fields and Codes"
	# and "New Field Definitions" - dropping irrelevant columns
	# The goal here is to see turkey hotspots & geographic distribution
	sip <- sip %>% dplyr::select(source,
	scientific_name,
	datetime,
	location,
	count,
	utm_easting,
	utm_northing,
	comments)

	# Make sf object
	# utm_zone col indicated that everything is 11N, i.e. EPSG:32611
	sip <- sf::st_as_sf(sip, crs = "EPSG:32611", coords = c("utm_easting", "utm_northing"))

	# Change CRS
	sip <- sf::st_transform(sip, "EPSG:3005")

	# Intersect with BC
	sip <- sf::st_intersection(sip, BC)

	# 01-3 Regional biologist reports ----

	# Similar to the birdwatching data, we can download directly via R
	url <- "https://www.env.gov.bc.ca/perl/soft/dl.pl/20221202132550-14-gp-4b15a990-2e54-45af-a2e2-fd2092ef?simple=y"
	temp <- tempfile()
	download.file(url, temp)
	br <- readxl::read_excel(file.path(temp),
	na = c("", "NA", "N/A", "#N/A", "na"),
	trim_ws = TRUE)
	unlink(temp)
	rm(url)

	# Clean up
	br <- janitor::clean_names(br)
	br %<>% dplyr::mutate_at(c("long", "lat"),
	stringr::str_trim)
	br %<>% dplyr::mutate_at(c("long", "lat"),
	as.numeric)

	# Adjust cols to match other two datasets
	names(br)[4] <- "reporting_method" # change colname from 'source' to 'reporting_method'
	br$source <- "Regional Biologist Report"
	br$scientific_name <- "Meleagris gallopavo"

	# Add datetime col
	br$datetime <- lubridate::as_datetime(br$date)

	# Drop irrelevant columns
	br <- br %>% dplyr::select(source,
	scientific_name,
	datetime,
	number_of_birds,
	reporting_method,
	long,
	lat)

	# Convert to sf object, transform, intersect
	br <- br %>%
	sf::st_as_sf(crs = "EPSG:4326", coords = c("long", "lat")) %>%
	sf::st_transform(crs = "EPSG:3005") %>%
	sf::st_intersection(BC)

	# 02 MERGE DATASETS ====

	turkeys <- dplyr::bind_rows(inat, sip, br)
	turkeys <- turkeys %>% dplyr::select(source:url, count:reporting_method, dplyr::everything())

	# 03 MAKE MAP ====

	library(ggplot2)
	library(plotly)
	library(ggsci)

	bb <- sf::st_bbox(BC)

	world <- rnaturalearth::ne_countries(scale = 50, returnclass = 'sf')
	world <- sf::st_transform(world, crs = "EPSG:3005")
	world <- sf::st_crop(world,
	xmin = bb[[1]] - 5000,
	ymin = bb[[2]] - 5000,
	xmax = bb[[3]] + 5000,
	ymax = bb[[4]] + 5000
	)

	bc_pal <- cartography::carto.pal(pal1 = "green.pal", n1 = 20, pal2 = "kaki.pal", n2 = 3)

	p <- ggplot() +
	geom_sf(data = world, fill = grey(0.9)) +
	geom_sf(data = BC,
	aes(fill = DISTRICT_NAME),
	alpha = 0.3,
	show.legend = F) +
	scale_fill_manual(values = bc_pal) +
	geom_sf(data = turkeys, aes(color = source)) +
	scale_color_viridis_d() +
	#scale_color_lancet(name = "Data source") +
	ggtitle("Wild turkey (Meleagris gallopavo) sightings in BC",
	subtitle = "1985-2022") +
	theme_light() +
	guides(fill = "none")

	ggplotly(p)