Chapitre 16 Exercices corrigés

16.1 Exercice 1 - manipuler des objets sf

Créer un objet contenant les points de p qui intersectent le polygone a, à partir du code suivant :

library(sf)
library(tidyverse)
# polygone (a - orange)
a_poly <- st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a <- st_sfc(a_poly)
# lignes (l - bleues)
l1 <- st_multilinestring(list(rbind(c(0.5, -1), c(-0.5, 1))))
l2 <- st_multilinestring(list(rbind(c(.9, -.9), c(.5, 0))))
l <- st_sfc(l1, l2)
# multi-points (p - noirs)
p_matrix <- matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi <- st_multipoint(x = p_matrix)
p <- st_cast(st_sfc(p_multi), "POINT")

library(sf)
library(tidyverse)
# polygone (a - orange)
a_poly <- st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a <- st_sfc(a_poly)
# lignes (l - bleues)
l1 <- st_multilinestring(list(rbind(c(0.5, -1), c(-0.5, 1))))
l2 <- st_multilinestring(list(rbind(c(.9, -.9), c(.5, 0))))
l <- st_sfc(l1, l2)
# multi-points (p - noirs)
p_matrix <- matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi <- st_multipoint(x = p_matrix)
p <- st_cast(st_sfc(p_multi), "POINT")

Résultat attendu :

# solution st_filter
res <- st_sf(p) %>% 
  st_filter(a)
# solution tidyverse
res <- st_sf(p) %>% 
  filter(st_intersects(., a, sparse = FALSE))
# solution crochets
res <- st_sf(p)[a, , op = st_intersects]

res

# A tibble: 2 × 1
        p
  <POINT>
1 (0.5 0)
2   (1 1)

16.2 Exercice 2 : exploitation des données DVF en API

Le but de cet exercice va être d’exploiter les données DVF sur les transactions immobilières dans l’ancien et la carte des quartiers de Nantes pour obtenir des indicateurs des transactions par quartier.
On va utiliser pour DVF l’API mise en place par Christian Quest : http://api.cquest.org/dvf.

## Activation des packages
library(httr)
library(jsonlite)
library(sf)
library(tidyverse)
get_dvf <- GET("http://api.cquest.org/dvf?code_commune=44109")
dvf_content <- content(get_dvf, "text", encoding = "UTF-8")
dvf_json <- fromJSON(dvf_content)$resultats %>%
  # On ne garde que les données avec une géolocalisation valide, un prix et une surface renseignés.
  filter(!is.na(lon), !is.na(lat), !is.na(valeur_fonciere), !is.na(surface_relle_bati))
dvf <- st_as_sf(dvf_json, coords = c("lon", "lat"), crs = 4326)

## Activation des packages
library(tidyverse)
library(httr)
library(jsonlite)
library(sf)
get_dvf <- GET("http://api.cquest.org/dvf?code_commune=44109")
dvf_content <- content(get_dvf, "text", encoding = "UTF-8")
dvf_json <- fromJSON(dvf_content)$resultats %>%
  # On ne garde que les données avec une géolocalisation valide, un prix et une surface renseignés.
  filter(!is.na(lon), !is.na(lat), !is.na(valeur_fonciere), !is.na(surface_relle_bati))
dvf <- st_as_sf(dvf_json, coords = c("lon", "lat"), crs = 4326)

• Contour des quartiers de Nantes, ils proviennent de Nantes Métropole Open Data :https://data.nantesmetropole.fr

quartier_nantes <- st_read("https://data.nantesmetropole.fr/explore/dataset/244400404_quartiers-communes-nantes-metropole/download/?format=geojson&disjunctive.libcom=true&refine.libcom=Nantes&timezone=Europe/Berlin&lang=fr")
quartier_nantes <- st_set_crs(quartier_nantes, 4326)

quartier_nantes <- st_read("https://data.nantesmetropole.fr/explore/dataset/244400404_quartiers-communes-nantes-metropole/download/?format=geojson&disjunctive.libcom=true&refine.libcom=Nantes&timezone=Europe/Berlin&lang=fr")
quartier_nantes <- st_set_crs(quartier_nantes, 4326)

On veut produire les infos suivantes par quartier et année :

• Volume de ventes (nb)
  
• Pourcentage de maisons dans les ventes
  
• Prix moyen au m2 par type de bien

# Jointure spatiale pour récupérer les ventes par quartiers ---- 
dvf_avec_quartier <- st_join(dvf, quartier_nantes %>% select(nom)) %>%
  rename(quartier = nom)

# Calcul indicateurs----
library(lubridate)
dvf_filtre <- dvf_avec_quartier %>%
  st_drop_geometry() %>%
  filter(
    nature_mutation == "Vente",
    type_local %in% c("Appartement", "Maison")
  ) %>%
  mutate(
    date_mutation = ymd(date_mutation),
    annee_mutation = year(date_mutation),
    nb_ventes = 1
  ) 

# Calculs volumes, surfaces, prix totaux par quartier, par type de bien et par année 
stat1 <- dvf_filtre %>%
  group_by(quartier, type_local, annee_mutation) %>%
  summarise(across(c(nb_ventes, valeur_fonciere, surface_relle_bati), sum, na.rm = TRUE), .groups = "drop")

# Calculs volumes, surfaces, prix totaux par quartier et par année, ensemble maisons + appartements
stat2 <- dvf_filtre %>%
  group_by(quartier, annee_mutation) %>%
  summarise(across(c(nb_ventes, valeur_fonciere, surface_relle_bati), sum, na.rm = TRUE), .groups = "drop") %>% 
  mutate(type_local = "Ensemble") 

stat <- bind_rows(stat1, stat2)

Résultat attendu :

# Calcul volume des ventes
indicateurs1 <- stat %>%
  filter(type_local == "Ensemble") %>%
  select(quartier, annee_mutation, nb_ventes)

# Calcul pourcentage de maison dans les ventes
indicateurs2 <- stat %>%
  select(quartier, annee_mutation, type_local, nb_ventes) %>%
  pivot_wider(names_from = type_local, values_from = nb_ventes, values_fill = 0) %>%
  mutate(pourcentage_maison = 100 * Maison / Ensemble) %>%
  select(quartier, annee_mutation, pourcentage_maison)

# Calcul des prix au m2
indicateurs3 <- stat %>%
  select(quartier, annee_mutation, type_local, valeur_fonciere, surface_relle_bati) %>%
  mutate(prix_m2 = valeur_fonciere / surface_relle_bati) %>%
  select(quartier, annee_mutation, type_local, prix_m2) %>%
  pivot_wider(names_from = type_local, values_from = prix_m2) %>%
  rename_with(.cols = c(Appartement, Maison, Ensemble), .fn = ~paste0("prix_m2_", tolower(.x)))

# Assemblage des tables d'indicateurs
indicateurs <- reduce(list(indicateurs1, indicateurs2, indicateurs3), left_join)

# Réintroduction géométries
indicateurs <- quartier_nantes %>%
  select(quartier = nom) %>%
  left_join(indicateurs)
  
indicateurs %>%
glimpse()

Rows: 66
Columns: 8
$ quartier            <chr> "Ile de Nantes", "Ile de Nantes", "Ile de Nantes",…
$ annee_mutation      <dbl> 2014, 2015, 2016, 2017, 2018, 2019, 2014, 2015, 20…
$ nb_ventes           <dbl> 255, 346, 405, 464, 475, 221, 212, 203, 196, 224, …
$ pourcentage_maison  <dbl> 1.960784, 2.312139, 2.962963, 1.293103, 1.263158, …
$ prix_m2_appartement <dbl> 4227.436, 2720.641, 3583.049, 3496.280, 3574.651, …
$ prix_m2_maison      <dbl> 3208.333, 2961.019, 19695.106, 3519.273, 4144.621,…
$ prix_m2_ensemble    <dbl> 4187.948, 2728.811, 4208.121, 3496.700, 3587.453, …
$ geometry            <POLYGON [°]> POLYGON ((-1.518491 47.2105..., POLYGON ((…

16.3 Exercice 3 ggplot chap 10 : Visualisation des données DVF en API

Avec les résultats de l’exercice 2, produire les cartes du nombre de ventes et du prix au m2 des maisons en 2019 par quartier de Nantes.

Résultats attendus :

library(tidyverse)
library(gouvdown)
ggplot() +
  geom_sf(
    data = indicateurs %>%
      filter(annee_mutation == 2019),
    aes(fill = nb_ventes)
  ) + 
  theme_gouv_map() +
  scale_fill_gouv_continuous()

ggplot() +
  geom_sf(
    data = indicateurs %>%
      filter(annee_mutation == 2019),
    aes(fill = prix_m2_maison)
  ) + 
  theme_gouv_map() +
  scale_fill_gouv_continuous()

16.4 Exercice 4 : Assemblage de cartes sur dvf

A partir des données dvf 2014 et 2017 de la région Pays de la Loire contenues dans le package {variousdata} et les fonds de carte de {COGiter}, produire :

• une carte régionale à l’EPCI comprenant :
  • un dégrade de couleur sur l’évolution des prix au m2 des maisons entre 2014 et 2017,
    
  • un rond sur le volume des prix au m2 des maisons,
    
• un zoom sur les communes des principaux EPCI, c’est à dire une carte à la commune par EPCI de type Métropole (ME) ou Communauté urbaine (CU).

Puis, assembler ces différentes cartes sur un même graphique.

## Activation des packages
library(tidyverse)
library(sf)
library(lubridate)
library(variousdata) 
library(cowplot)
library(stringr)
library(COGiter)
library(gouvdown)

## Préparation des données

### Fonds de carte
epci_geo_r52 <- COGiter::epci_geo %>%
  left_join(COGiter::epci, by = "EPCI") %>% 
  filter(grepl("52", REGIONS_DE_L_EPCI))
epci_ppaux_r52 <- filter(epci_geo_r52, NATURE_EPCI %in% c('ME', "CU")) %>% 
  pull(EPCI)
# communes des EPCI principaux 
com_epci_ppaux_r52 <- COGiter::communes_geo %>% 
  left_join(COGiter::communes, by = "DEPCOM") %>% 
  filter(EPCI %in% epci_ppaux_r52)

### dvf

data("dvf_r52")
dvf_r52 <- dvf_r52 %>% 
  select(-c(NOM_DEPCOM:NOM_REG)) %>% 
  passer_au_cog_a_jour(code_commune = DEPCOM, garder_info_supra = TRUE, aggrege = FALSE)

# On ne conserve que les données valides de ventes de maisons, et on les tronque, en filtrant les données à 98% pour lisser les moyennes

dvf_r52_maisons <- dvf_r52 %>%
  filter(nature_mutation == "Vente", type_local == "Maison") %>%
  filter(!is.na(valeur_fonciere), !is.na(surface_reelle_bati)) %>%
  mutate(prix_m2 = valeur_fonciere / surface_reelle_bati) %>%
  arrange(prix_m2) %>%
  filter(between(row_number(), n() * .01, n() * .99)) %>%
  select(-prix_m2)

Il faut comme toujours procéder par étape.

Etape 1 : Calcul de l’évolution des prix et du nombre de ventes

• A l’EPCI

prix_m2_maisons_epci <- dvf_r52_maisons %>%
  select(EPCI, NOM_EPCI, date_mutation, valeur_fonciere, surface_reelle_bati) %>%
  mutate(
    n = 1,
    annee = year(date_mutation)
  ) %>%
  select(-date_mutation) %>%
  group_by(EPCI, NOM_EPCI, annee) %>%
  summarise(across(everything(), sum), .groups = "drop") %>%
  group_by(EPCI, NOM_EPCI) %>%
  mutate(
    prix_m2 = valeur_fonciere / surface_reelle_bati,
    evo_prix_m2 = 100 * prix_m2 / lag(prix_m2) - 100
  ) %>%
  filter(annee == 2017) %>%
  ungroup()

• A la commune

prix_m2_maisons_com <- dvf_r52_maisons %>%
  select(EPCI, DEPCOM, date_mutation, valeur_fonciere, surface_reelle_bati) %>%
  mutate(
    n = 1,
    annee = year(date_mutation)
  ) %>%
  select(-date_mutation) %>%
  group_by(EPCI, DEPCOM, annee) %>%
  summarise(across(everything(), sum), .groups = "drop") %>%
  group_by(EPCI, DEPCOM) %>%
  mutate(
    prix_m2 = valeur_fonciere / surface_reelle_bati,
    evo_prix_m2 = 100 * prix_m2 / lag(prix_m2) - 100
  ) %>%
  filter(annee == 2017) %>%
  ungroup()

• Intégration des données calculées aux fonds de carte

prix_m2_maisons_epci_sf <- epci_geo_r52 %>% 
  left_join(prix_m2_maisons_epci, by = c("EPCI", "NOM_EPCI")) %>%
  mutate(n = coalesce(n, 0))
prix_m2_maisons_com_sf <- com_epci_ppaux_r52 %>% 
  left_join(prix_m2_maisons_com, by = c("DEPCOM", "EPCI")) %>%
  mutate(n = coalesce(n, 0)) %>% 
  filter()

Etape 2 : Datavisualisation

• Carte à l’EPCI de la région

p <- ggplot(prix_m2_maisons_epci_sf) +
  geom_sf(aes(fill = evo_prix_m2)) +
  scale_fill_gouv_continuous(palette = "pal_gouv_div1") +
  stat_sf_coordinates(aes(size = n), alpha = .5) +
  theme_gouv_map(plot_title_size = 16, subtitle_size = 12, plot_margin = margin(0, 0, 0, 0), 
                 plot_title_margin = 1, caption_margin = 1, subtitle_margin = 0) +
  guides(size = "none") +
  labs(
    fill = "En %", title = "Evolution du prix des maisons neuves en euros par m2",
    subtitle = "Entre 2014 et 2017",
    caption = "source : DVF"
  )
p

• Zooms à la commune

# Création des cartes zoom EPCI avec une fonction

creer_zoom <- function(code_epci = "244400404") {
  
  nom_epci <- filter(epci_geo_r52, EPCI == code_epci) %>% 
    pull(NOM_EPCI) %>% 
    str_wrap(20)
  
  prix_m2_maisons_com_sf %>% 
    filter(EPCI == code_epci) %>% 
    ggplot() +
    geom_sf(aes(fill = evo_prix_m2)) +
    scale_fill_gouv_continuous(palette = "pal_gouv_div1") +
    stat_sf_coordinates(aes(size = n), alpha = .5) +
    theme_gouv_map(plot_title_size = 11, plot_margin = margin(0, 0, 0, 0), plot_title_margin = 1, caption_margin = 1, subtitle_margin = 0) +
    guides(size = "none", fill = "none") +
    labs(title = nom_epci)
  
  }

# Réalisation des zooms

zooms <- map(.x = epci_ppaux_r52, .f = creer_zoom)
zooms[[1]]

• Assemblage

plot_grid(p, plot_grid(zooms[[1]], zooms[[2]], zooms[[4]], nrow = 2, ncol = 3), 
          nrow = 2)

16.5 Exercice 5 : cartes pour le web

Adapter la carte régionale à l’EPCI de l’exercice 4, pour le web :

• au survol d’un EPCI, afficher son nom, le prix au m2 et son évolution 2014-2017,
  
• au survol du rond d’un EPCI, afficher son nom, le nb de ventes 2017.

library(ggiraph)
library(mapfactory)
p_web <- prix_m2_maisons_epci_sf %>% 
  mutate(
    sign_evol = if_else(evo_prix_m2 >=0, "+", ""),
    tooltip_fill = paste0(NOM_EPCI, " :\n", format_fr(x = prix_m2, dec = 0)," €/m2 en 2017\n", sign_evol,
                          format_fr(x = evo_prix_m2, pourcent = TRUE), " depuis 2014"),
    tooltip_size = paste0(NOM_EPCI, " :\n", n," ventes en 2017")
    ) %>% 
  ggplot() +
  geom_sf_interactive(aes(fill = evo_prix_m2, tooltip = tooltip_fill)) +
  scale_fill_gouv_continuous(palette = "pal_gouv_div1") +
  geom_point_interactive(stat = "sf_coordinates", aes(geometry = geometry, size = n, tooltip = tooltip_size), alpha = .5) +
  theme_gouv_map(plot_title_size = 16, subtitle_size = 12, plot_margin = margin(3, 3, 3, 3), 
                 plot_title_margin = 3, caption_margin = 2, subtitle_margin = 2) +
  guides(size = "none") +
  labs(
    fill = "En %", title = "Evolution du prix au m2 des maisons neuves",
    subtitle = "Entre 2014 et 2017",
    caption = "source : DVF"
  )
ggiraph(ggobj = p_web)

16.6 Exercice 6 : exercices bonus (ODD)

Faire une carte à de ronds proportionnels à partir des données sur les ODD

• En taille des ronds : le taux de mortalité des mères à naissance
• En couleur : l’évolution de la mortalité des mères à la naissance entre 2000 et 2015
• Rajouter deux zooms sur les continents Africain et Sud Américain

Source :

• les données : sdg_indicators du package variousdata
• le fond de carte : la table World du package tmap

16.6.1 Chargement des données et datapréparation

library(tidyverse)
library(tmap)
library(variousdata)
library(sf)
library(cowplot)
data("World")
data("sdg_indicators")
World_centroid <- st_centroid(World, of_largest_polygon = T)
sdg_indicators_sf_centroid <- World_centroid %>%
  left_join(sdg_indicators)

sdf_indicators_evo_2000_2015_sf <- sdg_indicators_sf_centroid %>%
  filter(timeperiod %in% c("2000", "2015"), !is.na(sh_sta_mmr)) %>%
  group_by(iso_a3) %>%
  mutate(evo_sh_sta_mmr = 100 * (sh_sta_mmr / lag(sh_sta_mmr) - 1)) %>%
  select(timeperiod, iso_a3, geoareaname, sh_sta_mmr, evo_sh_sta_mmr) %>%
  ungroup() %>%
  filter(timeperiod == "2015")

sdf_indicators_evo_2000_2015_sf %>%
  glimpse()

Rows: 167
Columns: 6
$ timeperiod     <fct> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2…
$ iso_a3         <fct> AFG, AGO, ALB, ARE, ARG, ARM, AUS, AUT, AZE, BDI, BEL, …
$ geoareaname    <fct> "Afghanistan", "Angola", "Albania", "United Arab Emirat…
$ sh_sta_mmr     <dbl> 396, 477, 29, 6, 52, 25, 6, 4, 25, 712, 7, 405, 371, 17…
$ evo_sh_sta_mmr <dbl> -64.00000, -48.37662, -32.55814, -25.00000, -13.33333, …
$ geometry       <POINT [°]> POINT (66.00365 33.84035), POINT (17.49892 -12.27…

16.6.2 Datavisualisation

16.6.3 Carte monde

Résultat intermédiaire :

map1 <- ggplot(data = sdf_indicators_evo_2000_2015_sf) +
  geom_sf(data = World, fill = "white") +
  geom_sf(aes(color = evo_sh_sta_mmr, size = sh_sta_mmr)) +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "light blue"),
    legend.position = "right"
  ) +
  guides(size = F) +
  scale_color_gradient2() +
  labs(
    title = "Evolution de la mortalité de la mère à la naissance",
    subtitle = "Entre 2000 et 2015",
    color = "Evolution\nen %"
  )
map1

16.6.4 Zoom sur les deux continents

On va utiliser la bbox pour définir les bornes de notre carte zoomée.

Résultat intermédiaire :

bbox_south_america <- World %>%
  filter(continent == "South America") %>%
  st_bbox()

map2 <- ggplot(data = sdf_indicators_evo_2000_2015_sf) +
  geom_sf(data = World, fill = "white") +
  geom_sf(aes(color = evo_sh_sta_mmr, size = sh_sta_mmr)) +
  coord_sf(
    xlim = c(bbox_south_america[1], bbox_south_america[3]),
    ylim = c(bbox_south_america[2], bbox_south_america[4])
  ) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "light blue"),
    legend.position = "none"
  ) +
  scale_color_gradient2()
map2

Résultat intermédiaire :

bbox_africa <- World %>%
  filter(continent == "Africa") %>%
  st_bbox()
map3 <- ggplot(data = sdf_indicators_evo_2000_2015_sf) +
  geom_sf(data = World, fill = "white") +
  geom_sf(aes(color = evo_sh_sta_mmr, size = sh_sta_mmr)) +
  coord_sf(
    xlim = c(bbox_africa[1], bbox_africa[3]),
    ylim = c(bbox_africa[2], bbox_africa[4])
  ) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "light blue"),
    legend.position = "none"
  ) +
  scale_color_gradient2()
map3

16.6.5 Assemblage

Résultat attendu :

ggdraw() +
  draw_plot(map1, x = 0, y = .45, width = 1, height = .55) +
  draw_plot(map2, x = 0.2, y = 0, width = .25, height = .45) +
  draw_plot(map3, x = 0.45, y = 0, width = .25, height = .45)