library(data.table) # data manipulation

library(tidyverse) # data manipulation, date time parsing and plotting

library(sf) # spatial data handling and manipulation

library(here) # reproducible filepaths

library(lubridate)

species_code <- "RFB"

filepath <- here("Data", species_code)

filepattern <- "*.csv"

ls_filenames <- list.files(path = filepath, recursive = TRUE, pattern = filepattern)

IDstart <- 1 #start position of the ID in the filename

IDend <- 7 #end position of the ID in the filename

test <- fs::dir_ls(path = filepath, recurse = TRUE, type = "file", glob = filepattern)[1]

read_lines(test, n_max = 5) # change n_max to change the number of rows to read in

skiplines <- 0

date_formats <- c("dmY", "Ymd") #specify date formats (e.g. "dmY" works for 01-12-2022 and 01/12/2022)

datetime_formats <- c("dmY HMS", "Ymd HMS") #specify date & time format

trackingdatatimezone <- "GMT"

colnames <- TRUE

read_lines(test, n_max = 5, skip = skiplines)

user_delim <- ","

user_trim_ws <- TRUE # Should leading and trailing whitespace (ASCII spaces and tabs) be trimmed from each field before parsing it?

ID_type <- "ID"

df_combined <- fs::dir_ls(path = filepath, glob = filepattern, #use our defined filepath and pattern

type = "file", recurse = TRUE) %>% # recurse = T searches all sub-folders

purrr::set_names(nm = basename(.)) %>% # removing path prefix (makes filename more manageable)

purrr::map_dfr(read_delim, .id="filename", #read all the files in using filename as ID column

col_types = cols(.default = "c"), col_names = colnames,

skip = skiplines, delim = user_delim, trim_ws = user_trim_ws) %>%

mutate("{ID_type}" := str_sub(string = filename, start = IDstart, end = IDend), #substring ID from the filename (start to end of substring)

.after = filename) #position the new ID column after filename column

colnames(df_combined)

datetime_formats # a reminder of the datetime orders previously specified

datetime_colnames <- c("Date", "Time") # or c("DateTime")

df_slim <- data.frame(ID = as.character(df_combined$ID),

Date = df_combined$Date,

Time = df_combined$Time,

Y = as.numeric(df_combined$Latitude),

X = as.numeric(df_combined$Longitude))

str(df_slim);head(df_slim)

message("If you see any 'failed to parse' warnings below a date or time has not formatted (we will remove these NA's later)")

df_slim <- df_slim %>%

tidyr::unite(col = "DateTime_unparsed", all_of(datetime_colnames), sep = " ", remove = FALSE) %>%

mutate(DateTime = lubridate::parse_date_time(DateTime_unparsed, #use lubridate to parse DateTime

orders=datetime_formats, #using the datetime_formats object we made earlier

tz=trackingdatatimezone),

Date = lubridate::as_date(DateTime),

Year = lubridate::year(DateTime)) %>%

select(-DateTime_unparsed)

Fails <- df_slim %>% filter(is.na(DateTime)==T)

head(Fails)

df_raw <- df_slim %>%

arrange(across(all_of(c(ID_type, "DateTime")))) %>%

drop_na(DateTime) #remove NA's in datetime column

head(df_raw)

rm(list=ls()[!ls() %in% c("df_raw","date_formats","datetime_formats","trackingdatatimezone", "ID_type", "species_code")])

metadata_path <- paste0(species_code, "_Metadata.csv")

filepath_meta <- here("Data", metadata_path)

metadate_formats <- c("dmY", "Ymd") #specify date format used in metadata

metadatetime_formats <- c("dmY HMS", "Ymd HMS") #specify date & time format

metadatatimezone <- "Indian/Chagos" #specify timezone used for metadata

df_metadata <- readr::read_csv(filepath_meta)

names(df_metadata)

df_metadataslim <- data.frame(ID = as.character(df_metadata$BirdID), # compulsory column

TagID = as.character(df_metadata$TagID),

DeployID = as.character(df_metadata$DeployID),

DeployDate_local = df_metadata$DeploymentDate, # compulsory column (set to NA if not relevant)

DeployTime_local = df_metadata$DeploymentTime, # compulsory column (set to NA if not relevant)

RetrieveDate_local = df_metadata$RetrievalDate, # compulsory column (set to NA if not relevant)

RetrieveTime_local = df_metadata$RetrievalTime, # compulsory column (set to NA if not relevant)

CPY = df_metadata$NestLat,

CPX = df_metadata$NestLong,

Species = species_code,

Population = df_metadata$Population,

Age = df_metadata$Age,

BreedStage = df_metadata$BreedingStage)

df_metadataslim <- df_metadataslim %>%

mutate(Deploydatetime = parse_date_time(paste(DeployDate_local, DeployTime_local),#create deploy datetime

order=metadatetime_formats,

tz=metadatatimezone),

Retrievedatetime = parse_date_time(paste(RetrieveDate_local, RetrieveTime_local), #create retrieve datetime

order=metadatetime_formats,

tz=metadatatimezone)) %>%

select(-any_of(c("DeployDate_local","DeployTime_local", "RetrieveDate_local", "RetrieveTime_local"))) %>%

mutate(across(contains('datetime'), #return datetime as it would appear in a different tz

~with_tz(., tzone=trackingdatatimezone)))

df_temporalextents <- df_raw %>%

group_by(across(all_of(ID_type))) %>%

summarise(min_datetime = min(DateTime),

max_datetime = max(DateTime))

df_metadataslim <- df_metadataslim %>%

left_join(., df_temporalextents, by = ID_type) %>%

mutate(Deploydatetime = case_when(!is.na(Deploydatetime) ~ Deploydatetime,

is.na(Deploydatetime) ~ min_datetime),

Retrievedatetime = case_when(!is.na(Retrievedatetime) ~ Retrievedatetime,

is.na(Retrievedatetime) ~ max_datetime)) %>%

select(-c(min_datetime, max_datetime))

df_metamerged <- df_raw %>%

left_join(., df_metadataslim, by=ID_type)

head(df_metamerged)

rm(list=ls()[!ls() %in% c("df_metamerged", "species_code")]) #specify objects to keep

No_data_vals <- c(0, -999)

na_cols <- c("X", "Y", "DateTime", "ID")

df_clean <- df_metamerged %>%

drop_na(all_of(na_cols)) %>%

filter(!X %in% No_data_vals & !Y %in% No_data_vals) %>% # remove bad data values in Lat Lon columns

distinct(DateTime, ID, .keep_all = TRUE) %>% # this might be a problem for ACC data where we don't have milliseconds so beware if using it for this purpose

filter(case_when(!is.na(Retrievedatetime) ~ Deploydatetime < DateTime & DateTime < Retrievedatetime, # keep only data within deployment period

.default = Deploydatetime < DateTime)) # if retrieve date is NA (i.e., tags submit via satellite), only filter by deploy date

head(df_clean)

rm(list=ls()[!ls() %in% c("df_clean", "species_code")]) #specify objects to keep

tracking_crs <- 4326 # Only change if data are in a different coordinate system

meta_crs <- 4326 # Only change if data are in a different coordinate system

transform_crs <- 3857 # Suggest changing if using own data

df_diagnostic <- df_clean %>%

mutate(geometry_GPS = st_transform( # transform X/Y coordinates to `transform_crs` for distance calculations

st_as_sf(., coords=c("X","Y"), crs=tracking_crs), crs = transform_crs)$geometry) %>%

group_by(ID) %>% #back to grouping by ID for calculations per individual

mutate(dist = st_distance(geometry_GPS, lag(geometry_GPS), by_element = T), # distance travelled from previous fix, by_element = T means calculations are done by row

difftime = difftime(DateTime, lag(DateTime), units="secs"), # time passed since previous fix

netdisp = st_distance(geometry_GPS, geometry_GPS[1], by_element = F)[,1], # calculate distance between first location and current location, by_element = F returns dense matrix with all pairwise distances

speed = as.numeric(dist)/as.numeric(difftime), # calculate speed (distance/time)

dX = as.numeric(X)-lag(as.numeric(X)), #difference in longitude, relative to previous location

dY = as.numeric(Y)-lag(as.numeric(Y)), #difference in latitude, relative to previous location

turnangle = atan2(dX, dY)*180/pi + (dX < 0)*360) %>% #angle (in degrees) from previous to current location using formula theta = atan(y/x), where y = change along y axis & x = change along x axis

ungroup() %>% select(-c(geometry_GPS, dX, dY)) # ungroup and remove geometries

df_diagnostic <- st_coordinates(st_transform(st_as_sf(df_diagnostic, coords=c("X","Y"), crs=tracking_crs), crs = 4326)) %>%

as.data.frame() %>%

rename("Lon" = "X", "Lat" = "Y") %>%

cbind(df_diagnostic, .)

filepath_dfout <- here("DataOutputs","WorkingDataFrames")

dir.create(filepath_dfout)

filename_dfout <- paste0(species_code, "_diagnostic")

write_csv(df_diagnostic, file = here(filepath_dfout, paste0(filename_dfout,".csv")))

rm(list=ls()[!ls() %in% c("df_diagnostic", "species_code")]) #specify objects to keep

if (!require("shiny")) install.packages("shiny")

library(shiny)

runGitHub("ExMove", username = "ExMove",

ref = "master", subdir = "app")

filter_cutoff <- as.period(30, unit = "minutes")

filter_speed <- 20

filter_netdisp_dist <- 300

filter_netdist_units <- "km" # e.g., "m", "km"

filter_netdisp <- units::as_units(filter_netdisp_dist, filter_netdist_units)

df_filtered <- df_diagnostic %>%

filter(Deploydatetime + filter_cutoff < DateTime, # keep times after cutoff

speed < filter_speed, # keep speeds slower than speed filter

netdisp <= filter_netdisp) # keep distances less than net displacement filter

head(df_filtered)

rm(list=ls()[!ls() %in% c("df_diagnostic", "df_filtered", "species_code")]) #specify objects to keep

sampleRateUnits <- "mins"

grouping_factors_poplevel <- c("Species")

grouping_factors_indlevel <- c("ID")

se <- function(x) sqrt(sd(x, na.rm = T) / length(x[!is.na(x)]))

df_summary_ind <- df_filtered %>%

group_by(across(c(all_of(grouping_factors_poplevel), all_of(grouping_factors_indlevel)))) %>%

summarise(NoPoints = NROW(ID), # number of fixes

NoUniqueDates = length(unique(Date)), # number of tracking dates

FirstDate = as.Date(min(Date)), # first tracking date

LastDate = as.Date(max(Date)), # last tracking date

SampleRate = mean(as.numeric(difftime, units = sampleRateUnits), na.rm = T), # sample rate mean

SampleRate_se = se(as.numeric(difftime, units = sampleRateUnits))) # sample rate standard error

df_summary_pop <- df_summary_ind %>% # use the individual-level summary data

group_by(across(all_of(c(grouping_factors_poplevel)))) %>%

summarise(NoInds = length(unique(ID)), # number of unique individuals

NoPoints_total = sum(NoPoints), # total number of tracking locations

FirstDate = as.Date(min(FirstDate)), # first tracking date

LastDate = as.Date(max(LastDate)), # last tracking date

PointsPerBird = mean(NoPoints), # number of locations per individual: mean

PointsPerBird_se = se(NoPoints), # number of locations per individual: standard error

DatesPerBird = mean(NoUniqueDates), # number of tracking days per bird: mean

DatesPerBird_se = se(NoUniqueDates), # number of tracking days per bird: standard error

SampleRate_mean = mean(SampleRate), # sample rate mean

SampleRate_se = se(SampleRate)) # sample rate standard error

rm(list=ls()[!ls() %in% c("df_diagnostic", "df_filtered", "df_summary_ind", "df_summary_pop", "species_code")]) #specify objects to keep

filepath_filtered_out <- here("DataOutputs","WorkingDataFrames")

filepath_summary_out <- here("DataOutputs","SummaryDataFrames")

dir.create(filepath_filtered_out)

dir.create(filepath_summary_out)

filename_filtered_out <- paste0(species_code, "_filtered")

filename_summary_ind_out <- paste0(species_code, "_summary_ind")

filename_summary_pop_out <- paste0(species_code, "_summary_pop")

write_csv(df_filtered, file = here(filepath_filtered_out, paste0(filename_filtered_out,".csv")))

write_csv(df_summary_ind, file = here(filepath_summary_out, paste0(filename_summary_ind_out,".csv")))

write_csv(df_summary_pop, file = here(filepath_summary_out, paste0(filename_summary_pop_out,".csv")))

rm(list=ls()[!ls() %in% c("df_diagnostic", "df_filtered", "df_summary_ind", "df_summary_pop", "species_code")]) #specify objects to keep

device <- "tiff"

units <- "mm"

dpi <- 300

out_path <- here("DataOutputs","Figures")

dir.create(out_path)

topo_label = "Depth (m)"

message("If you see masking warning these are fine.

Watch out for packages that aren't installed yet")

library(rnaturalearth)

library(marmap)

library(plotly)

df_plotting <- df_filtered %>%

group_by(ID) %>%

mutate(diffsecs = as.numeric(difftime),

secs_elapsed = cumsum(replace_na(diffsecs, 0)),

time_elapsed = as.duration(secs_elapsed),

days_elapsed = as.numeric(time_elapsed, "days")) %>%

mutate(across(c(dist,speed, Lat, Lon), as.numeric))

countries <- ne_countries(scale = "medium", returnclass = "sf")

minlon <- min(df_plotting$Lon)

maxlon <- max(df_plotting$Lon)

minlat <- min(df_plotting$Lat)

maxlat <- max(df_plotting$Lat)

base_topography_map <- getNOAA.bathy(lon1 = minlon - 0.1, lon2 = maxlon + 0.1,

lat1 = minlat - 0.1, lat2 = maxlat + 0.1, resolution = 1)

base_topography_fort = fortify(base_topography_map)

map_base <- ggplot() +

geom_raster(data = base_topography_fort, aes(x=x, y=y, fill=z), alpha = 0.9) +

# add colour scheme for the fill

scale_fill_viridis_c(option="mako", name = topo_label) +

# add map of countries over the top

geom_sf(data = countries, aes(geometry = geometry), fill = NA) +

# set plot limits

coord_sf(xlim = c(minlon-0.1, maxlon+0.1),

ylim = c(minlat-0.1, maxlat+0.1), crs = 4326, expand = F) +

# add labels

labs(x = "Longitude", y = "Latitude") +

theme(axis.text=element_text(colour="black"),

axis.title.x = element_text(size = 15),

axis.text.x = element_text(hjust=0.7),

axis.title.y = element_text(angle=90, vjust = 0.4, size = 15),

axis.text.y = element_text(hjust=0.7, angle=90, vjust=0.3)) +

# set a theme

theme_light()

map_alllocs <- map_base +

# add GPS points

geom_point(data = df_plotting, aes(x = Lon, y = Lat), alpha = 0.8, size = 0.5, col = "violetred3")

map_individuals <- map_base +

# add GPS points and paths between them

geom_point(data = df_plotting, aes(x = Lon, y = Lat, col = speed),

alpha = 0.8, size = 0.5 ) +

geom_path(data = df_plotting, aes(x = Lon, y = Lat, col = speed),

alpha = 0.8, linewidth = 0.5 ) +

# colour birds using scale_colour_gradient2

scale_colour_gradient2(name = "Speed", low = "blue", mid = "white", high = "red",

midpoint = (max(df_plotting$speed,na.rm=TRUE)/2)) + # `midpoint` argument ensures an even transition of color across speed value

# facet for individual

facet_wrap(~ ID, ncol = round(sqrt(n_distinct(df_plotting$ID))))

ggsave(plot = map_alllocs,

filename = paste0(species_code, "_map_all_locs.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi,

)

ggsave(plot = map_individuals,

filename = paste0(species_code, "_map_individuals.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi,

)

speed_time_plot <- df_plotting %>% #speed over time

ggplot(data=., aes(x=days_elapsed, y=speed, group=ID)) +

# add line of speed over time

geom_line() +

# add axis labels

xlab("time elapsed (days)") + ylab("speed (m/s)") +

# facet by individual

facet_wrap(~ID, nrow= round(sqrt(n_distinct(df_plotting$ID)))) +

# set plotting theme

theme(axis.text=element_text(colour="black")) +

theme_light()

message("Warnings about 'non-finite' values for speed/step length plots are expected

non-finite values should be equal to number of individuals)")

ggsave(plot = speed_time_plot,

filename = paste0(species_code, "_speed_timeseries_plot.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi

)

speed_hist <- df_plotting %>%

ggplot(data=., aes(speed)) +

geom_histogram(binwidth=0.1, alpha=0.7) + # can adjust binwidth to suite your needs

geom_density(aes(y =0.1*after_stat(count))) +

# add plot labels

xlab("speed (m/s)") + ylab("count") +

# facet by individual

facet_wrap(~ID, nrow= round(sqrt(n_distinct(df_plotting$ID)))) +

# set plotting theme

theme(axis.text=element_text(colour="black"))+

theme_light()

ggsave(plot = speed_hist,

filename = paste0(species_code, "_speed_histogram.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi,

)

step_time_plot <- df_plotting %>% #step length over time

ggplot(data=., aes(x=days_elapsed, y=as.numeric(netdisp), group=ID)) +

geom_line() +

# add plot labels

xlab("time elapsed (days)") + ylab("Distance from first fix (m)") +

# facet by individual

facet_wrap(~ID, nrow= round(sqrt(n_distinct(df_plotting$ID)))) +

# set plotting theme

theme(axis.text=element_text(colour="black")) +

theme_light()

ggsave(plot = step_time_plot,

filename = paste0(species_code, "_step_time_plot.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi,

)

step_hist <- df_plotting %>% #step histogram

ggplot(data=., aes(as.numeric(dist))) +

geom_histogram(binwidth=5000, alpha=0.7) + # can adjust binwidth to suite your needs

# add plot labels

xlab("step length (m)") + ylab("count") +

# facet by individual

facet_wrap(~ID, nrow= round(sqrt(n_distinct(df_plotting$ID))))+

# set plotting theme

theme_light()

ggsave(plot = step_hist,

filename = paste0(species_code, "_step_hist.", device),

device = device,

path = out_path,

units = units, width = 200, height = 175, dpi = dpi,

)

rm(list=ls()[!ls() %in% c("df_diagnostic", "species_code")]) #specify objects to keep