---
title: "Untitled"
format: html
---

## Libraries

```{r}
library(tidyr)
library(dplyr)
library(terra)
library(mregions2)
library(biooracler)
library(stringr)
library(tibble)
library(catboost)
library(caret)
library(blockCV)
library(sf)
```

## Download

### Sea

Загружаем данные по Баренцеву морю

```{r}
barentsz_mrgid = 4247
geo = gaz_geometry(barentsz_mrgid, format = "sfc") |> vect()
```


Фиксируем охват

```{r}
bbox = ext(geo)
```

### Bio-Oracle

Фиксируем доступные в Bio-Oracle слои

```{r}
layers = list_layers()
```


Фиксируем слои, на которые нет прогнозных данных. Их мы не будем использовать в обучении и предсказании.

```{r}
# Нет прогнозных данных :/
removed_layers_ids = c(
    "par_mean_baseline_2000_2020_depthsurf",
    "kdpar_mean_baseline_2000_2020_depthsurf",
    "chl_baseline_2000_2018_depthmax",
    "chl_baseline_2000_2018_depthmean",
    "chl_baseline_2000_2018_depthmin"
)
```


Фиксируем слой с характеристиками поверхности дна. Они будут одинаковыми для всех временных срезов и сценариев.

```{r}
constant_layers_ids = c("terrain_characteristics")

constant_layers = layers |>
    filter(dataset_id %in% constant_layers_ids)
```


Фиксируем динамические слои. Они являются основой проводимого анализа. Разбиваем код на отдельные поля: переменная, сценарий, период, глубина.

```{r}
dynamic_layers = layers |>
    filter(! dataset_id %in% c(constant_layers_ids, removed_layers_ids)) |> 
    separate_wider_delim(dataset_id, delim = "_", names = c("var", "scenario", "year_star", "year_end", "depth"), cols_remove = FALSE)
```


Создаём заготовку для формирования срезов. Пока включаем только пространственную составляющую.

```{r}
lon_range = c(bbox$xmin, bbox$xmax)
lat_range = c(bbox$ymin, bbox$ymax)
constraints = list(
    longitude = lon_range,
    latitude = lat_range
)
```


Загружаем данные о поверхности дна.

```{r}
terrain_raster = download_layers(
    constant_layers$dataset_id[1],
    constraints = constraints,
    directory = "data/bio-oracle/terrain_characteristics"
)
# Переименовываем слои в полные названия, чтобы потом поля понятно назывались
names(terrain_raster) = longnames(terrain_raster)
```


Формируем функцию загрузки среза данных. Аргументами являются сценарий климатического развития и декада. Параметры охвата беруться из констант (см. `lon_range`, `lat_range`).

```{r}
download_slice = function(scenario_value, decade_start) {
    scenario_layers = dynamic_layers |> 
        filter(scenario == scenario_value)

    time_point = paste0(decade_start, "-01-01T00:00:00Z")

    slice_constraints = list(
        time = c(time_point, time_point),
        longitude = lon_range,
        latitude = lat_range
    )

    download_dir = file.path("./data/bio-oracle", scenario_value, decade_start)
    dir.create(download_dir, recursive = TRUE, showWarnings = FALSE)

    slice_rasters = sapply(
        scenario_layers$dataset_id,
        function(id) download_layers(
            id,
            constraints = slice_constraints,
            directory = download_dir
        ),
        simplify = TRUE
    )

    return(slice_rasters)
}
```


Выполняем загрузку среза на текущую декаду. Он будет использоваться в качестве контекстных данных для обучения.

```{r}
baseline_rasters = download_slice("baseline", 2010)
```

### IUCN
<!-- Publish dataset -->

```{r}
mammals_range = vect("data/iucn/MAMMALS_MARINE_ONLY.shp")
```


Загружаем данные об ареалах видов в пределах озвата Баренцевого моря.

#### Crop to the sea

```{r}
# mammals_range_barentsz = mammals_range |> 
#     crop(bbox)

# writeVector(mammals_range_barentsz, "data/iunc/mammals_barentsz.shp")

mammals_range_barentsz = vect("./data/iucn/mammals_barentsz.shp")
```

## Transform

Приводим загруженные данные к виду, пригодному для машинного обучения

### Bio-Oracle

Чтобы в df потом были нормальные названия колонок
```{r}
baseline_brick = rast(baseline_rasters)

baseline_brick_depths = baseline_brick |> 
    names() |> 
    str_extract("depth[:alpha:]+")
baseline_brick_longnames = baseline_brick |> longnames()
# baseline_brick_varnames = baseline_brick |> varnames() // краткая запись называний слоёв
# Человекочитаемые названия слоёв -> полей в датафрейме (см. следующий блок)
names(baseline_brick) = paste(baseline_brick_longnames, baseline_brick_depths)
```


Трансформируем растр в датафрейм. Сохраняем идентификаторы ячеек для надёжного установления соответствия между контекстными и целевыми данными.

```{r}
baseline_df = c(baseline_brick, terrain_raster) |> 
    as.data.frame(cells = TRUE)
```

### IUCN

Трансформируем данные об ареалах в растр с параметрами, аналогичными растрам параметров местообитаний.

```{r}
template_raster = baseline_rasters[[1]]
orca_range = mammals_range_barentsz |> 
    subset(mammals_range_barentsz$sci_name == "Orcinus orca")

orca_range_boundaries_buffer = orca_range |> 
    as.lines() |> 
    buffer(50000) |> 
    rasterize(template_raster, 1)

orca_range_raster = orca_range |> 
    rasterize(template_raster, "", background=0) |> 
    mask(orca_range_boundaries_buffer, inverse = TRUE)

plot(orca_range_raster)
```

```{r}
orca_df = orca_range_raster |> 
    as.data.frame(cells = TRUE) |> 
    mutate(target = factor(if_else(layer == 0, "0", "1"), levels = c("0", "1"))) |> 
    select(-layer)
```

```{r}
mammals_brick = mammals_range_barentsz$sci_name |> 
    sapply(function(sci_name) {
        mammals_range_barentsz |> 
            subset(mammals_range_barentsz$sci_name == sci_name) |> 
            rasterize(baseline_rasters[[1]], field="sci_name")
    }) |> 
    rast()
```

```{r}
plot(mammals_brick$`Orcinus orca`)
```


Трансформируем растр ареалов животных в датафрейм с идентификаторами ячеек. Так как параметры растров аналогичны, можно установить однозначное соответствие между ячейкой ареала и параметрами местообитания в данной ячейке.

```{r}
mammals_df = mammals_brick |> as.data.frame(cells = TRUE)
```

```{r}
m_r = ifel(is.na(mammals_brick), 0, 1)
mammals_df2 = m_r |> 
    as.data.frame(xy = TRUE, cells = TRUE)
```

```{r}
# mammals_sf = mammals_df2
#     st_as_sf(coords = c("x", "y"), crs = 4326)
# b = cv_spatial(
#     x = mammals_sf,
#     column = "Orcinus orca",
#     # r = baseline_brick,
#     # k = 3,
#     size = 500000,
#     # hexagon = TRUE,
#     # selection = "random"
# )
# train_cells = mammals_sf[b$folds_list[[1]][[1]], ]$cell
# test_cells = mammals_sf[b$folds_list[[1]][[2]], ]$cell
```

```{r}
train_cells = mammals_df2 |>
    select(cell, x, y) |> 
    filter(y > 75) |> 
    pull(cell)
test_cells = mammals_df2 |> 
    select(cell, x, y) |> 
    filter(y < 75) |> 
    pull(cell)
```

```{r}
cv_plot(
    cv = b, # a blockCV object
    x = mammals_sf, # sample points
)
```

<!-- Выбрать только виды, у которых присутсвие меньше 80% -->

## Learn

Learn what conditions does animals prefer

```{r}
learn = function(species_name, hyperparam) {
    # Подготовка данных
    species_df = mammals_df |> 
        select(cell, all_of(species_name)) |> 
        rename(target = species_name)

    baseline_species_df = baseline_df |> 
        left_join(species_df, by = "cell") |> 
        mutate(target = factor(if_else(is.na(target), "0", "1"), levels = c("0", "1")))# |> 
        #select(-cell)

    # Разделение выборки на обучающую и тестовую
    # set.seed(123)
    # train_index = createDataPartition(baseline_species_df$target, p = 0.8, list = FALSE)

    train_df = baseline_species_df |> 
        filter(cell %in% train_cells)
    train_features = train_df %>% select(-cell, -target) # параметры
    train_labels = train_df$target # наличие ареала
    train_pool = catboost.load_pool(data = train_features, label = train_labels)

    test_df = baseline_species_df |> 
        filter(cell %in% test_cells)
    test_features = test_df %>% select(-cell, -target)
    test_labels = test_df$target
    test_pool = catboost.load_pool(data = test_features, label = test_labels)


    # Обучение
    model = catboost.train(train_pool, test_pool = test_pool, params = hyperparam)

    return(model)
}
```

```{r}
baseline_orca_df = baseline_df |> 
    left_join(orca_df, by = "cell") |> 
    filter(!is.na(target))
```

```{r}
learn_orca = function(species_name, hyperparam) {
    baseline_species_df = baseline_df |> 
        left_join(orca_df, by = "cell") |> 
        filter(!is.na(target))

    train_df = baseline_species_df |> 
        filter(cell %in% train_cells)
    train_features = train_df %>% select(-cell, -target) # параметры
    train_labels = train_df$target # наличие ареала
    train_pool = catboost.load_pool(data = train_features, label = train_labels)

    test_df = baseline_species_df |> 
        filter(cell %in% test_cells)
    test_features = test_df %>% select(-cell, -target)
    test_labels = test_df$target
    test_pool = catboost.load_pool(data = test_features, label = test_labels)


    # Обучение
    model = catboost.train(train_pool, test_pool = test_pool, params = hyperparam)

    return(model)
}
```

```{r}
fit_params <- list(
  iterations = 100,
  learning_rate = 0.03,
  depth = 6,
  loss_function = 'Logloss',    # Standard for binary classification
  eval_metric = 'AUC',          # Good metric for imbalanced data
  random_seed = 42,
  verbose = 10                 # Print progress every 100 iterations
  # od_type = "Iter",           # Optional: Early stopping
  # od_wait = 50
)
```

```{r}
m_orca = learn("Orcinus orca", fit_params)
```

```{r}
m_orca_2 = learn_orca("Orcinus orca", fit_params)
```

```{r}
preds_prob <- catboost.predict(m_orca, pool_test, prediction_type = 'Probability')
preds_class <- ifelse(preds_prob > 0.5, 1, 0)
a = confusionMatrix(factor(preds_class), factor(labels_test))
```

```{r}
importance <- catboost.get_feature_importance(model, pool = pool_train)
```


## Predict

Predict changes in habitat ranges using different scenarios

```{r}
ssp585_2090_rasters = download_slice("ssp585", 2090)
```

```{r}
ssp585_2090_brick = rast(ssp585_2090_rasters)
ssp585_2090_brick_depths = ssp585_2090_brick |> 
    names() |> 
    str_extract("depth[:alpha:]+")
ssp585_2090_brick_longnames = ssp585_2090_brick |> longnames()
# baseline_brick_varnames = baseline_brick |> varnames() // коды longnames

names(ssp585_2090_brick) = paste(ssp585_2090_brick_longnames, ssp585_2090_brick_depths)
```

```{r}
ssp585_2090_df = c(ssp585_2090_brick, terrain_raster) |> 
    as.data.frame(cells = TRUE)
```

```{r}
ssp585_2090_features = ssp585_2090_df |> select(-cell)
```

```{r}
ssp585_2090_pool <- catboost.load_pool(data = ssp585_2090_features)
```

```{r}
preds_prob <- catboost.predict(model, ssp585_2090_pool, prediction_type = 'Probability')
preds_class <- ifelse(preds_prob > 0.5, 1, 0)
```

```{r}
ssp585_2090_prediction = ssp585_2090_df |> 
    mutate(prediction = preds_class) |> 
    select(cell, prediction)
```

```{r}
a = orca_df |> 
    left_join(ssp585_2090_prediction, by = "cell") |> 
    mutate(target = if_else(is.na(target), 0, 1)) |> 
    mutate(diff = prediction - target)
```

```{r}
hist(a$diff)
```

```{r}
# Select only your numeric predictors (environmental layers)
# Calculate correlation matrix

cor_matrix <- baseline_df |> 
    select(-cell, -landmass, -coastline) |> 
    cor()

# Find attributes that are highly corrected (ideal cutoff is debatable, 0.85 or 0.9 is common)
high_cor_features <- findCorrelation(cor_matrix, cutoff = 0.9)

# Print names of features suggested for removal
print(colnames(baseline_df)[high_cor_features])
```

```{r}
r = rast(baseline_brick)
r[a$cell] = a$diff
```

```{r}
values(r) = NA
```

```{r}
plot(r$`Long-term maximum AirTemperature depthsurf`)
```