arbres-quebec/arbres-gbif.Rmd

---
title: "Arbres-gbif.Rmd"
author: "François Pelletier"
date: "26 janvier 2016"
output: html_document
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/',
                      echo=TRUE, warning=FALSE, message=FALSE)
library("checkpoint")
checkpoint("2016-03-07")
# Chargement des libraries
devtools::install_github("hadley/multidplyr")
source("packages.R")
```
# 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}
# Chargement du jeu de données
dA <- read.csv("ARBRE.csv",sep = "|", dec=",", stringsAsFactors = FALSE) %>%
# Correction des dates qui sont dans un format numérique YYYYMMDDHHMMSS
  mutate(DATE_PLANTE = parse_date_time(as.character(DATE_PLANTE/1000000), orders="%Y%m%d"))

# Créer un cluster de processeurs pour l'exécution en parallèle
cluster <- create_cluster(7)
set_default_cluster(cluster)
```

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 <- dA %>% parSapply(cl=cluster,FUN=class) %>% parSapply(cl=cluster, FUN="[",1)
# Affichage
transpose_row <- function(x,col_names)
  {
    df <- data.frame(x %>% names,
                     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}
dAtable <- dA[,dAclass == "character"] %>% parSapply(cl=cluster,freq)

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

noms_tables_frequences <- names(tables_frequences) %>% strsplit(split=".",fixed = TRUE) %>% 
  parSapply(cl=cluster,"[",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, dev="Cairo_svg"}
range <- quantile(dA$DIAMETRE,probs=c(0.02,0.98), na.rm=TRUE)
ggplot(data = (dA %>% filter(range[1] <= DIAMETRE & DIAMETRE < range[2] & POS_MESURE != "")),
       mapping = aes(x = DIAMETRE, fill=POS_MESURE)) + geom_histogram() + facet_wrap("POS_MESURE")
```

Date de plantation

```{r, dev="Cairo_svg"}
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 [http://www.gbif.org](Système Mondial d'Informations sur la Biodiversité GFIB])

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()

#data_json_gfib <- t(sapply(nomsUrlUniques$"nom_url",get_url_gfib))
#save(data_json_gfib,file="data_json_gfib.RData")
```

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 %>% 
                         '$'(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) %>% 
  (function(x) do.call(smartbind,x)) %>% 
  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 <- dA %>% partition() %>% merge(json_content4,by=c("NOM_LAT")) %>% collect()
```

## 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()
  
#json_media <- t(sapply(disct_speciesKey$speciesKey,get_url_media))
#save(json_media,file="json_media.RData")
load("json_media.RData")

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

json_media1$speciesKey <- disct_speciesKey

json_media2 <- json_media1[lapply(json_media1$json_content1,class) == 'data.frame',] 

json_media3 <- json_media2 %>% 
  '$'(json_content1) %>% 
  sapply(as.data.frame) %>% 
  reshape2::melt() %>% 
  select(value, type, format, identifier, references, title, description, source, creator, publisher, license) %>% 
  distinct(value)

json_media4 <- cbind(json_media3,json_media2$speciesKey)
```
## Joindre les données médias

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

## Ajout du quartier et de l'arrondissement

```{r}
qrtqc <- readOGR("QUARTIERS/", layer="QUARTIER") %>% spTransform(CRS("+proj=longlat +datum=WGS84"))
arrqc <- readOGR("ARROND/", layer="ARROND") %>% spTransform(CRS("+proj=longlat +datum=WGS84"))

names(qrtqc@data) <- paste0(names(qrtqc@data),"_QRT")
names(arrqc@data) <- paste0(names(arrqc@data),"_ARR")

coordinates(dA3) = ~ LONGITUDE + LATITUDE
proj4string(dA3) = CRS("+proj=longlat +datum=WGS84")

dA4.1 <- dA3 %>% over(qrtqc) %>% cbind(dA3@data)
dA4 <- dA3 %>% over(arrqc) %>% cbind(dA4.1) %>% cbind(dA3@coords)
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()) %>% 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)

#fond de la carte avec OpenStreetMap
#quebec_map <- get_map(location = c(range_long[1],range_lat[0],range_long[0],range_lat[1]), zoom=10, source = "osm")
#save(quebec_map, file="quebec_map.RData")
load("quebec_map.RData")
#objet ggmap

QuebecMap <- ggmap(quebec_map, base_layer = ggplot(aes(x = LONGITUDE, y = LATITUDE), data = select_for_map))

#Ajout des zones de densité d'arbres
map1 <- QuebecMap + scale_fill_gradient(low = "blue", high = "red")
```

### Carte Simple
```{r, dev="Cairo_png"}
map1 + 
  stat_density2d(aes(x = LONGITUDE, y = LATITUDE, fill = ..level..), geom = "polygon", data = select_for_map)
```

### Carte Composée
```{r, dev="Cairo_png"}
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 = ..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}
qrtqc@data$id = rownames(qrtqc@data)
qrtqc.df <- as.data.frame(qrtqc)
qrtqc.fort = fortify(qrtqc, region="id")
qrtqc.line = join(qrtqc.fort, qrtqc.df, by="id")

arrqc@data$id = rownames(arrqc@data)
arrqc.df <- as.data.frame(arrqc)
arrqc.fort = fortify(arrqc, region="id")
arrqc.line = join(arrqc.fort, arrqc.df, by="id")

ggmap_quartiers <-   ggplot(qrtqc.line) + 
  aes(long,lat,group=group,fill=NOM_QRT) + 
  geom_polygon() +
  geom_path(color="white") +
  coord_equal()

ggmap_arrond <- ggplot(arrqc.line) + 
  aes(long,lat,group=group,fill=NOM_ARR) + 
  geom_polygon() +
  geom_path(color="white") +
  coord_equal()

plot_data_quartiers <- 
  ddply(dA4, .(NOM_QRT,order), summarise, freq=length(NOM_QRT)) %>% filter(!is.na(order) && freq>=500)

plot_data_arrond <- 
  ddply(dA4, .(NOM_ARR,order), summarise, freq=length(NOM_ARR)) %>% 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=order, labels=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=order, labels=NULL) + 
  xlab("Ordre") +
  ylab("Fréquence") + 
  ggtitle("Ordre par arrondissement")
```

```{r, dev="Cairo_svg"}
ggmap_quartiers
```

```{r, dev="Cairo_svg"}
gg_freq_ordre_quartier
```

```{r, dev="Cairo_svg"}
ggmap_arrond
```

```{r, dev="Cairo_svg"}
gg_freq_ordre_arrond
```