In pratikunterwegs/holeybirds: Link bird moult with movement.

Prepare data for selection analysis

Here we prepare the data for SSFs using the amt package.

Load libraries and prepare files

# libs for data
library(data.table)
library(lubridate)
library(glue)
library(dplyr)
library(tidyr)
library(purrr)

# for vectors
library(sf)

# for rasters
library(terra)

# library for SSF
library(amt)

# read data
data <- fread("data/results/data_patch_summary_ppa.csv")
points <- fread("data/results/data_patch_points_ppa.csv")

data[, time_median := as.POSIXct(
  time_median,
  origin = "1970-01-01", t = "Asia/Jerusalem"
)]

Prepare for SSF

Split by id and date

# make tibble and nest by id and date + other identifiers
setDF(data)

# select data cols for id
data <- select(
  data,
  id, patch, date, x_median, y_median, time_median
)

# nest the tibble,
tracks <- nest(
  data,
  data = -matches(c("id", "date"))
)

# prepare link point data
setDF(points)
points <- nest(
  points,
  data = -matches(c("id", "date"))
)

Make `amt` objects

# make amt objects from the nested data
# keep attribute columns
tracks <- mutate(
  tracks,
  data = map(data, function(df) {
    amt::make_track(
      tbl = df,
      .x = x_median,
      .y = y_median,
      .t = time_median,
      all_cols = T,
      crs = sp::CRS(
        st_crs(2039)$proj4string
      )
    )
  })
)

# count rows in data
tracks <- mutate(
  tracks,
  nrow = map_int(data, nrow)
)

# remove with fewer than 5 rows
tracks <- filter(
  tracks,
  nrow >= 5
)

# make steps
tracks <- mutate(
  tracks,
  data = lapply(data, steps, keep_cols = "both")
)

Prepare alternate steps

# define alternate steps
n_alt_steps <- 9

# prepare alternate steps and extract covariates
tracks <- mutate(
  tracks,
  data = imap(data, function(.x, .y) {
    # messages
    message(
      glue("operating row {.y}")
    )
    # get 19 random potential steps
    amt::random_steps(.x, n = n_alt_steps)
  })
)

Data still contains long term sparrow tracking and swallow data.

save(tracks, file = "data/processed/data_for_ssf.Rds")

# # unlist and view
# data = unnest(
#   tracks,
#   cols = "data"
# )
#
# # write
# fwrite(
#   data, file = "data/results/data_alternative_patches.csv"
# )
#
# # make geometry
# data = mutate(
#   data,
#   seg = pmap(list(x1_, x2_, y1_, y2_), function(x1, x2, y1, y2) {
#     sf::st_linestring(
#       x = matrix(c(x1, y1, x2, y2), nrow = 2, byrow = 2)
#     )
#   })
# )
#
# # make sf
# data_steps = st_sf(
#   data,
#   sf_column_name = "seg",
#   crs = 2039
# )
#
# # write
# st_write(
#   data_steps,
#   dsn = "data/spatial/data_patch_alt_steps.gpkg"
# )

Read patches and get cov

# load("data/processed/data_for_ssf.Rds")
ndvi <- "data/rasters/raster_hula_ndvi_2039.tif"
vis <- "data/rasters/raster_hula_visibility.tif"

ndvi <- terra::rast(ndvi)
vis <- terra::rast(vis)

# landcover
lc <- st_read("data/spatial/hula_lc_vector")
lc <- mutate(
  lc,
  Name = case_when(
    Name == "C" ~ "W",
    TRUE ~ Name
  )
) %>%
  rename(
    lc = Name
  )

# crop  ndvi to lc
ndvi <- terra::crop(
  ndvi,
  vect(lc)
)

# rasterise
lc <- terra::rasterize(
  vect(lc),
  ndvi,
  field = "lc"
)

# rename points data
points <- rename(
  points,
  points_data = data
)

# link points data to tracks
tracks <- left_join(
  tracks,
  points
)

rand_pts_x <- rnorm(15, 0, 20)
rand_pts_y <- rnorm(15, 0, 20)

step_covs <- pmap(
  tracks,
  function(id, date, data, nrow, points_data) {
    # print message
    # message(sprintf("working on row %i", .y))

    #### env values in true patches ####
    # link real patch movements with real patch locations
    # identify the true patches in point data
    pdf <- points_data %>%
      select(
        patch, x, y
      )
    # get ndvi and vis
    # get ndvi values
    mean_vis <- terra::extract(
      vis,
      as.matrix(pdf[, c("x", "y")]),
      fun = mean
    )
    # get visibility values
    mean_ndvi <- terra::extract(
      ndvi,
      as.matrix(pdf[, c("x", "y")]),
      fun = mean
    )
    lc_samp <- terra::extract(
      lc,
      as.matrix(pdf[, c("x", "y")])
    )

    # add to real patch points and get summary
    pdf$ndvi <- mean_ndvi$ndvi
    pdf$vis <- mean_vis$raster_hula_visibility
    pdf$lc <- lc_samp$lc

    #### env values of alternative steps ####
    sim_pts <- select(
      data,
      x2_, y2_, case_, patch_end, patch_start
    ) %>%
      filter(
        !case_
      )

    # sample around ends of alt steps
    sim_pts <- mutate(
      sim_pts,
      spts = pmap(
        list(x2_, y2_),
        function(x2_, y2_) {
          tibble(
            x2_sim = x2_ + rand_pts_x,
            y2_sim = y2_ + rand_pts_y
          )
        }
      )
    )
    # unnest samples
    sim_pts <- unnest(
      sim_pts,
      cols = "spts"
    )
    # get ndvi values
    mean_vis <- terra::extract(
      vis,
      as.matrix(sim_pts[, c("x2_sim", "y2_sim")]),
      fun = mean
    )
    # get visibility values
    mean_ndvi <- terra::extract(
      ndvi,
      as.matrix(sim_pts[, c("x2_sim", "y2_sim")]),
      fun = mean
    )
    lc_samp <- terra::extract(
      lc,
      as.matrix(sim_pts[, c("x2_sim", "y2_sim")])
    )

    # assign values to alternative steps
    sim_pts$ndvi <- mean_ndvi$ndvi
    sim_pts$vis <- mean_vis$raster_hula_visibility
    sim_pts$lc <- lc_samp$lc

    #### step start covariates ####
    # these are always from a real patch
    df <- left_join(
      data,

      # drop coords when joining for start locs
      select(
        ungroup(pdf),
        patch, ndvi, vis, lc
      ),
      by = c("patch_start" = "patch")
    )
    df <- distinct(df)
    df <- rename(
      df,
      vis_start = vis,
      ndvi_start = ndvi,
      lc_start = lc
    )

    #### step end covariates for alternative steps ####
    df <-
      left_join(
        df,
        sim_pts,
        by = c(
          "x2_", "y2_", "case_",
          "patch_end", "patch_start"
        )
      ) %>%
      rename(
        ndvi_end = ndvi,
        vis_end = vis,
        lc_end = lc
      )

    #### step end covariates for real steps ####
    pdf <- mutate(
      pdf,
      case_ = TRUE
    )
    df <- left_join(
      df,
      pdf,
      by = c(
        "patch_end" = "patch",
        "case_"
      )
    ) %>%
      mutate(
        ndvi_end = if_else(
          is.na(ndvi_end), ndvi, ndvi_end
        ),
        vis_end = if_else(
          is.na(vis_end), vis, vis_end
        ),
        lc_end = if_else(
          is.na(lc_end), lc, lc_end
        )
      )

    df <- mutate(df, log_sl = log(sl_ + 1e-5)) %>%
      filter(!is.na(log_sl), !is.infinite(log_sl))

    # summarise covariates across real and alternative patches
    # return summary for comparison
    setDT(df)
    df_s <- df[, lapply(.SD, mean),
      .SDcols = c("ndvi_end", "ndvi_start", "vis_end", "vis_start"),
      by = c(
        "x1_", "y1_", "x2_", "y2_",
        "step_id_", "patch_start", "patch_end",
        "case_"
      )
    ]
    df_s_lc <- df[, .N,
      by = c(
        "x1_", "y1_", "x2_", "y2_",
        "step_id_", "patch_start", "patch_end",
        "case_", "lc_end"
      )
    ]
    df_s_lc[, prop := N / sum(N),
      by = c(
        "x1_", "y1_", "x2_", "y2_",
        "step_id_", "patch_start", "patch_end",
        "case_"
      )
    ]
    df_s_lc <- dcast(
      df_s_lc,
      x1_ + y1_ + x2_ + y2_ + step_id_ + patch_start + patch_end + case_ ~ lc_end,
      value.var = "prop", fill = 0
    )

    # handle missing names
    lc_names <- c("B", "O", "R", "T")
    df_s_lc[, c(setdiff(lc_names, names(df_s_lc))) := 0]

    # merge with data and return
    df_s <- merge.data.table(
      df_s,
      df_s_lc,
      by = intersect(names(df_s), names(df_s_lc))
    )

    as_tibble(df_s)
  }
)

tracks$step_covs <- step_covs

# remove rasters
rm(ndvi, vis)

# unnest data
tracks <- unnest(
  select(
    tracks, -data, -points_data,
  ),
  step_covs
)

# remove cols
tracks <- select(
  tracks,
  id,
  date,
  matches(c("x", "y")),
  matches(c("ndvi", "vis")),
  case_, step_id_,
  B, O, R, `T`, W
)

# save
write.csv(
  tracks,
  file = "data/results/data_ssf_step_covs.csv"
)