sdm_1.qmd

---
title: "sdm_1"
format: html
editor_options: 
  chunk_output_type: console
---

## Goals

Create initial SDM using:

-   **Observations**\
    Species with most OBIS observations. Later try temporally subsetting to look for migratory patterns. And run for all.

-   **Environment**\
    Environmental data from `sdmpredictors`. Need to know whether demersal or pelagic. For now extract static climatic environmental predictors for model fitting and prediction. Later extract dynamic predictors synchronous in time for model fitting and predict with hind-/now-/fore-cast or climatic seasonal snapshot.

-   **Model**\
    Use Maxent.

## References

-   [Species distribution modeling --- R Spatial](https://rspatial.org/sdm/)
    -   old but more: [Species distribution modeling --- R Spatial (raster)](https://rspatial.org/raster/sdm/)

## Setup

```{r}
# libraries ----
options(java.parameters = "-Xmx8000m")
librarian::shelf(
  arrow, dismo, dplyr, DT, glue, here, leafem, 
  leaflet.extras, readr, rJava, robis, sdmpredictors, terra, sf, mapview,
  quiet = T)
options(readr.show_col_types = F)

# variables -----
dir_data <- "/Users/bbest/My Drive/projects/mbon-sdm/data"
obis_prq <- glue("{dir_data}/raw/obis.org/obis_20230726.parquet")

# functions -----
get_sp_occ_obis_prq <- function(
    aphia_id,
    obis_prq  = "/Users/bbest/My Drive/projects/mbon-sdm/data/raw/obis.org/obis_20230726.parquet",
    cols_keep = c(
      "id",
      "phylum",
      "class",
      "taxonRank",
      "scientificName",
      "AphiaID",
      "date_mid",
      "decimalLongitude",
      "decimalLatitude",
      "depth",
      "individualCount",
      "flags")) {
  # get species occurrences from OBIS parquet file
  
  # TODO: add caching per species request/args
  
  # return only observations with valid coordinates and year
  o <- open_dataset(obis_prq) |> 
    filter(
      !is.na(date_mid),
      !is.na(decimalLongitude),
      !is.na(decimalLatitude),
      AphiaID == !!aphia_id) |> 
    select(all_of(cols_keep)) |> 
    collect() |>
    mutate(
      date_mid =as.POSIXct(
        date_mid/1000, origin = "1970-01-01",tz = "GMT") |> 
        as.Date()) |> 
    st_as_sf(
      coords = c("decimalLongitude", "decimalLatitude"), 
      crs = 4326)
}
```

## Species candidates

```{r}
# species with most observations
sp_gull <- read_csv("data/obis_top-species.csv") |> 
  arrange(desc(n_obs)) |>
  slice(1)
  
# non-bird species with most observations
sp_herring <- read_csv("data/obis_top-species.csv") |> 
  filter(class != "Aves") |> 
  arrange(desc(n_obs)) |>
  slice(1)

# species with ~100 observations (and most in Class)
sp_jelly <- read_csv("data/obis_top-species.csv") |> 
  filter(n_obs > 240) |> 
  arrange(n_obs) |>
  slice(1)
```

* `sp_gull`: [_Larus fuscus_ lesser black-backed gull | WoRMS](https://www.marinespecies.org/aphia.php?p=taxdetails&id=137142)\
  n_obs: 2,162,536
  
* `sp_herring`: [_Clupea pallasii_ Pacific herring | WoRMS](https://www.marinespecies.org/aphia.php?p=taxdetails&id=151159)\
  n_obs: 1,961,525

* `sp_jelly`: [_Haliclystus auricula_ kaleidoscope jellyfish | WoRMS](https://www.marinespecies.org/aphia.php?p=taxdetails&id=135322)\
  n_obs: 242
  
* `sp_rwhale`: [_Eubalaena glacialis_ North Atlantic right whale | WoRMS](https://www.marinespecies.org/aphia.php?p=taxdetails&id=159023)\
  n_obs: 1,000
  * [_Eubalaena glacialis_ | OBIS](https://obis.org/taxon/159023)
    * [OBIS Mapper](https://mapper.obis.org/?taxonid=159023#)
  * [_Eubalaena glacialis_ | AquaMaps](https://aquamaps.org/preMap2.php?cache=1&SpecID=ITS-Mam-180537)

## obs: N. Atlantic right whale

N Atlantic right whale (_Eubalaena glacialis_)

```{r}
# right whale: surface, migratory, endangered
sp_rwhale <- tibble(
  AphiaID = 159023)

# prep single species observations ---
sp       <- sp_rwhale
aphia_id <- sp$AphiaID

# get observations
obs <- get_sp_occ_obis_prq(aphia_id)
obs

# plot observations
mapView(obs) %>%
  .@map |> 
  leaflet.extras::addFullscreenControl()
```

## env: `sdmpredictors`

```{r}
librarian::shelf(
  sdmpredictors, skimr, terra)

# see ../aquamaps-downscaled/index.qmd for creation
env_tif  <- here("../aquamaps-downscaled/data/bio-oracle.tif")

stopifnot(file.exists(env_tif))

env <- rast(env_tif)
names(env) <- c(
  "Temp",
  "Salinity",
  "PrimProd",
  "IceCon")
plot(env)

obs_env <- extract(env, obs, cells=T) |> 
  as_tibble() |> 
  mutate(
    presence = 1)

bk_cells <- setdiff(cells(env), obs_env$cell)
bk_env <- tibble(
  presence = 0,
  cell     = bk_cells) |> 
  bind_cols(
    values(env)[bk_cells,]) |> 
  distinct(pick(all_of(names(env))), .keep_all = T)
# rows: 6,181,644 -> 6,178,309 after distinct() and na.omit()

set.seed(42) 

d_env <- bind_rows(
  obs_env,
  bk_env |> 
    slice(sample(1:nrow(bk_env), nrow(obs_env)) ) ) |> 
  select(
    presence, cell,
    all_of(names(env))) |> 
  na.omit()

tail(d_env)

# table(d_env$presence)
#       0       1 
# 6178319    7355
#    0    1 
# 7359 7355 

d_env
```


## model: `predicts::maxent()`

### example

```{r}
librarian::shelf(
  fs, ggplot2, leaflet, predicts, 
  quiet = T)
#?predicts::MaxEnt
# MaxEnt()
# This is MaxEnt_model version 3.4.3

# get predictor variables
f <- system.file("ex/bio.tif", package="predicts")
preds <- rast(f)
plot(preds)

# file with presence points
occurence <- system.file("/ex/bradypus.csv", package="predicts")
occ <- read.csv(occurence)[,-1]

# witholding a 20% sample for testing 
fold     <- folds(occ, k=5)
occtest  <- occ[fold == 1, ]
occtrain <- occ[fold != 1, ]

# fit model
me <- MaxEnt(preds, occtrain)

# see the MaxEnt results in a browser:
# me
dir_copy(attr(me, "path"), here("me_bradypus"))
```

Results: [maxent.html](./data/me_bradypus/maxent.html)

```{r}
# plot showing importance of each variable
plot(me, main="me: Variable contribution")

# TODO: try categorical, not in preds
# # use "args"
# me2 <- MaxEnt(preds, occtrain, factors='biome', args=c("-J", "-P"))
# 
# # plot showing importance of each variable
# plot(me2, main="me2: Variable contribution")

# response curves

d <- tibble()
for (v in names(preds)){  # v = names(preds)[2]
  pr <- partialResponse(me, var=v)
  d <- bind_rows(
    d, 
    pr |> 
      rename(value = 1) |> 
      mutate(
        var = v))
  # plot(pr, type="l", las=1)
  # TODO: lattice ggplot
}
g <- d |> 
  ggplot(aes(x = value, y = p)) +
  geom_line() +
  facet_wrap(
    ~var, scales = "free") +
  theme_bw()
g
plotly::ggplotly(g)
# TODO: convert to function

# pr2 <- partialResponse2(me, var="bio1", var2 = "bio5")
# plot(pr2, type="l", las=1)

# predict to entire dataset
r <- predict(me, preds) 
plot(r)
points(occ)

# with some options:
r <- predict(me, preds, args=c("outputformat=raw"))
plot(r)
points(occ)


#testing ----

# background sample
bg <- backgroundSample(preds, 1000)

#simplest way to use 'evaluate'
e1 <- pa_evaluate(me, p=occtest, a=bg, x=preds)

# alternative 1
# extract values
pvtest <- data.frame(extract(preds, occtest))
avtest <- data.frame(extract(preds, bg))

e2 <- pa_evaluate(me, p=pvtest, a=avtest)

# alternative 2 
# predict to testing points 
testp <- predict(me, pvtest) 
head(testp)
testa <- predict(me, avtest) 

e3 <- pa_evaluate(p=testp, a=testa)
e3
threshold(e3)

plot(e3, 'ROC')
```

### rwhale

```{r}
# get predictor variables
preds <- env
plot(preds)

# witholding a 20% sample for testing 
occ <- obs |> 
  mutate(
    lon = st_coordinates(geometry)[,1],
    lat = st_coordinates(geometry)[,2]) |>
  st_drop_geometry() |> 
  select(lon, lat)
fold     <- folds(occ, k=5)
occtest  <- occ[fold == 1, ]
occtrain <- occ[fold != 1, ]

# fit model
system.time({
  me <- MaxEnt(preds, occtrain)
})
# Warning message:
# In .local(x, p, ...) :
#   3 (0.19%) of the presence points have NA predictor values


# see the MaxEnt results in a browser:
# see the MaxEnt results in a browser:
# me
dir_copy(attr(me, "path"), here("me_rwhale"))
```

Results: [maxent.html](./data/me_rwhale/maxent.html)

```{r}
# plot showing importance of each variable
plot(me, main="me: Variable contribution")

# TODO: try categorical, not in preds
# # use "args"
# me2 <- MaxEnt(preds, occtrain, factors='biome', args=c("-J", "-P"))
# 
# # plot showing importance of each variable
# plot(me2, main="me2: Variable contribution")

# response curves
system.time({
d <- tibble()
for (v in names(preds)){  # v = names(preds)[2]
  pr <- partialResponse(me, var=v)
  d <- bind_rows(
    d, 
    pr |> 
      rename(value = 1) |> 
      mutate(
        var = v))
  # plot(pr, type="l", las=1)
  # TODO: lattice ggplot
}
g <- d |> 
  ggplot(aes(x = value, y = p)) +
  geom_line() +
  facet_wrap(
    ~var, scales = "free") +
  theme_bw()
g
})
# plotly::ggplotly(g)
# TODO: convert to function

# pr2 <- partialResponse2(me, var="bio1", var2 = "bio5")
# plot(pr2, type="l", las=1)

# predict to entire dataset
system.time({
  r <- predict(me, preds)
})
plot(r)
points(occ)

# with some options:
# system.time({
#   r <- predict(me, preds, args=c("outputformat=raw"))
# })
# plot(r)
# points(occ)

plet(r, tiles = leaflet::providers$Esri.NatGeoWorldMap) |> 
  addCircleMarkers(
    data = occ, 
    lng = ~lon, 
    lat = ~lat,
    opacity     = 0.3,
    fillOpacity = 0.1,
    radius = 2,
    weight = 1)
```


## mdl: maxnet

```{r}
#| eval: false

librarian::shelf(
  maxnet, 
  quiet = T)

system.time({
  m <- maxnet(
    p    = d_env |> pull(presence), 
    data = d_env |> select(-presence, -cell) |> as.data.frame())
}) # present/absent:7,359/7,355: 50.149s

plot(m, type="cloglog")

predict(m, new)

env |> terra::add()

env
# 2160 * 4320 = 9331200, 5
d_env <- values(env, dataframe = T, na.rm=F)
# d_env <- d_env |> 
#   select(-pred)
dim(d_env) # 9,331,200       4

ncell(env) # 9,331,200

d_env_notna <- na.omit(d_env)
dim(d_env_notna) # 6,183,457

x <- d_env_notna

intersect(attr(x, "na.action"), attr(x, "row.names"))
class(attr(x, "row.names"))


p <- predict(m, d_env_notna, clamp = T, type = "logistic") #, clamp=T, type=c("link","exponential","cloglog","logistic"), ...)
length(p[,1]) # 6,183,457

p <- predict(m, raster::stack(env), clamp = F, type = "logistic") #, clamp=T,

env$pred <- NA
env$pred[attr(x, "row.names")] <- p[,1]

plot(env$pred)
setValues(env, p)

# why are these different lengths?
as.numeric(p)
ncell(env)

length(cells(env))
length(as.numeric(p))

plot(env$pred)
mod <- maxnet(p, data, maxnet.formula(p, data, classes="lq"))
plot(mod, "tmp6190_ann")
```


## vs Kernel Density Estimates (KDE)

* [Kernel density estimates for tidy and geospatial data in the eks package](https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html)


```{r}
librarian::shelf(
  eks, 
  quiet = T)
```