From 90eaccaccc4ec8feb3c529168ebb0587b13c0632 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 28 Nov 2023 16:31:41 +0100 Subject: [PATCH 01/14] GIS file loading now handles missing table file #1 and simulation hbGPS output is created to aid testing #1 --- R/hbGIS.R | 95 +++++++++++++++++++++++++++++++++++-------- inst/dev_code_hbGIS.R | 21 +++++++--- 2 files changed, 92 insertions(+), 24 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 1a8b064..ef5c561 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -37,6 +37,15 @@ hbGIS <- function(gisdir = "", Nlocations = length(locationNames) loca = vector("list", Nlocations) names(loca) = locationNames + # Initialise loca object which is a list that for each location + # has a list with four items: + # 1 - location shape + # 2 - location neighbourhood shape + # 3 - shape file to read 1 + # 4 - shape file to read 2 + # Sometimes neighbourhood is referred to as buffer + # For some location types such as parks some of these are missing, + # which for now is ok, those lists will be left empty. for (i in 1:Nlocations) { loca[[i]] = vector("list", 4) names(loca[[i]]) = c(locationNames[i], paste0(locationNames[i], "_nbh"), @@ -48,6 +57,7 @@ hbGIS <- function(gisdir = "", find_file = function(path, namelowercase) { allcsvfiles = dir(path, recursive = TRUE, full.names = TRUE) file_of_interest = allcsvfiles[which(tolower(basename(allcsvfiles)) == namelowercase)] + if (length(file_of_interest) == 0) file_of_interest = NULL return(file_of_interest) } # Find shape files @@ -55,33 +65,37 @@ hbGIS <- function(gisdir = "", findfile3 = find_file(path = gisdir, namelowercase = paste0(locationNames[jj], "_table.shp")) if (!is.null(findfile3)) { loca[[jj]][3] = findfile3 - } else { - stop(paste0("unable to find ", findfile3)) + loca[[jj]][[1]] = sf::read_sf(loca[[jj]][3]) #e.g. school or home } findfile4 = find_file(path = gisdir, namelowercase = paste0("loc_", locationNames[jj], "buffers.shp")) - if (!is.null(findfile3)) { + if (!is.null(findfile4)) { loca[[jj]][4] = findfile4 } else { - stop(paste0("unable to find ", findfile4)) + # If it is not a buffer file then it is space file, eg. public park + findfile4 = find_file(path = gisdir, namelowercase = paste0("loc_", locationNames[jj], ".shp")) + if (!is.null(findfile4)) { + loca[[jj]][4] = findfile4 + } + } + # Only load shape file if it exists + if (!is.null(loca[[i]][4][[1]])) { + loca[[jj]][[2]] = sf::read_sf(loca[[jj]][4]) #e.g. school_nbh, home_nbh, or park } - loca[[jj]][[1]] = sf::read_sf(loca[[jj]][3]) #home_nbh - loca[[jj]][[2]] = sf::read_sf(loca[[jj]][4]) #school_nbh } + # Force id numbers to be character( locationNames = names(loca) for (i in 1:length(loca)) { - if (locationNames[i] != "home") { - loc_id = paste0(groupinglocation, "_id") - } else if (locationNames[i] == "home") { # assumption that home is always the identifier for individuals - loc_id = "identifier" - } - for (j in 1:2) { - loca[[i]][j][[1]][[loc_id]] = as.character(loca[[i]][j][[1]][[loc_id]]) + if (!is.null(loca[[i]][3][[1]])) { # only consider locations with data + for (j in 1:2) { + loc_id = grep(pattern = "ID|identifier|_id", x = colnames(loca[[i]][j][[1]])) + if (length(loc_id) > 0) { + loca[[i]][j][[1]][[loc_id]] = as.character(loca[[i]][j][[1]][[loc_id]]) + } + } } } - - # >>> TO DO: Check that public locations that are not linked to an ID can be loaded in code above <<< - + #=============================================== # hbGPS output (PALMS output) #=============================================== @@ -95,6 +109,47 @@ hbGIS <- function(gisdir = "", palms_country_files <- list.files(path = palmsdir, pattern = "*.csv", full.names = TRUE) # skip the combined file that palms generates palms_country_files = grep(pattern = "combined.csv", x = palms_country_files, invert = TRUE, value = TRUE) + if (length(palms_country_files) == 0) { + # Simulate hbGPS output (only for code developement purposes) + Nmin = 500 + now = Sys.time() + dateTime = seq(now, now + ((Nmin - 1) * 60), by = 60) + example_object = loca[[1]][[2]][1,] + point_in_object = st_sample(x = example_object, size = 1) + xy = sf::st_coordinates(x = point_in_object) + + # latitude is for most of the time 1 lat degree away from location + # but for 30 minutes inside the location surround by 5 minute trips before and after + trip = seq(xy[1] - 1, xy[1] - 0.2, by = 0.2) + away = rep(xy[1] - 1, (Nmin/2) - 20) + lat = c(away, trip, rep(xy[1], 30), rev(trip), away) + lon = rep(xy[2], Nmin) # lon stays the same the entire time + tripNumber = c(rep(0, length(away)), rep(1, 5), rep(0, 30), rep(2, 5), rep(0, length(away))) + sedentaryBoutNumber = c(rep(1, length(away)), rep(0, 5), rep(2, 30), rep(0, 5), rep(3, length(away))) + tripType = rep(0, Nmin) + tripMOT = rep(0, Nmin) + tripType[which(diff(tripNumber) > 0)] = 1 + tripType[which(diff(tripNumber) < 0)] = 4 + tripMOT[which(tripNumber != 0)] = 3 + hbGPSout = data.frame(identifier = "sim1", + dateTime = dateTime, + dow = rep(5, Nmin), + lat = lat, + lon = lon, + fixTypeCode = rep(-1, Nmin), + iov = rep(2, Nmin), # all the time outdoor (indoor, outdoor, vehicle) + tripNumber = tripNumber, + tripType = tripType, + tripMOT = tripMOT, + activity = rep(0, Nmin), + activityIntensity = rep(0, Nmin), + activityBoutNumber = rep(0, Nmin), + sedentaryBoutNumber = sedentaryBoutNumber) + if (!dir.exists(palmsdir)) dir.create(palmsdir, recursive = TRUE) + palms_country_files = paste0(palmsdir, "/combined.csv") + write.csv(hbGPSout, file = palms_country_files, row.names = FALSE) + } + # read and combine palms csv output files csv_palms <- lapply(palms_country_files, FUN = readr::read_csv, col_types = list( identifier = readr::col_character(), @@ -120,16 +175,20 @@ hbGIS <- function(gisdir = "", write.csv(palms_reduced_cleaned, PALMS_reduced_file) palms = palmsplusr::read_palms(PALMS_reduced_file, verbose = FALSE) palms$datetime = as.POSIXct(palms$datetime, format = "%d/%m/%Y %H:%M:%S", tz = "") - + #=============================================== # Load linkage file and identify which PALMS ids and home/school # ids are missing, but allow for publiclocations that are not linked # to an ID #=============================================== + # >>> TO DO: + # >>> Allow for skipping of gislinkfile if it is missing + # >>> Allow for public locations that are not linked to an ID + participant_basis = read_csv(gislinkfile, show_col_types = FALSE) - # >>> TO DO: Allow for public locations that are not linked to an ID <<< + # Check for missing IDs ------------------------------------------------------------------------- withoutMissingId = check_missing_id(participant_basis, palmsplus_folder, dataset_name, palms, diff --git a/inst/dev_code_hbGIS.R b/inst/dev_code_hbGIS.R index b220f4e..df4d48e 100644 --- a/inst/dev_code_hbGIS.R +++ b/inst/dev_code_hbGIS.R @@ -17,13 +17,22 @@ library(purrr) library(geosphere) library(lwgeom) -hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/GIS", - palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/hbGPSoutput", - gislinkfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/Tables/participant_basis.csv", - outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010", - dataset_name = "NBBB2010", - configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") + +hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/loc_greenspace", + palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/hbGPSoutput", + gislinkfile = NULL, + outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023", + dataset_name = "JasperNov2023", + configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") kkk +# +# hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/GIS", +# palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/hbGPSoutput", +# gislinkfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/Tables/participant_basis.csv", +# outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010", +# dataset_name = "NBBB2010", +# configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") +# kkk # Belgium dataset # hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/GIS", From ab807e9e166937cf560db6a60d65e97e162a7dca Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 30 Nov 2023 15:23:35 +0100 Subject: [PATCH 02/14] adapting code to handle locations not linked to id --- R/build_days.R | 19 +++--- R/hbGIS.R | 156 +++++++++++++++++++++++++++++++++---------------- 2 files changed, 116 insertions(+), 59 deletions(-) diff --git a/R/build_days.R b/R/build_days.R index 3250f0a..49d9416 100644 --- a/R/build_days.R +++ b/R/build_days.R @@ -30,10 +30,10 @@ #' @export # Code modified from https://thets.github.io/palmsplusr/ build_days <- function(data = NULL, verbose = TRUE, - palmsplus_domains = NULL, - palmsplus_fields = NULL, - loca = NULL, - participant_basis = NULL) { + palmsplus_domains = NULL, + palmsplus_fields = NULL, + loca = NULL, + participant_basis = NULL) { # Note: # home, school, home_nbh, school_nbh (or similar) need to be present, # because the functions that are passed on assume that they exist @@ -41,8 +41,10 @@ build_days <- function(data = NULL, verbose = TRUE, Nlocations = length(loca) identifier = NULL for (i in 1:Nlocations) { - txt = paste0(names(loca[[i]])[1], " = loca[[i]][[1]]") - eval(parse(text = txt)) + for (j in 1:2) { + txt = paste0(names(loca[[i]])[j], " = loca[[i]][[j]]") + eval(parse(text = txt)) + } } duration = datetime = name = domain_field = NULL @@ -58,21 +60,20 @@ build_days <- function(data = NULL, verbose = TRUE, domain_args <- setNames("1", "total") %>% lapply(parse_expr) domain_args <- c(domain_args, setNames(domain_fields[[2]], domain_fields[[1]]) %>% lapply(parse_expr)) - data <- data %>% mutate(!!! domain_args) %>% mutate_if(is.logical, as.integer) fields <- palmsplus_fields %>% filter(domain_field == TRUE) %>% pull(name) - data <- data %>% + data <- data %>% st_set_geometry(NULL) %>% dplyr::select(identifier, datetime, any_of(domain_names), all_of(fields)) %>% mutate(duration = 1) %>% mutate_at(vars(-identifier,-datetime), ~ . * palms_epoch(data) / 60) %>% group_by(identifier, date = as.Date(datetime)) %>% dplyr::select(-datetime) - + x <- list() for (i in domain_names) { x[[i]] <- data %>% diff --git a/R/hbGIS.R b/R/hbGIS.R index ef5c561..a9a18aa 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -61,24 +61,28 @@ hbGIS <- function(gisdir = "", return(file_of_interest) } # Find shape files + locationNames_nbh = locationNames_table = NULL for (jj in 1:Nlocations) { findfile3 = find_file(path = gisdir, namelowercase = paste0(locationNames[jj], "_table.shp")) if (!is.null(findfile3)) { loca[[jj]][3] = findfile3 loca[[jj]][[1]] = sf::read_sf(loca[[jj]][3]) #e.g. school or home + locationNames_table = c(locationNames_table, locationNames[jj]) } findfile4 = find_file(path = gisdir, namelowercase = paste0("loc_", locationNames[jj], "buffers.shp")) if (!is.null(findfile4)) { loca[[jj]][4] = findfile4 + locationNames_nbh = c(locationNames_nbh, locationNames[jj]) } else { # If it is not a buffer file then it is space file, eg. public park findfile4 = find_file(path = gisdir, namelowercase = paste0("loc_", locationNames[jj], ".shp")) if (!is.null(findfile4)) { loca[[jj]][4] = findfile4 + locationNames_nbh = c(locationNames_nbh, locationNames[jj]) } } # Only load shape file if it exists - if (!is.null(loca[[i]][4][[1]])) { + if (!is.null(loca[[jj]][4][[1]])) { loca[[jj]][[2]] = sf::read_sf(loca[[jj]][4]) #e.g. school_nbh, home_nbh, or park } } @@ -128,7 +132,7 @@ hbGIS <- function(gisdir = "", sedentaryBoutNumber = c(rep(1, length(away)), rep(0, 5), rep(2, 30), rep(0, 5), rep(3, length(away))) tripType = rep(0, Nmin) tripMOT = rep(0, Nmin) - tripType[which(diff(tripNumber) > 0)] = 1 + tripType[which(diff(tripNumber) > 0) + 1] = 1 tripType[which(diff(tripNumber) < 0)] = 4 tripMOT[which(tripNumber != 0)] = 3 hbGPSout = data.frame(identifier = "sim1", @@ -172,34 +176,31 @@ hbGIS <- function(gisdir = "", # Write to csv and read using read_palms to format the object as expected from the rest of the code PALMS_reduced_file = normalizePath(paste0(palmsplus_folder, "/", stringr::str_interp("PALMS_${dataset_name}_reduced.csv"))) # if (verbose) cat(paste0("\nCheck PALMS_reduced_file: ", PALMS_reduced_file)) - write.csv(palms_reduced_cleaned, PALMS_reduced_file) + write.csv(palms_reduced_cleaned, PALMS_reduced_file, row.names = FALSE) palms = palmsplusr::read_palms(PALMS_reduced_file, verbose = FALSE) palms$datetime = as.POSIXct(palms$datetime, format = "%d/%m/%Y %H:%M:%S", tz = "") - #=============================================== # Load linkage file and identify which PALMS ids and home/school # ids are missing, but allow for publiclocations that are not linked # to an ID #=============================================== - # >>> TO DO: - # >>> Allow for skipping of gislinkfile if it is missing - # >>> Allow for public locations that are not linked to an ID - - participant_basis = read_csv(gislinkfile, show_col_types = FALSE) - - - - - # Check for missing IDs ------------------------------------------------------------------------- - withoutMissingId = check_missing_id(participant_basis, palmsplus_folder, dataset_name, palms, - loca, groupinglocation = groupinglocation, - verbose = verbose) - palms = withoutMissingId$palms - participant_basis = withoutMissingId$participant_basis - loca = withoutMissingId$loca - write.csv(participant_basis, paste0(palmsplus_folder, "/", stringr::str_interp("participant_basis_${dataset_name}.csv"))) # store file for logging purposes only - if (length(participant_basis) == 0 || nrow(participant_basis) == 0) { - stop("\nParticipant basis file does not include references for the expected recording IDs") + if (length(gislinkfile) > 0) { + participant_basis = read_csv(gislinkfile, show_col_types = FALSE) + + + # Check for missing IDs ------------------------------------------------------------------------- + withoutMissingId = check_missing_id(participant_basis, palmsplus_folder, dataset_name, palms, + loca, groupinglocation = groupinglocation, + verbose = verbose) + palms = withoutMissingId$palms + participant_basis = withoutMissingId$participant_basis + loca = withoutMissingId$loca + write.csv(participant_basis, paste0(palmsplus_folder, "/", stringr::str_interp("participant_basis_${dataset_name}.csv"))) # store file for logging purposes only + if (length(participant_basis) == 0 || nrow(participant_basis) == 0) { + stop("\nParticipant basis file does not include references for the expected recording IDs") + } + } else { + participant_basis = "" } #=============================================== @@ -234,41 +235,96 @@ hbGIS <- function(gisdir = "", CONF = read.csv(config, sep = ",") CONF$start_criteria = "" CONF$end_criteria = "" - # add standard location based fields to CONF object: + #===================================================== + # Expand CONF with standard location based fields + #===================================================== + # palmsplus_domain: + #------------------- + # ignore stored definition as we no longer use this + + CONF = CONF[which(CONF[,1] != "palmsplus_domain"), ] + element1 = ifelse(length(locationNames_table) > 0, yes = paste0("!", paste0("at_", locationNames_table, collapse = " & !"), " & "), no = "") + CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + "transport", + paste0(element1, "(pedestrian | bicycle | vehicle)"), + TRUE, NA, "", "") + CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + "other", + paste0(element1, "(!pedestrian & !bicycle & !vehicle)", # removed because theorectically possible + ifelse(test = length(locationNames_nbh) > 0, yes = " & ", no = ""), + paste0("!", paste0("at_", locationNames_nbh, "_nbh", collapse = " & !"))), + TRUE, NA, "", "") for (i in 1:Nlocations) { - if (locationNames[i] == "home") { - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + # palmsplus_domain: + #------------------- + if (locationNames[i] %in% locationNames_table) { + CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + locationNames[i], paste0("at_", locationNames[i]), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i],", identifier == i), identifier)"), - NA, "", "", "") - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", - paste0("at_", locationNames[i], "_nbh"), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i], "_nbh, identifier == i), identifier)"), - NA, "", "", "") + TRUE, + NA, "", "") + } else if (locationNames[i] %in% locationNames_nbh) { + # condition that only needs to be used if table element is present + at_table = ifelse(test = locationNames[i] %in% locationNames_table == TRUE, yes = paste0("!at_", locationNames[i], " &"), no = "") + CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + paste0(locationNames[i], "_nbh"), + paste0(at_table, " at_", locationNames[i], "_nbh", + " & (!vehicle)"), # removed !pedestrian & !bicycle & because unclear why these are not possible in a neighbourhood, e.g. park + TRUE, + NA, "", "") + } + # palmsplus_field: + #------------------- + if (!is.null(loca[[i]][[1]])) { + # only do this if there is table data (meaning that location is linked to participant basis file) + if (locationNames[i] == "home") { + CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + paste0("at_", locationNames[i]), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i],", identifier == i), identifier)"), + NA, "", "", "") + CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + paste0("at_", locationNames[i], "_nbh"), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i], "_nbh, identifier == i), identifier)"), + NA, "", "", "") + } else { + CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + paste0("at_", locationNames[i]), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i],",", locationNames[i], + "_id == participant_basis %>% filter(identifier == i) %>% pull(", + locationNames[i], "_id)))"), + NA, "", "", "") + CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + paste0("at_", locationNames[i], "_nbh"), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i], "_nbh,", locationNames[i], + "_id == participant_basis %>% filter(identifier == i) %>% pull(", + locationNames[i], "_id)))"), + NA, "", "", "") + } } else { - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", - paste0("at_", locationNames[i]), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i],",", locationNames[i], - "_id == participant_basis %>% filter(identifier == i) %>% pull(", - locationNames[i], "_id)))"), - NA, "", "", "") + # palmsplus_domain: + #------------------- + # locations not in linkagefile + # Note palms_in_polygon allows for aggregating polygons based by variable CONF[nrow(CONF) + 1, ] = c("palmsplus_field", paste0("at_", locationNames[i], "_nbh"), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i], "_nbh,", locationNames[i], - "_id == participant_basis %>% filter(identifier == i) %>% pull(", - locationNames[i], "_id)))"), + paste0("suppressMessages(any(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), NA, "", "", "") } + reference_location = ifelse("home" %in% locationNames, yes = "home", no = locationNames[1]) for (j in 1:Nlocations) { + #add nbh if table is missing + element2 = ifelse(test = is.null(loca[[j]][[1]]), yes = "_nbh", no = "") + # trajectory_location: + #------------------- CONF[nrow(CONF) + 1, ] = c("trajectory_location", - paste0(locationNames[i], "_", locationNames[i]), - paste0("at_", locationNames[i]), NA, "at_home", - paste0("at_", locationNames[i]), - paste0("at_", locationNames[j])) + paste0(locationNames[i], element2, "_", locationNames[i], element2), + paste0("at_", locationNames[i]), NA, paste0("at_", reference_location, element2), + paste0("at_", locationNames[i], element2), + paste0("at_", locationNames[j], element2)) } CONF = CONF[!duplicated(CONF),] } @@ -310,7 +366,7 @@ hbGIS <- function(gisdir = "", Nlocation_objects = NULL for (i in 1:Nlocations) { - Nlocation_objects = c(Nlocation_objects, length(loca[[i]][[1]]), length(loca[[i]][[2]])) + Nlocation_objects = c(Nlocation_objects, length(loca[[i]][[2]])) # at least a nbh object is expected #length(loca[[i]][[1]]), } if (verbose) cat("\n<<< building palmsplus...\n") if (length(palms) > 0 & length(palmsplus_fields) & From a6f130241b437748313321a4425a0e0034aac8a7 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 30 Nov 2023 15:26:58 +0100 Subject: [PATCH 03/14] tidying up syntax --- R/hbGIS.R | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index a9a18aa..4ee5269 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -27,6 +27,7 @@ hbGIS <- function(gisdir = "", publiclocation = "park" groupinglocation = "school" lon = identifier = palms = NULL # . = was also included, but probably wrong + #=============================================== # GIS files #=============================================== @@ -179,6 +180,7 @@ hbGIS <- function(gisdir = "", write.csv(palms_reduced_cleaned, PALMS_reduced_file, row.names = FALSE) palms = palmsplusr::read_palms(PALMS_reduced_file, verbose = FALSE) palms$datetime = as.POSIXct(palms$datetime, format = "%d/%m/%Y %H:%M:%S", tz = "") + #=============================================== # Load linkage file and identify which PALMS ids and home/school # ids are missing, but allow for publiclocations that are not linked @@ -228,9 +230,6 @@ hbGIS <- function(gisdir = "", config <- system.file("testfiles_hbGIS/config_hbGIS.csv", package = "hbGIS")[1] } - - # >>> TO DO: Add fields for public places <<< - # adding fields CONF = read.csv(config, sep = ",") CONF$start_criteria = "" @@ -254,6 +253,7 @@ hbGIS <- function(gisdir = "", ifelse(test = length(locationNames_nbh) > 0, yes = " & ", no = ""), paste0("!", paste0("at_", locationNames_nbh, "_nbh", collapse = " & !"))), TRUE, NA, "", "") + for (i in 1:Nlocations) { # palmsplus_domain: #------------------- @@ -273,6 +273,7 @@ hbGIS <- function(gisdir = "", TRUE, NA, "", "") } + # palmsplus_field: #------------------- if (!is.null(loca[[i]][[1]])) { @@ -328,6 +329,7 @@ hbGIS <- function(gisdir = "", } CONF = CONF[!duplicated(CONF),] } + palmsplusr_field_rows = which(CONF$context == "palmsplus_field") palmsplus_fields = tibble(name = CONF$name[palmsplusr_field_rows], formula = CONF$formula[palmsplusr_field_rows], @@ -337,24 +339,24 @@ hbGIS <- function(gisdir = "", palmsplus_domains = tibble(name = CONF$name[palmsplusr_domain_rows], formula = CONF$formula[palmsplusr_domain_rows], domain_field = CONF$domain_field[palmsplusr_domain_rows]) - # #============================= - # # trajectory_fields + #============================= + # trajectory_fields trajectory_field_rows = which(CONF$context == "trajectory_field") trajectory_fields = tibble(name = CONF$name[trajectory_field_rows], formula = CONF$formula[trajectory_field_rows], after_conversion = CONF$after_conversion[trajectory_field_rows]) - # #============================= - # # multimodal_fields + #============================= + # multimodal_fields multimodal_fields_rows = which(CONF$context == "multimodal_field") multimodal_fields = tibble(name = CONF$name[multimodal_fields_rows], formula = CONF$formula[multimodal_fields_rows]) - # #============================= - # # trajectory locations + #============================= + # trajectory locations trajectory_location_rows = which(CONF$context == "trajectory_location") trajectory_locations = tibble(name = CONF$name[trajectory_location_rows], start_criteria = CONF$start_criteria[trajectory_location_rows], end_criteria = CONF$end_criteria[trajectory_location_rows]) - # save(palms, loca, participant_basis, file = "~/projects/fontys/state_1_gui.RData") + # Run palmsplusr ---------------------------------------------------------- fns = c(paste0(palmsplus_folder, "/", dataset_name, "_palmsplus.csv"), paste0(palmsplus_folder, "/", dataset_name, "_days.csv"), @@ -368,6 +370,7 @@ hbGIS <- function(gisdir = "", for (i in 1:Nlocations) { Nlocation_objects = c(Nlocation_objects, length(loca[[i]][[2]])) # at least a nbh object is expected #length(loca[[i]][[1]]), } + if (verbose) cat("\n<<< building palmsplus...\n") if (length(palms) > 0 & length(palmsplus_fields) & all(Nlocation_objects > 0) & length(participant_basis) > 0) { @@ -382,6 +385,7 @@ hbGIS <- function(gisdir = "", } else { if (verbose) cat("skipped because insufficient input data>>>\n") } + if (verbose) cat("\n<<< building days...") if (length(palmsplus) > 0 & length(palmsplus_domains) > 0 & length(palmsplus_fields) & all(Nlocation_objects > 0) & length(participant_basis) > 0) { @@ -398,13 +402,11 @@ hbGIS <- function(gisdir = "", } else { if (verbose) cat(paste0(" WARNING: no days object produced.")) } - - # sf::st_write(palmsplus, dsn = paste0(palmsplus_folder, "/", dataset_name, "_palmsplus.shp"), append = FALSE) - } else { if (verbose) cat("skipped because insufficient input data>>>\n") } if (verbose) cat(">>>\n") + trajectory_locations = trajectory_locations[order(trajectory_locations$name),] if (verbose) cat("\n<<< building trajectories...\n") if (length(palmsplus) > 0 & length(trajectory_fields) > 0) { @@ -426,6 +428,7 @@ hbGIS <- function(gisdir = "", } else { if (verbose) cat("skipped because insufficient input data>>>\n") } + if (verbose) cat("\n<<< building multimodal...\n") if (length(palmsplus) > 0 & length(multimodal_fields) > 0 & length(trajectory_locations) > 0) { multimodal <- build_multimodal(data = trajectories, From 1e5eba717e8e9b661305005adbe3af6a0ecfd904 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 30 Nov 2023 16:46:47 +0100 Subject: [PATCH 04/14] fix dummy data, fix st_contains call --- R/hbGIS.R | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 4ee5269..166ec83 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -117,7 +117,7 @@ hbGIS <- function(gisdir = "", if (length(palms_country_files) == 0) { # Simulate hbGPS output (only for code developement purposes) Nmin = 500 - now = Sys.time() + now = as.POSIXct("2023-11-30 10:00:00 CET") dateTime = seq(now, now + ((Nmin - 1) * 60), by = 60) example_object = loca[[1]][[2]][1,] point_in_object = st_sample(x = example_object, size = 1) @@ -127,8 +127,8 @@ hbGIS <- function(gisdir = "", # but for 30 minutes inside the location surround by 5 minute trips before and after trip = seq(xy[1] - 1, xy[1] - 0.2, by = 0.2) away = rep(xy[1] - 1, (Nmin/2) - 20) - lat = c(away, trip, rep(xy[1], 30), rev(trip), away) - lon = rep(xy[2], Nmin) # lon stays the same the entire time + lon = c(away, trip, rep(xy[1], 30), rev(trip), away) + lat = rep(xy[2], Nmin) # lon stays the same the entire time tripNumber = c(rep(0, length(away)), rep(1, 5), rep(0, 30), rep(2, 5), rep(0, length(away))) sedentaryBoutNumber = c(rep(1, length(away)), rep(0, 5), rep(2, 30), rep(0, 5), rep(3, length(away))) tripType = rep(0, Nmin) @@ -312,7 +312,7 @@ hbGIS <- function(gisdir = "", # Note palms_in_polygon allows for aggregating polygons based by variable CONF[nrow(CONF) + 1, ] = c("palmsplus_field", paste0("at_", locationNames[i], "_nbh"), - paste0("suppressMessages(any(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), + paste0("suppressMessages(colSums(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), NA, "", "", "") } reference_location = ifelse("home" %in% locationNames, yes = "home", no = locationNames[1]) @@ -370,7 +370,6 @@ hbGIS <- function(gisdir = "", for (i in 1:Nlocations) { Nlocation_objects = c(Nlocation_objects, length(loca[[i]][[2]])) # at least a nbh object is expected #length(loca[[i]][[1]]), } - if (verbose) cat("\n<<< building palmsplus...\n") if (length(palms) > 0 & length(palmsplus_fields) & all(Nlocation_objects > 0) & length(participant_basis) > 0) { @@ -385,7 +384,6 @@ hbGIS <- function(gisdir = "", } else { if (verbose) cat("skipped because insufficient input data>>>\n") } - if (verbose) cat("\n<<< building days...") if (length(palmsplus) > 0 & length(palmsplus_domains) > 0 & length(palmsplus_fields) & all(Nlocation_objects > 0) & length(participant_basis) > 0) { From b4f15ffb915ddc4da4e87bf94b76ff25c581dfcc Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 30 Nov 2023 16:47:07 +0100 Subject: [PATCH 05/14] edits to development scripts --- inst/dev_code_hbGIS.R | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/inst/dev_code_hbGIS.R b/inst/dev_code_hbGIS.R index df4d48e..ead1284 100644 --- a/inst/dev_code_hbGIS.R +++ b/inst/dev_code_hbGIS.R @@ -18,30 +18,30 @@ library(geosphere) library(lwgeom) -hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/loc_greenspace", - palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/hbGPSoutput", - gislinkfile = NULL, - outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023", - dataset_name = "JasperNov2023", - configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") -kkk -# -# hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/GIS", -# palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/hbGPSoutput", -# gislinkfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/Tables/participant_basis.csv", -# outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010", -# dataset_name = "NBBB2010", -# configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") +# hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/loc_greenspace", +# palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023/hbGPSoutput", +# gislinkfile = NULL, +# outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/JasperNov2023", +# dataset_name = "JasperNov2023", +# configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") # kkk +# +hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/GIS", + palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/hbGPSoutput", + gislinkfile = NULL, # "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/Tables/participant_basis.csv", + outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010", + dataset_name = "NBBB2010", + configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/NBBB2010/config_palmsplusr.csv") +kkk # Belgium dataset -# hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/GIS", -# palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/hbGPSoutput", -# gislinkfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/Tables/participant_basis.csv", -# outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata", -# dataset_name = "BEtestdata", -# configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/config_palmsplusr.csv") - +hbGIS(gisdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/GIS/", + palmsdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/hbGPSoutput", + gislinkfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/Tables/participant_basis.csv", + outputdir = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata", + dataset_name = "BEtestdata", + configfile = "D:/Dropbox/Work/sharedfolder/DATA/Habitus/GPSprocessing/BEtestdata/config_palmsplusr.csv") +jkjj # hbGIS(gisdir = "/media/vincent/projects/Habitus/palmsplusr/testdata/GIS", # palmsdir = "/media/vincent/projects/Habitus/palmsplusr/testdata/PALMS_output/", # gislinkfile = "/media/vincent/projects/Habitus/palmsplusr/testdata/Tables/participant_basis.csv", From c81d9eccc8119220c0d0ac17f924a40c72a146ee Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 4 Dec 2023 16:43:03 +0100 Subject: [PATCH 06/14] fix trajectory locations, code did not correctly loop over combinations of locations --- R/hbGIS.R | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 166ec83..7236d11 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -242,14 +242,14 @@ hbGIS <- function(gisdir = "", # ignore stored definition as we no longer use this CONF = CONF[which(CONF[,1] != "palmsplus_domain"), ] - element1 = ifelse(length(locationNames_table) > 0, yes = paste0("!", paste0("at_", locationNames_table, collapse = " & !"), " & "), no = "") + element3 = ifelse(length(locationNames_table) > 0, yes = paste0("!", paste0("at_", locationNames_table, collapse = " & !"), " & "), no = "") CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", "transport", - paste0(element1, "(pedestrian | bicycle | vehicle)"), + paste0(element3, "(pedestrian | bicycle | vehicle)"), TRUE, NA, "", "") CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", "other", - paste0(element1, "(!pedestrian & !bicycle & !vehicle)", # removed because theorectically possible + paste0(element3, "(!pedestrian & !bicycle & !vehicle)", # removed because theorectically possible ifelse(test = length(locationNames_nbh) > 0, yes = " & ", no = ""), paste0("!", paste0("at_", locationNames_nbh, "_nbh", collapse = " & !"))), TRUE, NA, "", "") @@ -318,18 +318,20 @@ hbGIS <- function(gisdir = "", reference_location = ifelse("home" %in% locationNames, yes = "home", no = locationNames[1]) for (j in 1:Nlocations) { #add nbh if table is missing + element1 = ifelse(test = is.null(loca[[i]][[1]]), yes = "_nbh", no = "") element2 = ifelse(test = is.null(loca[[j]][[1]]), yes = "_nbh", no = "") + element4 = ifelse(test = is.null(loca[[reference_location]][[1]]), yes = "_nbh", no = "") # trajectory_location: #------------------- CONF[nrow(CONF) + 1, ] = c("trajectory_location", - paste0(locationNames[i], element2, "_", locationNames[i], element2), - paste0("at_", locationNames[i]), NA, paste0("at_", reference_location, element2), - paste0("at_", locationNames[i], element2), + paste0(locationNames[i], element1, "_", locationNames[j], element2), + paste0("at_", locationNames[i]), NA, + paste0("at_", reference_location, element4), + paste0("at_", locationNames[i], element1), paste0("at_", locationNames[j], element2)) } CONF = CONF[!duplicated(CONF),] } - palmsplusr_field_rows = which(CONF$context == "palmsplus_field") palmsplus_fields = tibble(name = CONF$name[palmsplusr_field_rows], formula = CONF$formula[palmsplusr_field_rows], @@ -408,15 +410,13 @@ hbGIS <- function(gisdir = "", trajectory_locations = trajectory_locations[order(trajectory_locations$name),] if (verbose) cat("\n<<< building trajectories...\n") if (length(palmsplus) > 0 & length(trajectory_fields) > 0) { - trajectories <- build_trajectories(data = palmsplus, trajectory_fields = trajectory_fields, trajectory_locations = trajectory_locations) if (length(trajectories) > 0) { write_csv(trajectories, file = fns[3]) - shp_file = paste0(palmsplus_folder, "/", dataset_name, "_trajecories.shp") + shp_file = paste0(palmsplus_folder, "/", dataset_name, "_trajectories.shp") if (file.exists(shp_file)) file.remove(shp_file) # remove because st_write does not know how to overwrite - sf::st_write(obj = trajectories, dsn = shp_file) if (verbose) cat(paste0(" N rows in trajectories object: ", nrow(trajectories))) } else { @@ -426,7 +426,6 @@ hbGIS <- function(gisdir = "", } else { if (verbose) cat("skipped because insufficient input data>>>\n") } - if (verbose) cat("\n<<< building multimodal...\n") if (length(palmsplus) > 0 & length(multimodal_fields) > 0 & length(trajectory_locations) > 0) { multimodal <- build_multimodal(data = trajectories, From ae95c88fc0da07c005bb0066c4d8754cb865be90 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 12:42:52 +0100 Subject: [PATCH 07/14] replace deprecate dependencies --- NAMESPACE | 1 - R/build_days.R | 6 +++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index dc7df75..78f74cf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,7 +19,6 @@ importFrom(geosphere,distGeo) importFrom(purrr,reduce) importFrom(readr,read_csv) importFrom(readr,write_csv) -importFrom(rlang,UQ) importFrom(rlang,parse_expr) importFrom(stats,as.formula) importFrom(stats,end) diff --git a/R/build_days.R b/R/build_days.R index 49d9416..3af30e2 100644 --- a/R/build_days.R +++ b/R/build_days.R @@ -24,7 +24,7 @@ #' #' @import dplyr #' @import sf -#' @importFrom rlang parse_expr UQ +#' @importFrom rlang parse_expr #' @importFrom purrr reduce #' #' @export @@ -77,8 +77,8 @@ build_days <- function(data = NULL, verbose = TRUE, x <- list() for (i in domain_names) { x[[i]] <- data %>% - filter(UQ(as.name(i)) > 0) %>% - dplyr::select(-one_of(domain_names), duration) %>% + filter(!!(as.name(i)) > 0) %>% + dplyr::select(-any_of(domain_names), duration) %>% summarise_all(~ sum(.)) %>% ungroup() %>% rename_at(vars(-identifier, -date), ~ paste0(i, "_", .)) From a7191ca071dcfaf5e863c9b80dd60dd84f98cc98 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 13:20:17 +0100 Subject: [PATCH 08/14] add count of visits to a domain per day --- R/build_days.R | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/R/build_days.R b/R/build_days.R index 3af30e2..29c96a5 100644 --- a/R/build_days.R +++ b/R/build_days.R @@ -86,5 +86,24 @@ build_days <- function(data = NULL, verbose = TRUE, result <- x %>% reduce(left_join, by = c("identifier" = "identifier", "date" = "date")) + + # Count the number of segments per domain per day per identifier + segmentcount = function(x) { + x = as.numeric(unlist(x)) # x is a tibble column, so first convert to numeric vector + return(length(which(rle(x)$values != 0))) + } + + for (dom in domain_names) { + if (dom != "total") { + result[, dom] <- NA + for (id in unique(result$identifier)) { + for (date in as.Date(unique(result$date[which(result$identifier == id)]))) { + CNT = segmentcount(data[which(data$identifier == id & data$date == date), dom]) + result[which(result$identifier == id & result$date == date), dom] = CNT + } + } + names(result)[which(names(result) == dom)] = paste0(dom, "_segmentcount") + } + } return(result) } From e5f42d01adf13410ff1469149a3261f7fdb14b7e Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 17:48:53 +0100 Subject: [PATCH 09/14] handle sublocations in GIS files, these are now included in day and palmsplus output file, but not in trajectory and multimodal file #1 --- R/hbGIS.R | 129 +++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 94 insertions(+), 35 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 7236d11..23f8bf2 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -24,8 +24,9 @@ hbGIS <- function(gisdir = "", dataset_name = "", configfile = "", verbose = TRUE) { - publiclocation = "park" groupinglocation = "school" + writeshp = FALSE + splitGIS = TRUE lon = identifier = palms = NULL # . = was also included, but probably wrong #=============================================== @@ -56,20 +57,25 @@ hbGIS <- function(gisdir = "", # Helper function to find shape files find_file = function(path, namelowercase) { - allcsvfiles = dir(path, recursive = TRUE, full.names = TRUE) - file_of_interest = allcsvfiles[which(tolower(basename(allcsvfiles)) == namelowercase)] + allshpfiles = dir(path, recursive = TRUE, full.names = TRUE, pattern = "shp") + file_of_interest = allshpfiles[which(tolower(basename(allshpfiles)) == namelowercase)] if (length(file_of_interest) == 0) file_of_interest = NULL return(file_of_interest) } # Find shape files locationNames_nbh = locationNames_table = NULL + ignore4trajectories = NULL for (jj in 1:Nlocations) { + # files with _table at the end of their name indicate that it is + # a location for which entries exist in the GIS-Person linkage file findfile3 = find_file(path = gisdir, namelowercase = paste0(locationNames[jj], "_table.shp")) if (!is.null(findfile3)) { loca[[jj]][3] = findfile3 loca[[jj]][[1]] = sf::read_sf(loca[[jj]][3]) #e.g. school or home locationNames_table = c(locationNames_table, locationNames[jj]) } + # files with "loc_" and "buffers" in their name are assumed to be the areas surrounding + # the table files (if exist) findfile4 = find_file(path = gisdir, namelowercase = paste0("loc_", locationNames[jj], "buffers.shp")) if (!is.null(findfile4)) { loca[[jj]][4] = findfile4 @@ -84,12 +90,43 @@ hbGIS <- function(gisdir = "", } # Only load shape file if it exists if (!is.null(loca[[jj]][4][[1]])) { - loca[[jj]][[2]] = sf::read_sf(loca[[jj]][4]) #e.g. school_nbh, home_nbh, or park + shp_dat = sf::read_sf(loca[[jj]][4]) #e.g. school_nbh, home_nbh, or park + # Check whether there are multiple polygons in the shapefile: + nshp = nrow(shp_dat) + loca[[jj]][[2]] = shp_dat # look at all sublocations combined either way + if (nshp > 1 & splitGIS == TRUE) { + collect_na = NULL + # Treat each polygon as a separate location + for (gi in 1:nshp) { + fn_4 = as.character(unlist(loca[[jj]][4])) + gi2 = Nlocations + gi + loca[[gi2]] = vector("list", 4) + objectname = as.character(st_drop_geometry(shp_dat[gi, "OBJECTID"])) + if (objectname == "NA") collect_na = c(collect_na, gi) + loca[[gi2]][[2]] = shp_dat[gi, ] + loca[[gi2]][[4]] = fn_4 + newname = paste0(locationNames[jj], objectname) + names(loca[[gi2]]) = c(newname, + paste0(newname, "_nbh"), + paste0(newname, "_tablefile"), + paste0(newname, "_locbufferfile")) + locationNames = c(locationNames, newname) + names(loca)[[gi2]] = newname + locationNames_nbh = c(locationNames_nbh, newname) + ignore4trajectories = c(ignore4trajectories, newname) + } + if (length(collect_na) > 0) { + shp_dat = shp_dat[-collect_na, ] + locationNames = locationNames[-collect_na] + loca = loca[-collect_na] + } + } } } + Nlocations = length(locationNames) - # Force id numbers to be character( locationNames = names(loca) + # Force id numbers to be character to standardise format: for (i in 1:length(loca)) { if (!is.null(loca[[i]][3][[1]])) { # only consider locations with data for (j in 1:2) { @@ -100,7 +137,6 @@ hbGIS <- function(gisdir = "", } } } - #=============================================== # hbGPS output (PALMS output) #=============================================== @@ -180,7 +216,6 @@ hbGIS <- function(gisdir = "", write.csv(palms_reduced_cleaned, PALMS_reduced_file, row.names = FALSE) palms = palmsplusr::read_palms(PALMS_reduced_file, verbose = FALSE) palms$datetime = as.POSIXct(palms$datetime, format = "%d/%m/%Y %H:%M:%S", tz = "") - #=============================================== # Load linkage file and identify which PALMS ids and home/school # ids are missing, but allow for publiclocations that are not linked @@ -237,6 +272,8 @@ hbGIS <- function(gisdir = "", #===================================================== # Expand CONF with standard location based fields #===================================================== + if (verbose) cat("\n<<< expand CONF...\n") + # palmsplus_domain: #------------------- # ignore stored definition as we no longer use this @@ -254,11 +291,15 @@ hbGIS <- function(gisdir = "", paste0("!", paste0("at_", locationNames_nbh, "_nbh", collapse = " & !"))), TRUE, NA, "", "") + cnt = nrow(CONF) + CONF[cnt + Nlocations,] = NA + cnt = cnt + 1 + for (i in 1:Nlocations) { # palmsplus_domain: #------------------- if (locationNames[i] %in% locationNames_table) { - CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + CONF[cnt, ] = c("palmsplus_domain", locationNames[i], paste0("at_", locationNames[i]), TRUE, @@ -266,72 +307,80 @@ hbGIS <- function(gisdir = "", } else if (locationNames[i] %in% locationNames_nbh) { # condition that only needs to be used if table element is present at_table = ifelse(test = locationNames[i] %in% locationNames_table == TRUE, yes = paste0("!at_", locationNames[i], " &"), no = "") - CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", + CONF[cnt, ] = c("palmsplus_domain", paste0(locationNames[i], "_nbh"), paste0(at_table, " at_", locationNames[i], "_nbh", " & (!vehicle)"), # removed !pedestrian & !bicycle & because unclear why these are not possible in a neighbourhood, e.g. park TRUE, NA, "", "") } - + cnt = cnt + 1 # palmsplus_field: #------------------- if (!is.null(loca[[i]][[1]])) { # only do this if there is table data (meaning that location is linked to participant basis file) if (locationNames[i] == "home") { - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + CONF[cnt, ] = c("palmsplus_field", paste0("at_", locationNames[i]), paste0("palms_in_polygon(datai, polygons = dplyr::filter(", locationNames[i],", identifier == i), identifier)"), NA, "", "", "") - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + cnt = cnt + 1 + CONF[cnt, ] = c("palmsplus_field", paste0("at_", locationNames[i], "_nbh"), paste0("palms_in_polygon(datai, polygons = dplyr::filter(", locationNames[i], "_nbh, identifier == i), identifier)"), NA, "", "", "") + cnt = cnt + 1 } else { - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + CONF[cnt, ] = c("palmsplus_field", paste0("at_", locationNames[i]), paste0("palms_in_polygon(datai, polygons = dplyr::filter(", locationNames[i],",", locationNames[i], "_id == participant_basis %>% filter(identifier == i) %>% pull(", locationNames[i], "_id)))"), NA, "", "", "") - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + cnt = cnt + 1 + CONF[cnt, ] = c("palmsplus_field", paste0("at_", locationNames[i], "_nbh"), paste0("palms_in_polygon(datai, polygons = dplyr::filter(", locationNames[i], "_nbh,", locationNames[i], "_id == participant_basis %>% filter(identifier == i) %>% pull(", locationNames[i], "_id)))"), NA, "", "", "") + cnt = cnt + 1 } } else { - # palmsplus_domain: - #------------------- # locations not in linkagefile - # Note palms_in_polygon allows for aggregating polygons based by variable - CONF[nrow(CONF) + 1, ] = c("palmsplus_field", + # note that colSums will ensure that sublocation are combined + CONF[cnt, ] = c("palmsplus_field", paste0("at_", locationNames[i], "_nbh"), paste0("suppressMessages(colSums(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), NA, "", "", "") + cnt = cnt + 1 } reference_location = ifelse("home" %in% locationNames, yes = "home", no = locationNames[1]) for (j in 1:Nlocations) { - #add nbh if table is missing - element1 = ifelse(test = is.null(loca[[i]][[1]]), yes = "_nbh", no = "") - element2 = ifelse(test = is.null(loca[[j]][[1]]), yes = "_nbh", no = "") - element4 = ifelse(test = is.null(loca[[reference_location]][[1]]), yes = "_nbh", no = "") - # trajectory_location: - #------------------- - CONF[nrow(CONF) + 1, ] = c("trajectory_location", - paste0(locationNames[i], element1, "_", locationNames[j], element2), - paste0("at_", locationNames[i]), NA, - paste0("at_", reference_location, element4), - paste0("at_", locationNames[i], element1), - paste0("at_", locationNames[j], element2)) + # Do not derive trajectories for sublocation combinations + if (all(locationNames[c(i, j)] %in% ignore4trajectories == FALSE)) { + #add nbh if table is missing + element1 = ifelse(test = is.null(loca[[i]][[1]]), yes = "_nbh", no = "") + element2 = ifelse(test = is.null(loca[[j]][[1]]), yes = "_nbh", no = "") + element4 = ifelse(test = is.null(loca[[reference_location]][[1]]), yes = "_nbh", no = "") + # trajectory_location: + #------------------- + CONF[cnt, ] = c("trajectory_location", + paste0(locationNames[i], element1, "_", locationNames[j], element2), + paste0("at_", locationNames[i]), NA, + paste0("at_", reference_location, element4), + paste0("at_", locationNames[i], element1), + paste0("at_", locationNames[j], element2)) + cnt = cnt + 1 + } } CONF = CONF[!duplicated(CONF),] } + if (verbose) cat(">>>\n\n") palmsplusr_field_rows = which(CONF$context == "palmsplus_field") palmsplus_fields = tibble(name = CONF$name[palmsplusr_field_rows], formula = CONF$formula[palmsplusr_field_rows], @@ -367,7 +416,6 @@ hbGIS <- function(gisdir = "", for (fn in fns) { if (file.exists(fn)) file.remove(fn) } - Nlocation_objects = NULL for (i in 1:Nlocations) { Nlocation_objects = c(Nlocation_objects, length(loca[[i]][[2]])) # at least a nbh object is expected #length(loca[[i]][[1]]), @@ -413,11 +461,20 @@ hbGIS <- function(gisdir = "", trajectories <- build_trajectories(data = palmsplus, trajectory_fields = trajectory_fields, trajectory_locations = trajectory_locations) + + + # library(mapview) + # mapview(list(trajectories, loca$green$green_nbh), col.regions=list("red","blue"), col = list("red", "blue"), + # layer.name = c("trajectories", "green spaces")) + + # browser() if (length(trajectories) > 0) { write_csv(trajectories, file = fns[3]) shp_file = paste0(palmsplus_folder, "/", dataset_name, "_trajectories.shp") - if (file.exists(shp_file)) file.remove(shp_file) # remove because st_write does not know how to overwrite - sf::st_write(obj = trajectories, dsn = shp_file) + if (writeshp == TRUE) { + if (file.exists(shp_file)) file.remove(shp_file) # remove because st_write does not know how to overwrite + sf::st_write(obj = trajectories, dsn = shp_file) + } if (verbose) cat(paste0(" N rows in trajectories object: ", nrow(trajectories))) } else { if (verbose) cat(paste0(" WARNING: no trajectories object produced.")) @@ -439,8 +496,10 @@ hbGIS <- function(gisdir = "", if (length(multimodal) > 0) { write_csv(multimodal, file = fns[4]) shp_file = paste0(palmsplus_folder, "/", dataset_name, "_multimodal.shp") - if (file.exists(shp_file)) file.remove(shp_file) # remove because st_write does not know how to overwrite - sf::st_write(obj = multimodal, dsn = shp_file) + if (writeshp == TRUE) { + if (file.exists(shp_file)) file.remove(shp_file) # remove because st_write does not know how to overwrite + sf::st_write(obj = multimodal, dsn = shp_file) + } if (verbose) cat(paste0(" N rows in multimodal object: ", nrow(multimodal))) } else { if (verbose) cat(paste0(" WARNING: no multimodal object produced.")) From 341584a8d782b17ee299b8e1129ec72e23ebca69 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 18:04:40 +0100 Subject: [PATCH 10/14] replace hyphen trip number values as it is interpretted by Excel as a date --- R/build_multimodal.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/build_multimodal.R b/R/build_multimodal.R index 17566eb..cfd6df8 100644 --- a/R/build_multimodal.R +++ b/R/build_multimodal.R @@ -152,9 +152,9 @@ build_multimodal <- function(data = NULL, group_by(identifier, mmt_number) %>% summarise(start_trip = first(tripnumber), end_trip = last(tripnumber), - trip_numbers = paste0(tripnumber, collapse = "-"), + trip_numbers = paste0(tripnumber, collapse = ">"), n_segments = n(), - mot_order = paste0(mot, collapse = "-"), + mot_order = paste0(mot, collapse = ">"), start = first(start), end = last(end), do_union = FALSE, .groups = 'keep') %>% From d758827e1275ff7ab1abdb4310852dda7c77265c Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 18:05:05 +0100 Subject: [PATCH 11/14] move candidate new arguments to top of function #1 --- R/hbGIS.R | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 23f8bf2..1c92e6a 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -24,9 +24,14 @@ hbGIS <- function(gisdir = "", dataset_name = "", configfile = "", verbose = TRUE) { + # Hard code arguments that may need to become function arguments to give control to user groupinglocation = "school" - writeshp = FALSE - splitGIS = TRUE + writeshp = FALSE # whether to wrist a shape file + splitGIS = TRUE # whether to split sublocations (TRUE) or union them + sublocid = "OBJECTID" # column name in GIS file to identify sublocation id + + + #------------------------------------------------------------ lon = identifier = palms = NULL # . = was also included, but probably wrong #=============================================== @@ -101,7 +106,7 @@ hbGIS <- function(gisdir = "", fn_4 = as.character(unlist(loca[[jj]][4])) gi2 = Nlocations + gi loca[[gi2]] = vector("list", 4) - objectname = as.character(st_drop_geometry(shp_dat[gi, "OBJECTID"])) + objectname = as.character(st_drop_geometry(shp_dat[gi, sublocid])) if (objectname == "NA") collect_na = c(collect_na, gi) loca[[gi2]][[2]] = shp_dat[gi, ] loca[[gi2]][[4]] = fn_4 @@ -492,7 +497,6 @@ hbGIS <- function(gisdir = "", multimodal_fields = multimodal_fields, trajectory_locations = trajectory_locations, verbose = verbose) - if (length(multimodal) > 0) { write_csv(multimodal, file = fns[4]) shp_file = paste0(palmsplus_folder, "/", dataset_name, "_multimodal.shp") From 39136341d9ab2c6a5a7c69adaf921e239b2f28ec Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 18:27:22 +0100 Subject: [PATCH 12/14] Update r_new.yml --- .github/workflows/r_new.yml | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/.github/workflows/r_new.yml b/.github/workflows/r_new.yml index 5c7a488..63a8f4a 100644 --- a/.github/workflows/r_new.yml +++ b/.github/workflows/r_new.yml @@ -6,7 +6,7 @@ on: pull_request: branches: [main, master] -name: R-CMD-check-new +name: R-CMD-check jobs: R-CMD-check: @@ -18,7 +18,7 @@ jobs: fail-fast: false matrix: config: - - {os: macOS-latest, r: 'release'} + - {os: macos-latest, r: 'release'} - {os: windows-latest, r: 'release'} - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - {os: ubuntu-latest, r: 'release'} @@ -29,10 +29,9 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v2 - + - uses: actions/checkout@v3 + - uses: r-lib/actions/setup-tinytex@v2 - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} @@ -41,7 +40,7 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::rcmdcheck, any::XML + extra-packages: any::rcmdcheck needs: check - uses: r-lib/actions/check-r-package@v2 From 1349b76aafda8c0b1b8bd4c16832c1b81f6cd648 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 19 Dec 2023 18:27:29 +0100 Subject: [PATCH 13/14] Update NEWS.Rd --- inst/NEWS.Rd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 881aac6..e1938be 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -1,8 +1,9 @@ \name{NEWS} \title{News for Package \pkg{HabitusGUI}} \newcommand{\cpkg}{\href{http://CRAN.R-project.org/package=#1}{\pkg{#1}}} -\section{Changes in version 0.0.1 (GitHub-only-release date: 09-11-2023)}{ +\section{Changes in version 0.0.1 (GitHub-only-release date: 20-12-2023)}{ \itemize{ \item First release, based on code from the HabitusGUI and palmsplusr libraries + \item Expansion with ability to handle public locations not tight to an individual } } From 394bfc0e7fb06639b8ee7e7a4406427e47a17db7 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Wed, 20 Dec 2023 08:40:09 +0100 Subject: [PATCH 14/14] fix broken functionality for datasets without public spaces --- R/hbGIS.R | 97 ++++++++++++++++++++++++++++--------------------------- 1 file changed, 50 insertions(+), 47 deletions(-) diff --git a/R/hbGIS.R b/R/hbGIS.R index 1c92e6a..d3d007e 100644 --- a/R/hbGIS.R +++ b/R/hbGIS.R @@ -74,7 +74,8 @@ hbGIS <- function(gisdir = "", # files with _table at the end of their name indicate that it is # a location for which entries exist in the GIS-Person linkage file findfile3 = find_file(path = gisdir, namelowercase = paste0(locationNames[jj], "_table.shp")) - if (!is.null(findfile3)) { + publiclocation = ifelse(test = is.null(findfile3), yes = TRUE, no = FALSE) + if (publiclocation == FALSE) { loca[[jj]][3] = findfile3 loca[[jj]][[1]] = sf::read_sf(loca[[jj]][3]) #e.g. school or home locationNames_table = c(locationNames_table, locationNames[jj]) @@ -99,22 +100,27 @@ hbGIS <- function(gisdir = "", # Check whether there are multiple polygons in the shapefile: nshp = nrow(shp_dat) loca[[jj]][[2]] = shp_dat # look at all sublocations combined either way - if (nshp > 1 & splitGIS == TRUE) { + if (publiclocation == TRUE) { + cat(paste0("\n", basename(as.character(unlist(loca[[jj]][4]))), + " => ", paste0(names(shp_dat), collapse = ", "), " (", nshp, " geoms)")) + } + if (nshp > 1 & splitGIS == TRUE & sublocid %in% names(shp_dat) & publiclocation) { collect_na = NULL # Treat each polygon as a separate location for (gi in 1:nshp) { fn_4 = as.character(unlist(loca[[jj]][4])) gi2 = Nlocations + gi loca[[gi2]] = vector("list", 4) + objectname = as.character(st_drop_geometry(shp_dat[gi, sublocid])) if (objectname == "NA") collect_na = c(collect_na, gi) loca[[gi2]][[2]] = shp_dat[gi, ] loca[[gi2]][[4]] = fn_4 newname = paste0(locationNames[jj], objectname) names(loca[[gi2]]) = c(newname, - paste0(newname, "_nbh"), - paste0(newname, "_tablefile"), - paste0(newname, "_locbufferfile")) + paste0(newname, "_nbh"), + paste0(newname, "_tablefile"), + paste0(newname, "_locbufferfile")) locationNames = c(locationNames, newname) names(loca)[[gi2]] = newname locationNames_nbh = c(locationNames_nbh, newname) @@ -129,7 +135,6 @@ hbGIS <- function(gisdir = "", } } Nlocations = length(locationNames) - locationNames = names(loca) # Force id numbers to be character to standardise format: for (i in 1:length(loca)) { @@ -195,7 +200,7 @@ hbGIS <- function(gisdir = "", palms_country_files = paste0(palmsdir, "/combined.csv") write.csv(hbGPSout, file = palms_country_files, row.names = FALSE) } - + # read and combine palms csv output files csv_palms <- lapply(palms_country_files, FUN = readr::read_csv, col_types = list( identifier = readr::col_character(), @@ -287,7 +292,7 @@ hbGIS <- function(gisdir = "", element3 = ifelse(length(locationNames_table) > 0, yes = paste0("!", paste0("at_", locationNames_table, collapse = " & !"), " & "), no = "") CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", "transport", - paste0(element3, "(pedestrian | bicycle | vehicle)"), + paste0(element3, "(pedestrian | bicycle | vehicle)"), TRUE, NA, "", "") CONF[nrow(CONF) + 1, ] = c("palmsplus_domain", "other", @@ -299,25 +304,25 @@ hbGIS <- function(gisdir = "", cnt = nrow(CONF) CONF[cnt + Nlocations,] = NA cnt = cnt + 1 - + for (i in 1:Nlocations) { # palmsplus_domain: #------------------- if (locationNames[i] %in% locationNames_table) { CONF[cnt, ] = c("palmsplus_domain", - locationNames[i], - paste0("at_", locationNames[i]), - TRUE, - NA, "", "") + locationNames[i], + paste0("at_", locationNames[i]), + TRUE, + NA, "", "") } else if (locationNames[i] %in% locationNames_nbh) { # condition that only needs to be used if table element is present at_table = ifelse(test = locationNames[i] %in% locationNames_table == TRUE, yes = paste0("!at_", locationNames[i], " &"), no = "") CONF[cnt, ] = c("palmsplus_domain", - paste0(locationNames[i], "_nbh"), - paste0(at_table, " at_", locationNames[i], "_nbh", - " & (!vehicle)"), # removed !pedestrian & !bicycle & because unclear why these are not possible in a neighbourhood, e.g. park - TRUE, - NA, "", "") + paste0(locationNames[i], "_nbh"), + paste0(at_table, " at_", locationNames[i], "_nbh", + " & (!vehicle)"), # removed !pedestrian & !bicycle & because unclear why these are not possible in a neighbourhood, e.g. park + TRUE, + NA, "", "") } cnt = cnt + 1 # palmsplus_field: @@ -326,42 +331,42 @@ hbGIS <- function(gisdir = "", # only do this if there is table data (meaning that location is linked to participant basis file) if (locationNames[i] == "home") { CONF[cnt, ] = c("palmsplus_field", - paste0("at_", locationNames[i]), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i],", identifier == i), identifier)"), - NA, "", "", "") + paste0("at_", locationNames[i]), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i],", identifier == i), identifier)"), + NA, "", "", "") cnt = cnt + 1 CONF[cnt, ] = c("palmsplus_field", - paste0("at_", locationNames[i], "_nbh"), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i], "_nbh, identifier == i), identifier)"), - NA, "", "", "") + paste0("at_", locationNames[i], "_nbh"), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i], "_nbh, identifier == i), identifier)"), + NA, "", "", "") cnt = cnt + 1 } else { CONF[cnt, ] = c("palmsplus_field", - paste0("at_", locationNames[i]), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i],",", locationNames[i], - "_id == participant_basis %>% filter(identifier == i) %>% pull(", - locationNames[i], "_id)))"), - NA, "", "", "") + paste0("at_", locationNames[i]), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i],",", locationNames[i], + "_id == participant_basis %>% filter(identifier == i) %>% pull(", + locationNames[i], "_id)))"), + NA, "", "", "") cnt = cnt + 1 CONF[cnt, ] = c("palmsplus_field", - paste0("at_", locationNames[i], "_nbh"), - paste0("palms_in_polygon(datai, polygons = dplyr::filter(", - locationNames[i], "_nbh,", locationNames[i], - "_id == participant_basis %>% filter(identifier == i) %>% pull(", - locationNames[i], "_id)))"), - NA, "", "", "") + paste0("at_", locationNames[i], "_nbh"), + paste0("palms_in_polygon(datai, polygons = dplyr::filter(", + locationNames[i], "_nbh,", locationNames[i], + "_id == participant_basis %>% filter(identifier == i) %>% pull(", + locationNames[i], "_id)))"), + NA, "", "", "") cnt = cnt + 1 } } else { # locations not in linkagefile # note that colSums will ensure that sublocation are combined CONF[cnt, ] = c("palmsplus_field", - paste0("at_", locationNames[i], "_nbh"), - paste0("suppressMessages(colSums(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), - NA, "", "", "") + paste0("at_", locationNames[i], "_nbh"), + paste0("suppressMessages(colSums(st_contains(", locationNames[i], "_nbh, datai, sparse = FALSE)))"), + NA, "", "", "") cnt = cnt + 1 } reference_location = ifelse("home" %in% locationNames, yes = "home", no = locationNames[1]) @@ -448,30 +453,28 @@ hbGIS <- function(gisdir = "", loca = loca, participant_basis = participant_basis, verbose = verbose) - if (length(days) > 0) { - if (verbose) cat(paste0(" N rows in days object: ", nrow(days))) + if (verbose) cat(paste0("\n N rows in days object: ", nrow(days))) write_csv(x = days, file = fns[2]) } else { - if (verbose) cat(paste0(" WARNING: no days object produced.")) + if (verbose) cat(paste0("\n WARNING: no days object produced.")) } } else { if (verbose) cat("skipped because insufficient input data>>>\n") } if (verbose) cat(">>>\n") - trajectory_locations = trajectory_locations[order(trajectory_locations$name),] + if (verbose) cat("\n<<< building trajectories...\n") if (length(palmsplus) > 0 & length(trajectory_fields) > 0) { trajectories <- build_trajectories(data = palmsplus, trajectory_fields = trajectory_fields, trajectory_locations = trajectory_locations) - # library(mapview) # mapview(list(trajectories, loca$green$green_nbh), col.regions=list("red","blue"), col = list("red", "blue"), # layer.name = c("trajectories", "green spaces")) - + # browser() if (length(trajectories) > 0) { write_csv(trajectories, file = fns[3])