# /!\/!\/!\ à exécuter si nécessaire
# install.packages("tidyverse", type = "win.binary")

# 0 - Nettoyer son env de travail
rm(list=ls())

# Charger librairies
library(dplyr)
library(sf)
library(leaflet)
library(ggplot2)
library(ggspatial)
library(osmdata)
library(tidyr)
library(readxl)

# 0 -Fixer le repertoire de travail
setwd("~/git/MASTER_1_R/data")

# 1 - Ouvrir les donnees
communes <- st_read("communes_vulnerables_polygn_5490.shp", quiet=T)

# 2 - Sélectionner une commune
communes <- st_read("communes_vulnerables_polygn_5490.shp", quiet=T) %>% 
  slice(1)

1 - Visualisation interactive avec Leaflet

# 3 - Afficher les communes sur une carte
leaflet() %>% 
  addProviderTiles("Esri.WorldImagery") %>% 
  addPolygons(data = communes %>% st_transform(4326), color='red', fillColor = "red", fillOpacity = .2)

2 - Créer une grille

# 4 - Récuperer la bounding box de la zone d'étude
bbox <- st_bbox(communes)

grid <- communes %>%
  st_make_grid(cellsize = 200) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id_grid = row_number()) %>% 
  st_intersection(.,communes[,c("NOM","geometry")]) %>% 
  select(-NOM)

# 3 - Afficher les communes sur une carte
leaflet() %>% 
  addProviderTiles("Esri.WorldImagery") %>% 
  addPolygons(data = grid %>% st_transform(4326), color='red', fillColor = "red", fillOpacity = .2)

3 - Télécharger les surfaces bâties

# Depuis OSM (library osmdata)
db_buildings <- opq(bbox=st_bbox(communes %>% st_transform(4326))) %>%
  add_osm_feature(key = "building") %>%
  osmdata_sf()

# Ne garder que les polygones
buildings <- db_buildings$osm_polygons
rm(db_buildings)

# Réduire la dimension de la table
buildings <- buildings %>%
  select(osm_id, geometry)

# Afficher le résultats
map1 <- ggplot(communes) + 
  geom_sf() + 
  geom_sf(data = buildings)+
  annotation_scale(style='ticks',width_hint=0.1)+
  annotation_north_arrow(style = north_arrow_fancy_orienteering,
                         height = unit(.75,"cm"),
                         width = unit(.75,"cm"),
                         location = "tl", which_north = "true")+
  annotate('text', x = Inf, y = -Inf, label = '\u00a9 OpenStreetMap\n\u00a9 INSEE',
         hjust = 1, vjust = -1, color = 'black', size = 3) +
  theme_void()
map1

4 - Calculer des superficies

# Calucler les superficies
buildings$area <- sf::st_area(buildings)

# Calucler les centroides
buildings_centroids <- sf::st_centroid(buildings)

5 - Réaliser une jointure spatiale

# Reprojeter les centroides en 5490
buildings_centroids <- buildings_centroids %>% 
  st_transform(5490)

# Joindre l'ID des cellules de la grille à chaque centroide
buildings_centroids <- buildings_centroids %>% 
  st_join(., grid)

buildings_centroids %>% head()

Simple feature collection with 6 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 638324.4 ymin: 1763639 xmax: 639648.7 ymax: 1764195
Projected CRS: RGAF09 / UTM zone 20N
             osm_id            area id_grid                 geometry
41380356   41380356 480.87172 [m^2]      26 POINT (639379.9 1763671)
41380357   41380357 221.21241 [m^2]      26 POINT (639323.2 1763660)
41380358   41380358  19.22201 [m^2]      21 POINT (638324.4 1763648)
125187318 125187318 101.56157 [m^2]      62 POINT (638975.3 1764195)
125187319 125187319  74.56967 [m^2]       8 POINT (639442.2 1763639)
125187323 125187323  66.12055 [m^2]      28 POINT (639648.7 1763819)

5 - Calcul des statistiques

# Grouper puis sommer les superficies
stats_byID <- buildings_centroids %>% 
  st_drop_geometry() %>% 
  group_by(id_grid) %>% 
  summarize(built_up=sum(area))

stats_byID %>% head()

# A tibble: 6 × 2
  id_grid built_up
    <int>    [m^2]
1       6    4302.
2       7    2257.
3       8     445.
4      21    1432.
5      22     213.
6      23    1909.

6 - Afficher les résultats

# Joindre les statistiques au spatialDataFrame
grid <- grid %>% 
  merge(stats_byID)

# Reprojeter les données
grid_4326 <- grid %>%
  st_transform(4326)

# Discretiser la série statistique
percentile <- quantile(as.numeric(grid$built_up),c(.25,.50,.75))
bins <- c(0, percentile[1],percentile[2],percentile[3], Inf)
pal <- colorBin("YlOrRd", domain = grid$built_up, bins = bins)

# Afficher les résultats
map2 <- leaflet() %>% 
  addProviderTiles("Esri.WorldImagery") %>% 
  addPolygons(data= grid_4326,
    fillColor = ~pal(built_up),
    weight = 2,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7)
map2