---
title: "Arbres-gbif"
author: "François Pelletier"
format: pdf
---

```{r}
#| include: false
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/',
                      echo=TRUE, warning=FALSE, message=FALSE)
library(tidyverse)
library(sf)
sf_use_s2(FALSE)
library(httr)
library(jsonlite)
library(lubridate)
library(pander)
library(scales)
```

# Une forêt de données:
## Enrichir ses données à l'aide des catalogues de données ouvertes et des interfaces de programmation publiques

Ce projet utilise des données provenant du [catalogue de données de la Ville de Québec](http://donnees.ville.quebec.qc.ca/catalogue.aspx).
Ces données sont disponibles sous la [version 4.0 de la licence Creative Commons Attribution](https://creativecommons.org/licenses/by/4.0/deed.fr).

L'objectif de ce projet est de tirer un maximum d'information du jeu de données [Arbres répertoriés](http://donnees.ville.quebec.qc.ca/donne_details.aspx?jdid=26).

## Chargement du jeu de données

```{r}
#| echo: false
dA <- read.csv("ARBRE.csv", sep = "|", dec = ",", stringsAsFactors = FALSE) %>%
  mutate(DATE_PLANTE = parse_date_time(as.character(DATE_PLANTE / 1000000),
                                       orders = "%Y%m%d"))
```

Le jeu de données contient `r nrow(dA)` observations de `r ncol(dA)` variables.

## Classe de chaque variable

Le tableau suivant décrit la classe attribuée par R à chaque variable.

```{r}
#| results: asis
#| echo: false
dAclass <- sapply(dA, class) %>% sapply("[", 1)

transpose_row <- function(x, col_names) {
  df <- data.frame(
    names(x),
    x %>% unname %>% unlist
  )
  names(df) <- col_names
  df
}

pandoc.table(transpose_row(dAclass, c("Variable", "Classe")),
             style = "rmarkdown")
```

# Exploration de base du jeu de données

## Description des variables qualitatives

```{r}
#| results: asis
#| echo: false
freq <- function(x) sort(table(x), decreasing = TRUE)
dAtable <- dA[, dAclass == "character"] %>% sapply(freq)

tables_frequences <- dAtable %>%
  sapply(function(x) head(x[order(-x)], 10)) %>%
  lapply(transpose_row, c("Valeur", "Fréquence"))

noms_tables_frequences <- names(tables_frequences) %>%
  strsplit(split = ".", fixed = TRUE) %>%
  sapply("[", 1)

dump <- mapply(pandoc.table, t = tables_frequences,
               caption = noms_tables_frequences, style = "rmarkdown")
```

## Description des variables quantitatives

```{r}
#| results: asis
#| echo: false
dA %>%
  summarise(
    moyDIA = mean(DIAMETRE, na.rm = TRUE),
    minDIA = min(DIAMETRE, na.rm = TRUE),
    maxDIA = max(DIAMETRE, na.rm = TRUE),
    minLON = min(LONGITUDE, na.rm = TRUE),
    minLAT = min(LATITUDE, na.rm = TRUE),
    maxLON = max(LONGITUDE, na.rm = TRUE),
    maxLAT = max(LATITUDE, na.rm = TRUE)
  ) %>%
  transpose_row(c("Valeur", "Fréquence")) %>%
  pandoc.table(style = "rmarkdown")

dA %>%
  summarise(
    minDATE = min(DATE_PLANTE, na.rm = TRUE),
    maxDATE = max(DATE_PLANTE, na.rm = TRUE)
  ) %>%
  pandoc.table(style = "rmarkdown")
```

On remarque entre autres qu'il y a des erreurs dans les données de diamètres.
Certaines valeurs semblent inscrites en millimètres.

## Distributions

Diamètre des arbres

```{r}
range_vals <- quantile(dA$DIAMETRE, probs = c(0.02, 0.98), na.rm = TRUE)
ggplot(
  data = dA %>% filter(range_vals[1] <= DIAMETRE & DIAMETRE < range_vals[2] & POS_MESURE != ""),
  mapping = aes(x = DIAMETRE, fill = POS_MESURE)
) + geom_histogram() + facet_wrap("POS_MESURE")
```

Date de plantation

```{r}
ggplot(
  data = dA %>% filter(TYPE_ARBRE != "NON DISPONIBLE"),
  mapping = aes(x = DATE_PLANTE, fill = TYPE_ARBRE)
) + geom_histogram() + facet_wrap("TYPE_ARBRE")
```

# Enrichissement des variétés d'arbres

J'utilise les données provenant du [Système Mondial d'Informations sur la Biodiversité GFIB](http://www.gbif.org).

J'extrais d'abord les noms latins des espèces présentes dans la table avec `build_name` pour construire les URL de requêtes.
Puis, je vais les requêtes en lot avec un `mapply` sur la fonction `get_url_gfib`.

```{r}
get_url_gfib <- function(x)
  httr::GET(url = paste0("http://api.gbif.org/v1/species/match/?name=", x))

build_name <- function(x)
  gsub(pattern = " ",
      replacement = "+",
      x %>%
        strsplit("'") %>%
        unlist %>%
        '['(1) %>%
        trimws())

## Noms uniques (incluant la variété locale)
nomsUniques <-
  dA %>%
  select(NOM_LAT) %>%
  distinct() %>%
  mutate(nom_url = sapply(NOM_LAT, build_name) %>% tolower())

## Noms uniques pour construire les URL (excluant la variété locale qui ne se trouve pas dans GFIB)
nomsUrlUniques <-
  nomsUniques %>%
  select(nom_url) %>%
  distinct()
```

Je transforme ensuite les données recueillies dans le format JSON en une table que je pourrai joindre aux données source.

```{r}
load("data_json_gfib.RData")
json_content <- sapply(data_json_gfib %>%
                         as.data.frame %>%
                         pull(content), rawToChar)
json_content2 <- data.frame(
  nom_url = names(json_content),
  json_content, row.names = NULL, stringsAsFactors = FALSE
)
json_content3 <- json_content2$json_content %>%
  lapply(fromJSON, flatten = TRUE) %>%
  lapply(as.data.frame) %>%
  bind_rows() %>%
  cbind(nom_url = json_content2$nom_url)
json_content4 <- merge(json_content3, nomsUniques, by = c("nom_url"))
```

Je peux maintenant joindre ces nouvelles informations aux données source

```{r}
dA2 <- merge(dA, json_content4, by = c("NOM_LAT"))
```

## Médias par espèces

```{r}
get_url_media <- function(x)
  httr::GET(url = paste0("http://api.gbif.org/v1/species/", x, "/media"))

disct_speciesKey <- dA2 %>% select(speciesKey) %>% filter(!is.na(speciesKey)) %>% distinct()

load("json_media.RData")

json_media1 <- data.frame(
  json_content = sapply(json_media %>%
                          as.data.frame %>%
                          pull(content), rawToChar),
  stringsAsFactors = FALSE
) %>%
  mutate(json_content1 = json_content %>%
           lapply(fromJSON, flatten = TRUE) %>%
           sapply('[', "results"))

json_media1$speciesKey <- disct_speciesKey$speciesKey

json_media2 <- json_media1[sapply(json_media1$json_content1, function(x) is.data.frame(x) && nrow(x) > 0), ]

# Combine all media data.frames with their speciesKey
media_dfs <- list()
for (i in seq_len(nrow(json_media2))) {
  df <- json_media2$json_content1[[i]]
  if (is.data.frame(df) && nrow(df) > 0) {
    df$speciesKey <- json_media2$speciesKey[i]
    media_dfs[[length(media_dfs) + 1]] <- df
  }
}
json_media3 <- bind_rows(media_dfs)

# Keep only columns that actually exist
keep_cols <- intersect(
  c("value", "type", "format", "identifier", "references", "title",
    "description", "source", "creator", "publisher", "license", "speciesKey"),
  names(json_media3)
)
json_media3 <- json_media3 %>% select(all_of(keep_cols))
```

## Joindre les données médias

```{r}
dA3 <- merge(dA2, json_media3, all.x = TRUE)
```

## Ajout du quartier et de l'arrondissement

```{r}
qrtqc <- sf::read_sf("QUARTIERS/QUARTIER.shp") %>% sf::st_transform(4326)
arrqc <- sf::read_sf("ARROND/ARROND.shp") %>% sf::st_transform(4326)

# Rename columns to avoid conflicts
names(qrtqc)[names(qrtqc) != "geometry"] <- paste0(
  names(qrtqc)[names(qrtqc) != "geometry"], "_QRT"
)
names(arrqc)[names(arrqc) != "geometry"] <- paste0(
  names(arrqc)[names(arrqc) != "geometry"], "_ARR"
)

# Convert points to sf and perform spatial join
dA3_sf <- sf::st_as_sf(dA3, coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
dA4.1 <- sf::st_join(dA3_sf, qrtqc)
dA4 <- sf::st_join(dA4.1, arrqc)

# Add coordinates back and remove geometry
dA4 <- dA4 %>%
  mutate(
    LONGITUDE = sf::st_coordinates(.)[, 1],
    LATITUDE = sf::st_coordinates(.)[, 2]
  ) %>%
  sf::st_drop_geometry()

save(dA4, file = "dA4.RData")
readr::write_csv(dA4, "arbres-augmented.csv")
```

## Arbre le plus courant par quartier (ayant une photo disponible)

```{r}
#| results: asis
count_arbre_arr <- dA4 %>%
  filter(identifier != "" & !is.na(NOM_QRT)) %>%
  select(NOM_QRT, scientificName, identifier) %>%
  group_by(NOM_QRT, scientificName, identifier) %>%
  summarise(freq = n(), .groups = "drop") %>%
  group_by(NOM_QRT) %>%
  top_n(n = 1)

pandoc.table(count_arbre_arr %>% select(NOM_QRT, scientificName, freq))
```

```{r}
#| results: asis
count_arbre_arr %>%
  mutate(image = paste0("[", scientificName, "](", identifier, ")")) %>%
  select(NOM_QRT, image) %>%
  t %>%
  pandoc.table()
```

## Localisation sur une carte

```{r}
select_for_map <- dA4 %>%
  select(LONGITUDE, LATITUDE, order)
range_long <- range(select_for_map$LONGITUDE)
range_lat <- range(select_for_map$LATITUDE)

# Without ggmap (not available), use a clean ggplot2 background
QuebecMap <- ggplot(data = select_for_map, aes(x = LONGITUDE, y = LATITUDE))
map1 <- QuebecMap + scale_fill_gradient(low = "blue", high = "red")
```

### Carte Simple

```{r}
map1 +
  stat_density2d(aes(x = LONGITUDE, y = LATITUDE, fill = after_stat(level)),
                 geom = "polygon", data = select_for_map)
```

### Carte Composée

```{r}
table_order <- with(select_for_map, table(order))
select_for_map2 <- table_order %>%
  as.data.frame() %>%
  merge(select_for_map, all.y = TRUE) %>%
  filter(Freq > 500)
# Garder seulement les niveaux actifs
select_for_map2$order2 <- factor(select_for_map2$order)

map1 +
  stat_density2d(aes(x = LONGITUDE, y = LATITUDE, fill = after_stat(level)),
                 geom = "polygon", data = select_for_map2) +
  facet_wrap(facets = "order2")
```

## Arbres recensés par quartier

Source: [plotting polygon shapefiles](https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles)

```{r}
ggmap_quartiers <- ggplot(qrtqc) +
  aes(fill = NOM_QRT) +
  geom_sf(color = "white")

ggmap_arrond <- ggplot(arrqc) +
  aes(fill = NOM_ARR) +
  geom_sf(color = "white")

plot_data_quartiers <-
  dA4 %>%
  group_by(NOM_QRT, order) %>%
  summarise(freq = n(), .groups = "drop") %>%
  filter(!is.na(order) & freq >= 500)

plot_data_arrond <-
  dA4 %>%
  group_by(NOM_ARR, order) %>%
  summarise(freq = n(), .groups = "drop") %>%
  filter(!is.na(order) & freq >= 500)

gg_freq_ordre_quartier <- ggplot(data = plot_data_quartiers,
                                 aes(x = order, y = freq, fill = order)) +
  geom_bar(position = "stack", stat = "identity") +
  facet_wrap(facets = "NOM_QRT", ncol = 4) +
  scale_x_discrete(breaks = NULL) +
  xlab("Ordre") +
  ylab("Fréquence") +
  ggtitle("Ordre par quartier")

gg_freq_ordre_arrond <- ggplot(data = plot_data_arrond,
                               aes(x = order, y = freq, fill = order)) +
  geom_bar(position = "stack", stat = "identity") +
  facet_wrap(facets = "NOM_ARR") +
  scale_x_discrete(breaks = NULL) +
  xlab("Ordre") +
  ylab("Fréquence") +
  ggtitle("Ordre par arrondissement")
```

```{r}
ggmap_quartiers
```

```{r}
gg_freq_ordre_quartier
```

```{r}
ggmap_arrond
```

```{r}
gg_freq_ordre_arrond
```