## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center"
)
set.seed(0)
library(dplyr)
library(ggplot2)
library(terra)
library(patchwork)
library(BayesianTools)
library(skytrackr)


## ----label = profiles, echo = FALSE, message = FALSE, warning = FALSE, fig.cap="Modelled and observed illuminance values in log(lux) for three days in early July 2021. Observed values are shown as light grey dots, modelled values for the exact location are shown as a full red line. Longitudinal offset data is shown as the blue full line, while latitudinal offset modelled data for the same period is shown as a full green line.", fig.width=7, out.width="80%"----
library(skylight)

cc876 <- skytrackr::cc876 |>
  dplyr::filter(
    date == "2021-08-02",
    measurement == "lux"
  ) |>
  tidyr::pivot_wider(
    names_from = "measurement",
    values_from = "value"
  )

# sequence of hours
dates <- seq(ISOdate(2021,08,01), ISOdate(2021,08,03), by = "hour")

# longitude offset
location_1 <- data.frame(
  date = dates,
  longitude = 3.74,
  latitude = 51.07
) |>
  skylight()

location_2 <- data.frame(
  date = dates,
  longitude = 20,
  latitude = 51.07
) |>
  skylight()

location_3 <- data.frame(
  date = dates,
  longitude = 3.74,
  latitude = 0
) |>
  skylight()

p <- ggplot() +
  geom_point(
    data = cc876,
    aes(
      date_time,
      log(lux)
    ),
    colour = "grey"
  ) +
  geom_line(
    data = location_1,
    aes(
     date,
     log(total_illuminance)
    ),
    colour = "red"
  ) +
    geom_line(
    data = location_2,
    aes(
     date,
     log(total_illuminance)
    ),
    colour = "blue"
  ) +
  geom_line(
    data = location_3,
    aes(
     date,
     log(total_illuminance)
    ),
    colour = "green"
  ) +
  labs(
    title = "Total illuminance early July 2021",
    y = "log(lux)",
    x = "Date"
  ) +
  theme_bw()

plot(p)

## ----label = ssf, echo = FALSE, fig.cap="A step selection function showing the probability of moving a given distance from the previous location.", fig.width=7, out.width="80%"----

# define a step selection distribution
ssf <- function(x, shape = 0.9, scale = 100, tolerance = 1500){
  # normalize over expected range with km increments
  norm <- sum(stats::dgamma(1:tolerance, shape = shape, scale = scale))
  prob <- stats::dgamma(x, shape = shape, scale = scale) / norm
  return(prob)
}

plot(1:1500, ssf(1:1500), xlab = "distance (km)", ylab = "probability", type ='l')


## ----label = light-profile, fig.cap="Light profile of a data logger", fig.width=7, out.width="80%"----
# read in the lux file
df <- stk_read_lux(
      system.file("extdata/cc876.lux", package="skytrackr")
    )

# plot the data as a static plot
df |> stk_profile()

## ----label = false-twilight, fig.cap = "Data with a false twilight", fig.width=7, out.width="80%"----
# introduce a false twilight
# by setting some twilight values very low
df <- df |> 
  mutate(
    value = ifelse(
      date_time > "2021-08-15 05:00:00" & date_time < "2021-08-15 12:00:00",
      0.1,
      value)
  )

# plot false twilight
df |> stk_profile()

## ----label = screen-profile, fig.cap="Screened light values, with offending days removed", fig.width=7, out.width="80%"----
# screen for twilight and other disturbing factors
df <- df |> stk_screen_twl(filter = TRUE)

# plot the updated dataset
df |> stk_profile()

## ----eval=FALSE---------------------------------------------------------------
# # visualize interactively
# df |> stk_profile(plotly = TRUE)
# 
# # subset to a particular date range
# df <- df |>
#    filter(
#    date >= "2019-08-31" & date <= "2020-04-15"
#    )

## ----label = land-mask, fig.cap = "Land mask to constrain location estimates"----
# generate a mask
mask <- stk_mask(
  bbox  =  c(-20, -40, 60, 60), #xmin, ymin, xmax, ymax
  buffer = 150, # in km
  resolution = 0.5
)

# plot the mask
plot(mask)

## ----eval = FALSE-------------------------------------------------------------
# # set a random seed
# # for reproducibility reasons
# set.seed(1)
# 
# # estimate locations
# locations <- df |>
#     skytrackr(
#       start_location = c(51.08, 3.73),
#       tolerance = 1500,
#       scale = log(c(0.00001,50)), # defaults
#       range = c(0.09, 140), # defaults
#       control = list(
#         sampler = 'DEzs',
#         settings = list(
#           burnin = 250, # default
#           iterations = 3000, # default
#           message = FALSE
#         )
#       ),
#       step_selection = ssf,
#       mask = mask,
#       plot = TRUE
#     )

## ----label = "preview", echo = FALSE, fig.cap="Dynamic preview during location estimation"----
knitr::include_graphics("skytrackr_preview.png")

## ----eval = FALSE-------------------------------------------------------------
# # map the output
# locations |> stk_map(bbox = c(-20, -40, 60, 60))

## ----label = "final-plot", echo = FALSE, fig.cap="Overview plot after location estimation"----
knitr::include_graphics("skytrackr_final_plot.png")

## ----eval = FALSE-------------------------------------------------------------
# if(!require(remotes)){install.packages("remotes")}
# remotes::install_github("bluegreen-labs/skytrackr")