Calculate spatio-temporal distance between points
Johannes Krietsch
2023-08-15
#==============================================================================================================Data and code from “Mutual mate guarding and limited sexual conflict
in a sex-role reversed shorebird” Contributor: Johannes Krietsch
📍 This script runs relative to the project’s root directory and
describes how I calculated the spatio-temporal distance between
points.
#==============================================================================================================
# Packages
sapply( c('data.table', 'magrittr', 'sdb', 'ggplot2', 'auksRuak', 'sf','foreach', 'doFuture', 'knitr'),
require, character.only = TRUE)## data.table magrittr sdb ggplot2 auksRuak sf foreach doFuture knitr
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# Functions
source('./R/0_functions.R')
# Lines to run to create html output
opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# rmarkdown::render('./R/3_spatio_temporal_distance.R', output_dir = './OUTPUTS/R_COMPILED')
# Projection
PROJ = '+proj=laea +lat_0=90 +lon_0=-156.653428 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 '
# Data
d = fread('./DATA/NANO_TAGS.txt', sep = '\t', header = TRUE, nThread = 20) %>% data.table
d = d[filtered == TRUE]
d[, datetime_ := as.POSIXct(as.character(datetime_), tz = 'UTC')]
st_transform_DT(d)Subset ID’s with overlapping data
# start and end of the data
d[, ID_year := paste0(ID, '_', year_)]
d[, start := min(datetime_), by = ID_year]
d[, end := max(datetime_), by = ID_year]
dID = unique(d[, .(ID_year, year_, ID, sex, start, end)], by = 'ID_year')
# all pairwise combinations
dpu = CJ(ID_year1 = dID[, ID_year], ID_year2 = dID[, ID_year], unique = TRUE)
# merge with start and end
dpu = merge(dpu, dID[, .(ID_year1 = ID_year, year_1 = year_, ID1 = ID, sex1 = sex, start1 = start, end1 = end)],
by = 'ID_year1')
dpu = merge(dpu, dID[, .(ID_year2 = ID_year, year_2 = year_, ID2 = ID, sex2 = sex, start2 = start, end2 = end)],
by = 'ID_year2')
# overlapping intervals
dpu[, overlap := DescTools::Overlap(c(start1, end1), c(start2, end2)), by = 1:nrow(dpu)]
dpu = dpu[overlap > 0] # exclude non-overlapping data
dpu = dpu[ID_year1 != ID_year2] # exclude within-individual dataDistance to all closest positions
# register cores
registerDoFuture()
plan(multisession)
# loop for all pairwise closest data
dp = foreach(i = 1:nrow(dpu), .combine = 'rbind', .packages = c('data.table')) %dopar% {
# subset pair
dpus = dpu[i, ]
d1 = d[year_ == dpus[1, year_1] & ID == dpus[1, ID1],
.(ID1 = ID, sex1 = sex, datetime_ = datetime_, lat1 = lat, lon1 = lon, gps_speed1 = gps_speed,
altitude1 = altitude)]
d2 = d[year_ == dpus[1, year_2] & ID == dpus[1, ID2],
.(ID2 = ID, sex2 = sex, datetime_ = datetime_, lat2 = lat, lon2 = lon, gps_speed2 = gps_speed,
altitude2 = altitude)]
setkeyv(d1, 'datetime_')
setkeyv(d2, 'datetime_')
# join with closest date
dt = d2[, datetime_2 := (datetime_)][d1, roll = 'nearest']
dt = dt[, .(ID1, ID2, datetime_1 = datetime_, datetime_2)]
dt[, time_btw := abs(as.numeric(difftime(datetime_1, datetime_2, units = 'min')))]
# subset closest position in time
dt[, time_btw_min := min(time_btw), by = .(ID1, datetime_2)]
dt = dt[time_btw == time_btw_min]
dt[, time_btw_min := NULL]
# merge with other data
dt = merge(dt, d1[, .(datetime_1 = datetime_, sex1, lat1, lon1, gps_speed1, altitude1)],
by = 'datetime_1', all.x = TRUE)
dt = merge(dt, d2[, .(datetime_2 = datetime_, sex2, lat2, lon2, gps_speed2, altitude2)],
by = 'datetime_2', all.x = TRUE)
# calculate distance
dt[, distance_pair := sqrt(sum((c(lon1, lat1) - c(lon2, lat2))^2)) , by = 1:nrow(dt)]
dt
}
# remove rows with time difference bigger 10 min
dp = dp[time_btw < 10]
# set colour order
setcolorder(dp, c('ID1', 'ID2', 'sex1', 'sex2', 'datetime_1', 'datetime_2', 'time_btw', 'lat1', 'lon1', 'lat2', 'lon2',
'distance_pair', 'gps_speed1', 'gps_speed2', 'altitude1', 'altitude2'))
# check data
ggplot(data = dp) +
geom_histogram(aes(time_btw))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 5% 50% 95%
## 0.2666667 2.6000000 4.9666667## [1] 0## [1] 15637# save data
# fwrite(dp, './DATA/PRODUCTS/PAIR_WISE_DIST_CLOSEST.txt', quote = TRUE, sep = '\t', row.names = FALSE)
# version information
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8 LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.2 doFuture_0.12.2 future_1.31.0 ggnewscale_0.4.8 windR_1.0.0 stringr_1.5.0 knitr_1.42 sf_1.0-13
## [9] foreach_1.5.2 auksRuak_0.1.0 viridis_0.6.2 viridisLite_0.4.1 anytime_0.3.9 ggplot2_3.4.0 sdb_2019.12 magrittr_2.0.3
## [17] data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] sfext_0.1.0 lubridate_1.9.2 bit64_4.0.5 httr_1.4.4 rprojroot_2.0.3 bslib_0.4.2 tools_4.2.2
## [8] utf8_1.2.3 R6_2.5.1 KernSmooth_2.23-20 rgeos_0.6-1 DBI_1.1.3 colorspace_2.1-0 raster_3.6-14
## [15] withr_2.5.0 sp_1.6-0 tidyselect_1.2.0 gridExtra_2.3 Exact_3.2 RMySQL_0.10.25 bit_4.0.5
## [22] compiler_4.2.2 cli_3.6.0 expm_0.999-7 labeling_0.4.2 sass_0.4.5 scales_1.2.1 classInt_0.4-9
## [29] mvtnorm_1.1-3 proxy_0.4-27 askpass_1.1 digest_0.6.31 rmarkdown_2.20 pkgconfig_2.0.3 htmltools_0.5.4
## [36] parallelly_1.34.0 filenamr_0.1.0 fastmap_1.1.0 rlang_1.0.6 RMariaDB_1.2.2 readxl_1.4.1 rstudioapi_0.14
## [43] farver_2.1.1 jquerylib_0.1.4 generics_0.1.3 jsonlite_1.8.4 dplyr_1.1.0 Matrix_1.5-1 Rcpp_1.0.10
## [50] DescTools_0.99.47 munsell_0.5.0 fansi_1.0.4 lifecycle_1.0.3 terra_1.7-3 stringi_1.7.12 yaml_2.3.7
## [57] snakecase_0.11.0 MASS_7.3-58.1 rootSolve_1.8.2.3 grid_4.2.2 parallel_4.2.2 listenv_0.9.0 lmom_2.9
## [64] lattice_0.20-45 hms_1.1.3 pillar_1.8.1 boot_1.3-28 gld_2.6.6 codetools_0.2-18 glue_1.6.2
## [71] evaluate_0.20 vctrs_0.5.2 cellranger_1.1.0 gtable_0.3.1 purrr_1.0.1 cachem_1.0.6 xfun_0.40
## [78] janitor_2.1.0 e1071_1.7-13 class_7.3-20 ggspatial_1.1.7 tibble_3.1.8 iterators_1.0.14 cliExtras_0.1.0
## [85] units_0.8-2 timechange_0.2.0 globals_0.16.2