3_assign_grid.Rmd

---
title: "Assign grid to occurrences"
author:
- Damiano Oldoni
- Peter Desmet
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    number_sections: true
---

In this document we assign grid to occurrences.

# Setup

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

Load libraries:

```{r load_libraries}
library(tidyverse)      # To do datascience
library(sp)             # To work with geospatial data
library(here)           # To find files
library(glue)           # To write queries
library(RSQLite)        # To interact with SQlite databases
```

# Get geographic coordinates and coordinate uncertainty

Name and path of `.sqlite` file:

```{r name_path}
key <- "0031758-231002084531237"
sqlite_file <- paste(key, "occurrence.sqlite", sep = "_")
sqlite_path <- here::here("data", "interim", sqlite_file)
```

Table name:

```{r define_table_name}
table_name <- "occ"
```

Open connection to database:

```{r open_connection_to_db}
sqlite_occ <- dbConnect(SQLite(), dbname = sqlite_path)
```

Retrieve geographic coordinates, `decimalLatitude` and `decimalLongitude` and coordinate uncertainty, `coordinateUncertaintyInMeters`:

```{r get_geodata}
query <- glue_sql("SELECT {`cols`*} FROM {table}",
              cols = c("decimalLatitude", 
                       "decimalLongitude",
                       "coordinateUncertaintyInMeters"),
              table = table_name,
              .con = sqlite_occ
)
geodata_df <- 
  dbGetQuery(sqlite_occ, query) %>%
  dplyr::as_tibble() %>%
  dplyr::mutate(coordinateUncertaintyInMeters = as.double(.data$coordinateUncertaintyInMeters))
```

Number of occurrences:

```{r n_occs}
nrow_geodata_df <- nrow(geodata_df)
nrow_geodata_df
```

Preview:

```{r preview_geodata}
geodata_df %>% head(10)
```

Number of occurrences per each value of `coordinateUncertaintyInMeters`:

```{r n_occ_per_uncertainty}
geodata_df %>%
  group_by(coordinateUncertaintyInMeters) %>%
  count() %>%
  arrange(desc(n))
```

Number of occurrences without coordinate uncertainty or set to zero:

```{r n_occs_without_uncertainty}
geodata_df %>%
  filter(is.na(coordinateUncertaintyInMeters) | 
           coordinateUncertaintyInMeters == 0) %>%
  count()
```

We assign 1000 meters to these occurrences:

```{r assign_fix_uncertainty}
geodata_df <- 
  geodata_df %>%
  mutate(coordinateUncertaintyInMeters = if_else(
    is.na(coordinateUncertaintyInMeters) | coordinateUncertaintyInMeters == 0,
    1000.0,
    coordinateUncertaintyInMeters
  ))
```

We apply the same changes in SQLite table. To speed up the process we deactivate database synchronization and wrap the execution of the query in a transaction:

```{r assign_1000_as_uncertainty_on_db}
# speed up by deactivating db synchronization
query <- glue_sql(
  "PRAGMA synchronous = OFF",
  .con = sqlite_occ
)
dbExecute(sqlite_occ, query)

# Start transaction
dbBegin(sqlite_occ)
query <- glue_sql(
            "UPDATE {table} SET {`column`} = 1000 WHERE 
            {`column`} = '' OR {`column`} = 0",
            table = table_name,
            column = "coordinateUncertaintyInMeters",
            .con = sqlite_occ)
dbExecute(sqlite_occ, query)
# Commit transaction
dbCommit(sqlite_occ)
```

Get all possible values of coordinate uncertainty from SQLite database:

```{r get_unique_values_uncertainty}
query <- glue_sql("SELECT DISTINCT {`column`} FROM {table}",
              table = table_name,
              column = "coordinateUncertaintyInMeters",
              .con = sqlite_occ
)
uncertainty_df <- 
  dbGetQuery(sqlite_occ, query) %>%
  as_tibble()
uncertainty_df
```

Are there occurrences with 0 uncertainty left?

```{r zero_in_uncertainty_df}
0 %in% uncertainty_df
```

Are there occurrences without uncertainty?

```{r zero_in_uncertainty_df}
NA %in% uncertainty_df
```

# Assign grid to occurrences

We use the official grid of Belgium at 1x1km resolution as provided by EEA.

## Project geographic coordinates, assign occurrence within uncertainty circle and assign grid cell

In this section we do the following:
1. Project latitude and longitude by using the projection of the grid. We transform GBIF data which have coordinate reference system equal to EPSG code 4326 to Lambert projection with EPSG code 3035.
2. Assign the occurrences randomly within the circle with radius equal to `coordinateUncertaintyInMeters`
3. Define the grid cell the occurrence belongs to. 

Due to memory issues, we perform this three operations in chunks on a temporary tsv file containing coordinates and coordinate uncertainties.

First, we create such file:

```{r create_temp_file}
temp_file_coords <- here::here(
  "data", 
  "interim", 
  paste0(key, "_coordinates_and_uncertainty_epsg_4326.tsv")
)
col_names_geodata_df <- names(geodata_df)
write_tsv(geodata_df, temp_file_coords, na = "")
remove(geodata_df)
```

Set random number generator seed (this helps reproducibility). We use the unique identifier of the reserved [Zenodo dataset's DOI](https://doi.org/10.5281/zenodo.10058400) which the occurrence cube will be published to:

```{r set_seed}
set.seed(10058400)
```

We define the function to apply to each chunk:

```{r transform_to_3035_assign_pts_in_circle}
reproject_assign <- function(df, pos){
  
  # Step 1: reprojection
  nrow_df <- nrow(df)
  coordinates(df) <- ~decimalLongitude+decimalLatitude
  proj4string(df) <- CRS("+init=epsg:4326")
  df <- spTransform(df, CRS("+init=epsg:3035"))
  colnames(df@coords) <- c("x", "y")
  
  # Step 2: assign occurrence within uncertainty circle
  df@data <-
    df@data %>%
    mutate(random_angle = runif(nrow_df, 0, 2*pi))
  df@data <-
    df@data %>%
    mutate(random_r = sqrt(runif(
      nrow_df, 0, 1)) * coordinateUncertaintyInMeters)
  df@data <-
  df@data %>%
    mutate(x = df@coords[, "x"],
           y = df@coords[, "y"])
  df@data <-
    df@data %>%
    mutate(x = x + random_r * cos(random_angle),
           y = y + random_r * sin(random_angle)) %>%
    select(-c(random_angle, random_r))
  
  # Step 3: Find grid cell the occurrence belongs to
  df@data <-
    df@data %>%
    mutate(eea_cell_code = paste0(
    "1km", 
    "E", floor(x/1000), 
    "N", floor(y/1000))) %>%
    select(x, y, coordinateUncertaintyInMeters, eea_cell_code)
  return(df@data)
}
```

And finally we apply reprojection to all occurrences:

```{r reproject_assign_coords_by_chunk}
chunk_size <- 1000000
geodata_df <- read_tsv_chunked(
  file = temp_file_coords, 
  callback = DataFrameCallback$new(reproject_assign), 
  chunk_size = chunk_size,
  col_types = cols(.default = col_double()),
  na = ""
)
```

Preview:

```{r preview_geodata_df_data}
geodata_df %>% head(n = 10)
```

Before proceeding, we can delete the temporary file created at the beginning of this section:

```{r delete_temp_file}
file.remove(temp_file_coords)
```

## Add grid cell code to sqlite file

We can now add the column `eea_cell_code` to the table `occ` of sqlite file. We first create the new column `eaa_cell_code` in the table:

```{r add_eaa_cellcode_to_sqlite}
new_col <- "eea_cell_code"
query <- glue_sql("ALTER TABLE {table} ADD COLUMN {colname} {type}",
                  table = table_name,
                  colname = new_col,
                  type = "CHARACTER",
                  .con = sqlite_occ
)
dbExecute(sqlite_occ, query)
```

And then we populate it with the values in `geodata_df$eaa_cell_code`. This step can take long:

```{r add_values_cellcode}
# Start transaction
dbBegin(sqlite_occ)
dbExecute(
  sqlite_occ,
  glue_sql(
    "UPDATE {table} SET {`column`} = :eea_cell_code WHERE _ROWID_ = :id",
    table = table_name,
    column = new_col,
    .con = sqlite_occ),
  params = data.frame(
  eea_cell_code = geodata_df$eea_cell_code,
  id = rownames(geodata_df))
)
# Commit transaction
dbCommit(sqlite_occ)
```

Preview:

```{r preview}
query <- glue_sql("SELECT * FROM {table} WHERE _ROWID_ <= 10",
                  table = table_name,
                  .con = sqlite_occ)
dbGetQuery(sqlite_occ, query) %>%
  select(eea_cell_code , everything())
```

Close connection:

```{r close_connection}
dbDisconnect(sqlite_occ)
```