Add Admin 1-4 boundary information to a gps point dataset

add-adm-to-gps.R

# load libraries required for analysis
library(tidyverse)
library(sf)
library(qxl)
library(fs)

# path to the GIS Centre Palestine geodatabase file on my computer
# containing boundaries for admin 1-4
path_adm <-  path(
  Sys.getenv("SHAREPOINT_PATH"),
  "OutbreakTools - GeoBase",
  "PSE",
  "resources",
  "shp",
  "GeoMSF",
  "MasterData",
  "MasterData.gdb"
)

# path to excel file with incident gps data on my computer
path_gps <- "~/Downloads/Palestine Incidents GPS Coordinates.xlsx"

# read boundary data for admin 1-4, only keeping name and pcode columns
adm1 <- st_read(path_adm, layer = "main_owner_polbnd_adm1_a_msf") |> 
  select(adm1_name = name_en, adm1_pcode = pcode)

adm2 <- st_read(path_adm, layer = "main_owner_polbnd_adm2_a_msf") |> 
  select(adm2_name = name_en, adm2_pcode = pcode)

adm3 <- st_read(path_adm, layer = "main_owner_polbnd_adm3_a_msf") |> 
  select(adm3_name = name_en, adm3_pcode = pcode)

adm4 <- st_read(path_adm, layer = "main_owner_polbnd_adm_other_a_msf") |> 
  select(adm4_name = name_en, adm4_pcode = pcode)

# read GPS incidents excel file
gps_raw <- qxl::qread(path_gps)

# split lat and lon into 2 separate columns
# convert coordinates to a point geometry type
# spatial join admin boundary information based on what polygon 
# each point falls in for each admin level
# replace missing values with '(Unknown)' - these are points that did
# not fall within any polygons at the given admin level 
gps <- gps_raw |> 
  separate(Coordinates, into = c("lat", "lon"), sep = ", | ") |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |> 
  st_join(adm1) |> 
  st_join(adm2) |> 
  st_join(adm3) |> 
  st_join(adm4) |> 
  mutate(across(starts_with("adm"), ~replace_na(.x, "(Unknown)")))

# create a summary table of incidents by admin level
df_summary <- gps |> 
  st_drop_geometry() |> 
  count(
    adm1_name,
    adm2_name,
    adm3_name,
    adm4_name,
    name = "incidents",
    sort = TRUE
  ) 

# join admin data to original dataset and save as new excel file
df_out <- gps_raw |> 
  left_join(st_drop_geometry(gps), by = "ID")

qxl::qxl(
  list(data = df_out, summary = df_summary),
  file = "~/Downloads/palestine-incidents-gps-with-admin-levels.xlsx"
)