Aide à l’utilisation de la BD CASSMIR

Base de données spatialisées

L’ensemble des données spatiales est stocké dans un fichier commun au format .gpkg. Il est déposé dans un dossier local nommé “CASSMIR_Outputs”. Le fichier geopackage est composé des trois jeux de données produits au niveau communal, du carroyage 1km et du carroyage 200m. L’import des données dans la console se réalise avec la fonction st_read du package sf, en sélectionnant le niveau géographique souhaité (couche).

library(sf)
library(ggplot2)
library(readxl)
library (openxlsx)
library(tidyverse)
library(SpatialPosition)
library(cartography)
library(tmap)
library(linemap)
library(tidyverse)
library(osmdata)
library(devtools)
# install_github('riatelab/mapinsetr') # voir : https://github.com/riatelab/mapinsetr
library(mapinsetr)
library(mapview)
library(reshape2)
library(ggpubr)

# Import
## Objets disponibles dans le fichier
st_layers("CASSMIR_Outputs/CASSMIR_SpatialDataBase.gpkg")
## Import des couches
cassmir_Grid200m <- st_read("CASSMIR_Outputs/CASSMIR_SpatialDataBase.gpkg", quiet = TRUE, layer = "Grid200m")
cassmir_Grid1km<- st_read("CASSMIR_Outputs/CASSMIR_SpatialDataBase.gpkg", quiet = TRUE, layer = "Grid1km")
cassmir_Communes <- st_read("CASSMIR_Outputs/CASSMIR_SpatialDataBase.gpkg", quiet = TRUE, layer = "Communes")

ls(cassmir_Communes)

Avant toute manipulation de la donnée, il est conseillé de se reporter à la documentation relative aux métadonnées de la base CASSMIR

Selon les objectifs de l’analyse, l’utilisateur de la base CASSMIR est amené à utiliser des filtres et à transformer le format du tableau de données. Chaque utilisation doit être réalisée avec une parfaite connaissance de la composition de la base CASSMIR.

knitr::kable(cassmir_Grid200m[1:5,1:5], format="markdown")
knitr::kable(cassmir_Grid1km[1:5,1:5], format="markdown")
knitr::kable(cassmir_Communes[1:5,1:5], format="markdown")

Base de données sur les groupes de population

L’ensemble des données sur les groupes de population est stocké dans un fichier commun au format .csv. Ce fichier est composé d’un jeu de données unique. L’import des données dans la console se réalise avec la fonction read.csv2.

# Import
cassmir_GroupesPop <-  read.csv2("CASSMIR_Outputs/CASSMIR_GroupesPopDataBase.csv", stringsAsFactors=FALSE)

ls(cassmir_GroupesPop)

Sans restriction, pour faciliter l’utilisation de la base CASSMIR sur les groupes de population, il est conseillé de procéder de la manière suivante :

  • Couper le tableau en trois tableau distincts en fonction des Types de groupe

  • Si l’utilisateur souhaite étudié uniquement les acquéreurs/vendeurs, un filtre doit être appliqué en conséquence.

  • Pour travailler sur une série de variables particulières, la sélection des variables doit s’effectuer avec vigilance en prenant en compte toutes les variables à la nomenclature majuscule commune.

Chaque utilisation doit être réalisée avec une parfaite connaissance de la composition de la base CASSMIR.

Elaboration d’un modèle cartographique

Objectif général

Pour récolter ces données géométriques, l’utilisateur peut suivre les liens individuels de téléchargement présentés ci-dessous ou télécharger directement l’ensemble des données stockées dans le dossier CASSMIR_OpenDataRaws.zip en suivant ce lien. Ce dossier peut ensuite être déposé à l’emplacement du projet CASSMIR.

Dans un objectif de création d’un template cartographique, aux données CASSMIR sont associées des couches d’informations géographiques utiles pour permettre leur contextualisation géographique (fleuves, routes, forêts etc).

On vise à produire les couches géographiques suivantes :

Pour le fonds de carte régional :

Des couches sur le maillage territorial régional : - Les IRIS de l’Île-de-France (couche nommée iris) - Les communes de l’Île-de-France (couche com) - Les départements de l’Île-de-France (couche dep) - La limite régionale de l’Île-de-France (couche reg) - Les carreaux de 200m et d’1km de l’Île-de-France (couches grid200 et grid1km) - Une sélection sur les carreaux de 200m et d’1km de l’Île-de-France avec des transactions immobilières observées dans la période étudiées, correspondant à “l’espace du marché” (couches MarketSpaceGrid200 et MarketSpaceGrid1km)

A cela s’ajoute des couches OSM sur des éléments d’habillage cartographique : - La Seine (couche nommée Seine) - La Marne (couche nommée Marne) - Les principaux axes routiers (couche nommée Autoroutes) - Les axes ferroviaires de longue portée (couche nommée TrainsGL) - Les axes ferroviaires à portée régionale longue portée (couche nommée TrainsRegion) - Le principal aéroport international, Roissy Charles-de-Gaulle (couche nommée Aeroports) - Les trois principales forêts (en foncton de la surface occupée) de la région (couche nommée Forets) - Les principales villes (en fonction de population) de la région (couche nommée Villes)

Pour le fonds de carte centré sur la Métropole du Grand Paris (MGP) : Des couches sur le maillage territorial : - Limites de la MGP (couche nommée MGP_box) - Les IRIS (couche nommée IrisMgp_box) - Les communes (couche ComMgp_box) - Les carreaux de 200m et d’1km (couches Grid200Mgp_box et Grid1kMgp_box) - Une sélection sur les carreaux de 200m et d’1km de l’Île-de-France avec des transactions immobilières observées dans la période étudiées, correspondant à “l’espace du marché” (couches MarketSpaceGrid200Mgp_box et MarketSpaceGrid1kmMgp_box)

A cela s’ajoute des couches sur des éléments d’habillage cartographique : - La Seine (couche nommée SeineMgp) - La Marne (couche nommée MarneMgp)

Origine des données conventionnelles sur les mailles et limites géographiques

Deux sources de données est utilisée pour récolter les géométries de référence : le zonage IRIS 2017 et les grilles de 200m de côté et de 1km de côté (version 2015), construits par l’INSEE et délivré par l’IGN et l’INSEE. En ce qui concerne les entités communales, la base de données GEofla produite par l’IGN est mobilisée.

INSEE : IRIS et grilles en entrée

C’est la dernière version des IRIS connue en février 2020 qui est utilisée ici. Il s’agit de la publication Contours…IRIS édition 2019 mis à disposition par l’IGN sur son portail dédié à la livraison des jeux de données en OpenData. La date de référence de ce découpage géographique est le 17 janvier 2020.

Les niveaux supérieurs de la hiérarchie administrative française sera reconstituée à partir de ce découpage en Iris. Il s’agira en particulier de reconstituer et d’extraire les entités communales, les limites des départements et de la région Île-de-France.

2 grilles de référence pour la France sont aussi utilisées :

Les données inclues dans ces trois grilles sont issues du millésime 2015 Filosofi. Le millésime Filosofi de l’année 2015 est élaboré à partir des revenus perçus en 2015 qui ont été déclarés en 2016 et de la taxe d’habitation au 1er janvier 2016. En revanche aucune donnée attributaire n’est associée au découpage IRIS.

IGN : Geofla pour les limites communales

C’est la version 2016 des communes qui est utilisée ici. Il s’agit de la publication Geofla®Communes édition 2016 France Métropolitaine mis à disposition par l’IGN sur son portail dédié à la livraison des jeux de données en OpenData. La date de référence de ce découpage géographique est le 1er juin 2016.

APUR et EPT Grand-Orly Seine Bièvre

Il s’agit ici des communes de la Métropole du Grand Paris au 1er janvier 2016, suite à la loi NOTRe. Les données et la documentation sont disponibles via ce lien. Elles ont été mises en ligne en décembre 2018 et sont livrées en Opendata. Ce fichier a été produit par l’APUR et l’EPT Grand-Orly Seine Bièvre

Au total,il s’agit de 5 ressources volumineuses : 114 Mo pour le fond géographique IRIS pour la France métropolitaine (.shp) ; 141 Mo pour la grille kilométrique (.shp) et 1,2 Go pour la grille de 200 m (.shp) ; 98 Mo pour la couche des communes de france métropolitaine (.shp) et la documentation associée ; 1.52 Mo pour la couche (.shp) des communes de la métropole du Grand Paris.

Les données sont importées sans modification.

Préparation des objets géométriques de référence du projet

Les objets sont importés du dossier local “CASSMIR_OpenDataRaws” dans leur format d’origine et sont exportés dans le dossier local “CASSMIR_Geom” en format Geopackage.

# - Import des données en local pour fonds de carte

### Couche IRIS 
 iris <- st_read("CASSMIR_OpenDataRaws/CONTOURS-IRIS_2-1__SHP__FRA_2020-01-01/CONTOURS-IRIS/1_DONNEES_LIVRAISON_2020-01-00139/CONTOURS-IRIS_2-1_SHP_LAMB93_FXX-2019", quiet = TRUE)  
# iris <- st_read("CASSMIR_OpenDataRaws/CONTOURS-IRIS_2-1__SHP__FRA_2019-01-01/CONTOURS-IRIS/1_DONNEES_LIVRAISON_2020-01-00139/CONTOURS-IRIS_2-1_SHP_LAMB93_FXX-2020", quiet = TRUE)

iris <- st_transform(iris, 2154)


### Couche Communes 
com <- st_read("CASSMIR_OpenDataRaws/GEOFLA_2-2_COMMUNE_SHP_LAMB93_FXX_2016-06-28/GEOFLA/1_DONNEES_LIVRAISON_2016-06-00236/GEOFLA_2-2_SHP_LAMB93_FR-ED161/COMMUNE/COMMUNE.shp", quiet = TRUE) 
  
com <- st_transform(com, 2154)

### Grilles de référence
grid200<- st_read("CASSMIR_OpenDataRaws/carreaux200/Filosofi2015_carreaux_200m_metropole.shp", quiet = TRUE) 

grid1k <- st_read("CASSMIR_OpenDataRaws/carreaux1000/Filosofi2015_carreaux_1000m_metropole.shp", quiet = TRUE)

## Métropole du Grand Paris 

### Communes
comMgp <- st_read("CASSMIR_OpenDataRaws/communes-mgp/Communes_MGP.shp", quiet = TRUE) %>% st_transform(., crs = "+init=epsg:2154")


# Agrégation Région
iris$DEP <- substr(iris$INSEE_COM,1,2) # création variable département à partir des codes communes

# On garde uniquement les entités franciliennes et tranformation crs2154
irisIDF <- iris %>%
   filter(DEP =="75" |DEP =="77" |DEP =="78" |DEP =="91" |DEP =="92" |DEP =="93" | DEP =="94" |DEP =="95") %>%
  st_transform(., crs = "+init=epsg:2154")
  
# Obtenir la limite régionale
RegIDF <- iris %>%
   filter(DEP =="75" |DEP =="77" |DEP =="78" |DEP =="91" |DEP =="92" |DEP =="93" | DEP =="94" |DEP =="95") %>%
  summarise(geom = st_union(geometry)) %>%
  st_transform(., crs = "+init=epsg:2154")

# Agrégation départements
Dep <- com %>%
  group_by (CODE_DEPT) %>%
  summarise() %>%
  st_transform(., crs = "+init=epsg:2154") # Obtenir les départements

## Création d'une Bbox
bb <- st_bbox(RegIDF)

bboxRegIDF <- st_as_sfc(st_bbox(bb + c(-40000,-7000,5000,7000), crs = 2154)) # On tire volontairement la bbox vers la gauche du plan en prévision d'un emplacement pour le carton
bboxRegIDF <-  st_transform(bboxRegIDF, crs = "+init=epsg:2154")
## découpage de la couche des départements français en faisant l'intersection avec la Bbox
DepIDF <- st_intersection(Dep, st_geometry(bboxRegIDF))
  
# Sélection communale
comIDF <- com %>%
  filter(CODE_REG == "11") %>%
  st_transform(.,crs = "+init=epsg:2154")
  # Obtenir les communes

##  Extraction grilles (sélection des entités par intersection des couches)
grid1k<- st_transform(grid1k,crs = "+init=epsg:2154")
grid1k$IN_IDF <- st_intersects(x = grid1k, y = RegIDF, sparse = FALSE)
grid1kIDF <- grid1k[grid1k$IN_IDF == "TRUE",]

grid200<- st_transform(grid200,crs = "+init=epsg:2154")
grid200$IN_IDF <- st_intersects(x = grid200, y = RegIDF, sparse = FALSE)
grid200IDF <- grid200[grid200$IN_IDF == "TRUE",]

## Limites métropole
MGP <- comMgp %>%
summarise() %>%
 st_transform(.,crs = "+init=epsg:2154")

Génération des fonds de carte métropole du Grand Paris

Une fois que les principales géométries sont préparées, la prochaine étape consiste à créer un fonds de carte centré sur la métropole du Grand Paris (MGP) et d’en faire un carton cartographique qui s’ajouterait aux cartes produites à l’échelle de la région Île-de-France, afin de faciliter la lecture des cartes pour la zone correspondant au centre urbain métropolitain.

# Créer un carton de zoom cartographique sur le Grand Paris
# Niveau du zoom et déplacement du carton dans la bbox, pour l'ensemble des couches.

## Création et jointure du carton à la couche originale 
bgMapWZoom <- function(sf){
  
  # Bbox autour de la métropole
  bb <- st_bbox(MGP)
  bboxMGP <- st_as_sfc(st_bbox(bb + c(-1000,-1000,1000,1000), crs = 2154))
  bboxMGP <- st_transform(bboxMGP, crs = "+init=epsg:2154")
  
  # Limites du sf : Placement sur les coordonnées et agrandissement * 1.5
  zoom <- move_and_resize(x = sf, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5)
  
  result <- inset_rbinder(list(sf, zoom))
  
  return(result)
  
}

MGP_box <- bgMapWZoom(sf = MGP)
DepMgp_box <- bgMapWZoom(sf = DepIDF)
ComMgp_box <- bgMapWZoom(sf = comIDF)
IrisMgp_box <- bgMapWZoom(sf = irisIDF)
Grid200Mgp_box <- bgMapWZoom(sf = grid200IDF)
Grid1kMgp_box <- bgMapWZoom(sf = grid1kIDF)

Les données Open Street Map pour l’habillage cartographique

Tous les éléments qui habillent les cartes sont issus des données Open Street Map. Ces données Open Street Map sont directement importées sur R via l’API Overpass, facilitant le travail de reproductibilité. Un lien renvoie l’utilisateur à la documentation WIki de ces données importées, afin de prendre connaissance des métadonnées.

Dans un premier temps on récupère les données et on transforme les objets souhaités en format Simple Feature (sf), système de coordonnées en Lambert 93, système de coordonnées du projet. Dans un second temps, les données géométriques préparées sont exportées dans l’objet au format geopackage “geom_IDF.gpkg”, crée précédemment. Dans un troisième temps, des éléments d’habillage sont ajoutés au carton sur la Métropole du Grand Paris, les données sont exportées dans l’objet au format geopackage “geom_WithBoxMgp.gpkg”, créé précédemment.

#1 Import et mise en forme es données OSM 


# available_features() : Parcourir les thèmes des objets géographiques disponibles et  sélectionner un thème pour parcourir ses tags : available_tags ("")
# Mise en garde : ! beaucoup de thèmes sont vides ; Ne pas hésiter à vérifier les infos sur le Wiki OpenStreetMap ; ! Si erreur relative à la mémoire, voir issue https://github.com/ropensci/osmdata/issues/97 et https://github.com/ropensci/osmdata/issues/147
#  Astuce : ! recourir à un VPN permet de ne pas être contraint pas des quotas de téléchargement pour chaque adresse IP
# Augmenter la mémoire
memory.limit(size = 20000)

## Principaux cours d'eau

###La Seine

Seine <- getbb("Île-De-France") %>%
      opq(timeout = 50) %>% # ! Augmenter le temps de transmission des données (ex : de 25 (default) à 50)
       add_osm_feature(key = 'name', value = "La Seine") %>%
  osmdata_sf() # Voir doc : https://wiki.openstreetmap.org/wiki/France/Cours_d%27eau

Seine <- Seine$osm_multilines # Transformation sf

Seine <- Seine %>%
  st_transform(., crs = "+init=epsg:2154") %>% # Système de coordonnées en lambert 93
  st_intersection(., st_geometry(bboxRegIDF)) # Garder uniquement les informations contenues dans la Bbox de la région

### La Marne

Marne <- getbb("Île-De-France") %>%
      opq(timeout = 50) %>%
       add_osm_feature(key = 'name', value = "La Marne")%>%
  osmdata_sf() # Voir doc : https://wiki.openstreetmap.org/wiki/France/Cours_d%27eau

Marne <- Marne$osm_multilines

Marne <- Marne %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF))

## Réseaux de communication

### Réseau routier

Autoroutes <-getbb("Île-De-France") %>%
      opq() %>%
      add_osm_feature(key = 'highway', value = "motorway")%>%
  osmdata_sf()
# Voir doc : https://wiki.openstreetmap.org/wiki/FR:France_roads_tagging

Autoroutes <- Autoroutes$osm_lines

Autoroutes <- Autoroutes %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) 

# ! Astuce : enlever les objets géographiques qui viendraient se superposer à la bbox :
Autoroutes <- Autoroutes %>%
  st_difference(., zoom_Dep %>%
                   summarise()) %>%
    filter(osm_id %in% .$osm_id)

## Réseau ferroviaire

#### Grandes lignes

# Import et passage en sf
Trains <-getbb("Île-De-France") %>%
      opq(timeout = 50) %>%
       add_osm_feature(key = 'railway', value = "rail")%>%
  osmdata_sf() # Voir doc : https://wiki.openstreetmap.org/wiki/FR:Key:railway

Trains <- Trains$osm_lines

# Système de projection, intersection et filtrage sur attribut
TrainsGL <- Trains %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) %>%
  filter(importance == "international" | importance == "national")

# ! Idem, enlever les objets géographiques qui viendraient se superposer à la bbox :

TrainsGL <- TrainsGL %>%
  st_difference(., zoom_Dep %>%
                   summarise()) %>%
    filter(osm_id %in% .$osm_id)


#### Lignes régionales (RER)

TrainsRegion <- Trains %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) %>%
  filter(importance == "regional")


TrainsRegion <- TrainsRegion %>%
  st_difference(., zoom_Dep %>%
                   summarise()) %>%
    filter(osm_id %in% .$osm_id)

### Aéroports

Aeroports <-getbb("Île-De-France") %>%
      opq(timeout = 50) %>%
       add_osm_feature(key = 'aeroway', value = "aerodrome")%>%
  osmdata_sf() # Voir doc : https://wiki.openstreetmap.org/wiki/Key:aeroway

Aeroports <- Aeroports$osm_polygons

Aeroports <- Aeroports %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) %>%
  filter(aerodrome == "international") %>%
   select(osm_id, name, geometry) 

## Eléments forestiers

Forets <-getbb("Île-De-France") %>%
      opq(timeout = 100) %>%
       add_osm_feature(key = 'landuse', value = "forest")%>%
  osmdata_sf(stringsAsFactors = T)  #  Voir doc : https://wiki.openstreetmap.org/wiki/FR:Tag:landuse%3Dforest

Forets <- Forets$osm_polygons
Forets$area <- as.numeric(st_area(Forets))# Calcul de l'aire de chaque polygone pour sélection

Forets <- Forets %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) %>%
  filter(rank(desc(area)) <= 25) %>% # On sélectionne les plus grandes forets (polygones)
  select (osm_id, name, geometry)

## Villes
### Principales villes (population >100 000)

GrdeVilles <-getbb("Île-De-France") %>%
      opq(timeout = 100) %>%
       add_osm_feature(key = 'place', value = "city")%>%
  osmdata_sf(stringsAsFactors = T)  #  Voir doc : https://wiki.openstreetmap.org/wiki/FR:Key:place

GrdeVilles <- GrdeVilles$osm_points

GrdeVilles <- GrdeVilles %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF)) %>%
  select(osm_id, name, geometry)

### Villes moyennes (population [10 000 ; 100 000])

MoyVilles <-getbb("Île-De-France") %>%
      opq() %>%
       add_osm_feature(key = 'place', value = "town")%>%
  osmdata_sf()

MoyVilles <- MoyVilles$osm_points

MoyVilles <- MoyVilles %>%
  st_transform(., crs = "+init=epsg:2154") %>%
  st_intersection(., st_geometry(bboxRegIDF))

#### Choisir les deux principales villes pour les départements de la grande couronne

MoyVilles$population <-as.numeric(MoyVilles$population)

MoyVilles$DEP <- substr(MoyVilles$addr.postcode,1,2)

# ! La sélection sur le rang ne fonctionne pas pour le département 77 informations lacunaires), on sélectionne la ville de Meaux et de Fontainebleau manuellement
Villes77 <- MoyVilles %>%
   ungroup %>%
  group_by(DEP) %>%
  filter (name == "Meaux" | name == "Fontainebleau")

Villes <- MoyVilles %>%
  group_by(DEP) %>%
  filter (DEP== "78" | DEP== "91" | DEP== "95", rank(desc(population)) <= 2) %>%
  ungroup(.)%>%
  select(osm_id, name, geometry, -DEP) %>%
  rbind(GrdeVilles,.) 

#Variables sur les coordonnées X,Y pour affichage des labels (noms des principales communes)
coords <- st_coordinates(Villes)
Villes$X <- coords[,1]
Villes$Y <- coords[,2]

# Recodage quand écriture utf-8 ne passe pas, exemple avec la commune d'Evry
Villes <- Villes %>% mutate(name =recode(name, 
                         "Évry" = "Evry"))


# 3 ajout d'éléments sur le carton

##  Seine et Marne
# zoom_seine <- move_and_resize(x = Seine, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 
zoom_seine <- bgMapWZoom(sf = Seine)
zoom_marne <- bgMapWZoom(sf = Marne)
# zoom_marne <- move_and_resize(x = Marne, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 
##
# SeineMgp <- inset_rbinder(list(Seine, zoom_seine))

# MarneMgp <- inset_rbinder(list(Marne, zoom_marne))

##  Aeroport
zoom_Aeroports <- move_and_resize(x = Aeroports, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 

AeroportsMgp <- inset_rbinder(list(Aeroports, zoom_Aeroports))

## Réseau routier 
zoom_Autoroutes <- move_and_resize(x = Autoroutes, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 

AutoroutesMgp <- inset_rbinder(list(Autoroutes, zoom_Autoroutes))

## Réseau ferroviaire
### Longue portée
zoom_TrainsGL <- move_and_resize(x = TrainsGL, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 

TrainsGLMgp <- inset_rbinder(list(TrainsGL, zoom_TrainsGL))

### Portée régionale
zoom_TrainsRegion <- move_and_resize(x = TrainsRegion, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 

TrainsRegionMgp <- inset_rbinder(list(TrainsRegion, zoom_TrainsRegion))

## Villes
zoom_villes <- move_and_resize(x = Villes, mask = bboxMGP, xy = c(547000, 6774000), k = 1.5) 

VillesMgp <- inset_rbinder(list(Villes, zoom_villes))

Production cartographique densité de transactions

# Cartographie - exemple : Creation d'une carte de densité avec linemap

## Exemple sur les indicateurs issus de BIEN
DensityPrices <- cassmir_Communes %>%
  group_by(INSEE_COM)%>%
  summarise (Total_transactions = sum(B_AC_TOT, na.rm = T)) %>%
  select( INSEE_COM, Total_transactions)


DensityPrices <- getgrid(x = DensityPrices, cellsize = 1000, var = "Total_transactions" )

# Format et emplacement de stockage de l'image
png(file = "fig/DensityTransacsBIENMap.png", width = 2000, height = 1250, res = 150)

opar <- par(mar=c(0,0,0,0), bg = "ivory2")

plot(st_geometry(RegIDF), col="ivory1", border = NA)

linemap(x = DensityPrices, var = "Total_transactions", k = 4, threshold = 40,
        col = "ivory1", border = "ivory4", lwd = 0.6, add = TRUE)

mtext(text = "Les transactions immobilières\nen Ile-de-France (1996-2018)", side = 3, line = -4, 
      col = "black", font = 4, adj = 0.05, cex = 2 )
mtext(text = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN - ANR WISDHOM - UMR Géographie-Cités, UMS RIATE, 2020", side= 1, 
      line = - 1, col = "black", font = 3, adj = 0.98, cex = 0.8 )

par(opar)

dev.off()

## graphique sur le nombre d'observations annuelles

## Theme
theme_tmd <- function(){
  tmd <-  theme(axis.text.x = element_text(size = 10, angle= 45,),
                axis.text.y = element_text(size = 10),
                legend.title = element_text(size = 0, colour = NA),
                legend.position = 'bottom',
                legend.direction = "horizontal",
                legend.background = element_rect(fill = NA),
                panel.background = element_rect(fill = "white"), 
                plot.background = element_rect(fill = "gray92"),
                panel.grid.major = element_line(colour = "gray92", size = 1),
                legend.text = element_text(size = 12),
                plot.margin = unit(x = c(0.25, 0.25, 0.25, 0.5), units = "cm"),
                plot.caption= element_text(size = 10),
                legend.box.spacing = unit(x = 0.25, units = "cm"),
                legend.margin=margin(t = 0, r = 0, b = -0.2, l = 0,unit = "cm"),
                plot.title=element_text(size =12,face="bold",hjust = 0.5, vjust = 0.5),
                axis.title = element_text(size =12),
                plot.subtitle = element_text(size =12))
}


TotPopBIEN <- as.data.frame(cassmir_Communes) %>%
  group_by(annee)%>%
  summarise (Total_transactions = sum(B_AC_TOT, na.rm = T)) %>%
  select( annee, Total_transactions)
TotPopBIEN$Base <- "Observations issues de la base BIEN"

TotPopPTZ <- as.data.frame(cassmir_Communes) %>%
  group_by(annee)%>%
  summarise (Total_transactions = sum(P_AC_TOT, na.rm = T)) %>%
  select( annee, Total_transactions)
TotPopPTZ$Base <- "Observations issues de la base PTZ"

TotPop <- rbind(TotPopBIEN, TotPopPTZ)

png(file = "fig/NobservationsPlot.png", width = 250, height = 150, units = "mm", res = 150)
options(scipen=999)
TotPop%>%
  ggplot(., aes(x=annee, y=Total_transactions)) +
   geom_line()+
  facet_wrap(~Base,scales = "free")+
   scale_x_continuous(breaks = c(1996,1997, 1998, 1999,2000, 2001,2002, 2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016, 2017,2018))+
  theme_tmd() +
  labs(title = "Nombre annuel d'observations pour chaque échantillon préparé" , x= "Année" , y= "Nombre d'obsevations")+ 
  labs(caption = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN et PTZ, 2020") 

dev.off()

Production cartographique et graphique sur les prix

# Cartographie : cartographie choroplète à la commune sur les données des prix en ajoutant un carton pour un zoom localisé sur la MGP. 

## Creation d'une fonction générique adaptée à la production d'une planche de cartes choroplètes avec un carton
createFacetsMapWithBox <- function(.data, varname, Colors, cuts, titleLeg,FacetsVar, NfacetCols, Mask, maptitle){
  
  tm_shape(.data, projection = 2154, unit = "km") +
    tm_polygons(varname,
                breaks = cuts,
                palette = Colors, 
                border.col = "white", 
                border.alpha = 0.5,
                title = titleLeg) +
    tm_layout(frame = FALSE)+
    tm_facets(by = FacetsVar, 
              drop.NA.facets = TRUE,
              free.coords = FALSE, ncol = NfacetCols)+
    tm_layout(main.title = maptitle,
              main.title.size = 1.1,
              main.title.position = c("left", "top"),
              between.margin = 0,
              legend.outside = TRUE,
              legend.position = c("left","bottom")) +
    tm_scale_bar(color.dark = "gray60",
                 position = c(0.8, 0.1),
                 breaks = c(0,10,20,30),
                 lwd= 0.15, 
                 text.size = 0.6) + 
    tm_compass(type = "arrow", size = 1, text.size = 0.5,
               color.dark = "gray60", text.color = "gray60",
               position = c(0.01,0.77)) +
    tm_shape(Mask) +
    tm_borders(col = NA, alpha=0.1) +
    tm_layout(frame = FALSE)+
    tm_shape(Marne) +
    tm_lines (col = "blue", 
              lwd = 1, 
              lty = "solid", 
              alpha = NA) +
    tm_layout(frame = FALSE)+
    tm_shape(Seine) +
    tm_lines (col = "blue", 
              lwd = 1, 
              lty = "solid", 
              alpha = NA) +
    tm_layout(frame = FALSE)+
    tm_shape(DepMgp_box) +
    tm_borders(col = "black", 
               alpha=0.8)+
    tm_layout(frame = FALSE) +
    tm_credits("CC - Base de données CASSMIR, indicateur construit à partir des données BIEN - 2020", position=c("left", "bottom"))
  
}



## Sélection de la variable dans la couche commune de la base. On sélectionne uniquement 4 années pour faciliter la visualisation
MeanPricesBox<- as.data.frame(cassmir_Communes) %>%
  filter (B_PX_M >=100, annee == 1999| annee == 2007 | annee == 2012 | annee == 2018) %>%
  select( INSEE_COM, annee, B_PX_M, -geom)

## Discrétisation
Breaks <- getBreaks(MeanPricesBox$B_PX_M, nclass = 10, method = "quantile")

## Couleurs
Col <-  carto.pal(pal1 = "blue.pal", pal2 = "red.pal",  n1 = 5, n2 = 5)

## Jointure avec les données géométriques "vierges"
MeanPricesBox$INSEE_COM<-as.character(MeanPricesBox$INSEE_COM)
MeanPricesBox <- full_join(ComMgp_box, MeanPricesBox, by = "INSEE_COM")   

## Format et emplacement de stockage de l'image
png(file = "fig/4YearsfacetPricesMap.png", width = 250, height = 150, units = "mm", res = 150)

## Exécution de la fonction
createFacetsMapWithBox(.data = MeanPricesBox, Colors = Col,
                       varname = "B_PX_M", 
                       cuts = Breaks,
                       Mask = ComMgp_box,
                       FacetsVar = "annee",
                       NfacetCols = 2,
                       maptitle = "Dynamique des prix immobiliers en île-de-France",
                       titleLeg = "Prix moyen à la commune (quantile)")

dev.off()


# Production graphique

png(file = "fig/facetMedianPricesGroupesPopPlot.png", width = 250, height = 150, units = "mm", res = 150)
options(scipen=999)
cassmir_GroupesPop%>%
  group_by(Groupes,annee,Parties)%>%
  ggplot(., aes(x=annee, y=B_PX_APP_Q2, fill=Groupes )) +
   geom_line()+
  geom_path(aes(color = Groupes))+
  facet_grid(TypeGroupe~Parties,scales = "free")+
  theme_tmd() +
  labs(title = "Prix d'achat et de vente annuel médian des acquéreurs et vendeurs pour chaque groupe de population" , x= "Année" , y= "Prix (en euros)")+ 
  labs(caption = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN - 2020")

dev.off()


# Regression plot
# cassmir_Communes$INSEE_COM<- as.numeric(as.character(cassmir_Communes$INSEE_COM))
Cassmir_AllForFacetReg<- bind_rows(
  cassmir_Communes%>% 
    filter (annee == 1999 | annee == 2018) %>%
    select(ID = INSEE_COM, annee,B_PM_APP_M) %>%
    spread(annee, B_PM_APP_M)%>%
    mutate(Groupe = "Communes")  ,
   cassmir_Grid1km%>% 
    filter (annee == 1999 | annee == 2018) %>%
    select(ID = Carreau_ID, annee,B_PM_APP_M) %>%
    spread(annee, B_PM_APP_M)%>%
    mutate(Groupe = "Carreaux1km") ,
   cassmir_Grid200m%>% 
    filter (annee == 1999 | annee == 2018) %>%
    select(ID = Carreau_ID, annee,B_PM_APP_M) %>%
    spread(annee, B_PM_APP_M)%>%
    mutate(Groupe = "Carreaux200m"))

options(scipen = 1)
options(digits = 2)
png(file = "fig/AllForFacetRegPrices.png", width = 250, height = 150, units = "mm", res = 150)
# The plot
ggplot(data = Cassmir_AllForFacetReg, aes(`2018`, `1999` )) +
  # Plot the points
  geom_point() +
  # Add the regression line (without S.E. shading)
  geom_smooth(method = "lm") +
  # Change the theme
  theme_tmd() +
  # Stats
    stat_cor(method = "pearson")+
  # Plot by group
  facet_wrap(~Groupe, nrow = 2, ncol = 2) +
  ggtitle("Evolution de la hiérarchie spatiale des prix moyens au m² des appartements entre 1999 et 2018",
    subtitle = "Régression linéaire et coefficient de Pearson sur les prix des 3 agrégats spatiaux de CASSMIR") + 
  labs(caption = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN - 2020",
        x= "Prix moyen en 2018" , y= "Prix moyen en 1999")

dev.off()

Production cartographique et graphique sur les acquéreus-vendeurs

# Cartographie -  Combinaison d'une carte mixte à partir de la réalisation de 2 cartes simples sur les acquéreurs : l'une choro, l'autre cercle proportionnelle. Objectif : pour l'exemple on crée une carte sur le pourcentage d'ouvriers acquéreurs et le nombre total d'ouvriers acquéreurs en 2012.


## Sélection des variables dans la couche commune de la base. 
OuvDoubleMapBox <- as.data.frame(cassmir_Communes) %>%
  filter (annee == 2012) %>%
  select(INSEE_COM, annee, B_AC_TOT,  B_AC_PPH_PCS6, -geom) %>%
  mutate(N_Ouvriers = (B_AC_TOT*B_AC_PPH_PCS6)/100)

## Discrétisation
Breaks <- getBreaks(OuvDoubleMapBox$B_AC_PPH_PCS6, nclass = 4, method = "quantile")

##  Couleurs
Col <-  carto.pal(pal1 = "red.pal", n1 = 4)

## Jointure avec les données géométriques "vierges"
OuvDoubleMapBox$INSEE_COM<-as.character(OuvDoubleMapBox$INSEE_COM)
OuvDoubleMapBox <- full_join(ComMgp_box, OuvDoubleMapBox, by = "INSEE_COM")   

## Map choro
MapPourcOuv2012<- tm_shape(OuvDoubleMapBox, projection = 2154) +
  tm_polygons("B_AC_PPH_PCS6",
              breaks = Breaks,
              palette = Col, 
              border.col = NA, 
              border.alpha = 0,
              legend.reverse = TRUE,
              title ="Pourcentage d'acquéreurs ouvrier") +
  tm_shape(ComMgp_box) +
  tm_borders(col = NA, alpha = 0.1) +
  tm_shape(DepMgp_box %>% filter(CODE_DEPT %in% c("75", "77", "78", "91", "92", "93", "94", "95"))) +
  tm_borders(col = "gray20",
             lwd = 1.5,
             alpha = 0.8)+
  tm_layout(frame = FALSE) +
  tm_scale_bar(color.dark = "gray60",
               position = c(0.8, 0.1),
               breaks = c(0,10,20,30),
               lwd = 0.15, 
               text.size = 0.6) +
  tm_credits(text = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN, 2020", position = c(0.3, 0.02), size = 0.5)

## Map stock
MapStockOuv2012 <- tm_shape(OuvDoubleMapBox, projection = 2154) + 
  tm_bubbles(size = "N_Ouvriers",
             alpha = 0.8, col = "white", 
             border.lwd = 0.8,
             title.size = "Nombre d'acquéreurs ouvrier",
            style= "cont") +
  tm_layout(frame = FALSE)


##  Format et emplacement de stockage de l'image
png(file = "fig/WorkersDoubleInsideMap_rv.png", width = 250, height = 150, units = "mm", res = 300)
##  visualiser toutes les infos sur une seule carte 
MapPourcOuv2012 + MapStockOuv2012
dev.off()


# Production graphique
# str(cassmir_GroupesPop, list.len=ncol(cassmir_GroupesPop))
Contingence_ac_vendeur<- cassmir_GroupesPop%>% 
  filter(TypeGroupe == "Generationnel", Parties == "Acquereurs",  annee == 2012) %>%
  select (Groupes, B_AC_ACHA_AGE1, B_AC_ACHA_AGE2, B_AC_ACHA_AGE3,B_AC_ACHA_PCS7)

row.names(Contingence_ac_vendeur) <- Contingence_ac_vendeur$Groupes
Contingence_ac_vendeur <-  Contingence_ac_vendeur%>%
  select (-Groupes) %>%
 rename(`Personne active entre 18 et 30 ans` = `B_AC_ACHA_AGE1`,
        `Personne active entre 30 et 50 ans` = `B_AC_ACHA_AGE2`,
        `Personne active de plus de 50 ans` =  `B_AC_ACHA_AGE3`,
       `Personne à la retraite` = B_AC_ACHA_PCS7)

Khi2test<-chisq.test(data.matrix(Contingence_ac_vendeur))
residus_test_Khi2<-Khi2test$residuals

residus_test_Khi2<-melt(residus_test_Khi2, varnames = c("Acquereurs", "Vendeurs"),value.name = "residus")
Contingence_ac_vendeur<- Contingence_ac_vendeur %>%
  group_by(rownames(.))%>%
  gather("Vendeurs", "effectif", c(1:4)) %>%
  rename("Acquereurs" =`rownames(.)`)%>%
left_join(., residus_test_Khi2, by = c("Acquereurs", "Vendeurs"))

#plot
 TablePlotAcqVenGeneration <- ggplot(Contingence_ac_vendeur, aes(Acquereurs,Vendeurs, fill= residus)) +
  geom_tile()+
  geom_text(aes(label=round(effectif, digits = 0))) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-20,20), space = "Lab",
                       name= "Variation des résidus" )+
    # theme_tmd() +
   theme(axis.text.x = element_text(size = 10, angle= 10,),
         axis.text.y = element_text(size = 10))+
  labs(title = "Les relations acquéreurs-vendeurs par groupes d'âges")+ 
  labs(subtitle = "Heat map sur les résidus du test du Khi2 de Pearson.\nLes valeurs répertorient le nombre d'échanges entre parties de la vente.\nLes couleurs représentent les résidus du Khi².\nLecture: En 2012, 619 logements ont été achetés par une population âgée entre 18 et 30 ans\nà une population dans la même tranche d'âge.\nCette valeur est au-dessus de la valeur théorique en cas d'indépendance statistique des relations.",
       x= "Vendeurs", y= "Acquéreurs")+
   scale_x_discrete(labels=c("AGE1" = "Personne active entre 18 et 30 ans", 
                             "AGE2" = "Personne active entre 30 et 50 ans",
                              "AGE3" = "Personne active de plus de 50 ans",
                             "PCS7" = "Personne à la retraite"))+
  labs(caption = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN, 2020")
 
 ggsave("TablePlotAcqVenGeneration.png",plot= TablePlotAcqVenGeneration, path = "fig" ,  device = "png", width = 250, height = 150, units = "mm", dpi = 150)

Production cartographique et graphique sur les régimes d’achat et types de mutation

# Production graphique

png(file = "fig/PlotDurCredSocial.png", width = 250, height = 150, units = "mm", res = 150)
options(scipen=999)
cassmir_GroupesPop%>%
  filter(Parties == "Acquereurs", TypeGroupe== "Social", Groupes == "PCS3" | Groupes == "PCS4" | Groupes == "PCS5" | Groupes == "PCS6") %>%
  select(Groupes, annee, P_AC_PPH_DURPRET_Q1, P_AC_PPH_DURPRET_M, P_AC_PPH_DURPRET_Q2, P_AC_PPH_DURPRET_Q3 )%>%
   rename(`Premier quartile` = `P_AC_PPH_DURPRET_Q1`,
        `Moyenne` = `P_AC_PPH_DURPRET_M`,
        `Médiane` =  `P_AC_PPH_DURPRET_Q2`,
       `Troisième quartile` = P_AC_PPH_DURPRET_Q3)%>%
  gather("Indicateurs", "Value", 3:6 ) %>%
  mutate(Groupes = case_when(
    Groupes == "PCS3" ~ "Cadres et professions intellectuelles supérieures",
         Groupes == "PCS4" ~ "Professions intermédiaires",
         Groupes == "PCS5" ~ "Employés",
         Groupes == "PCS6" ~ "Ouvriers"))%>%
   group_by(Groupes,annee)%>%
  ggplot(., aes(x=annee, y=Value, fill=Indicateurs )) +
   geom_line()+
  geom_path(aes(color = Indicateurs)) +
   facet_wrap(~ Groupes,scales = "fix")+
  scale_x_continuous(breaks = c(1996,1997, 1998, 1999,2000, 2001,2002, 2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016)) +
   theme_tmd() +
  labs(title = "Evolution de la durée de remboursement du crédit pour les acquéreurs des principales catégories sociales" , x= "Année" , y= "Durée (en mois)")+ 
  labs(subtitle = "Il s'agit de la durée de remboursement du crédit principal pour des primo-accédants aidés par un Prêt à Taux Zéro.",
       x= "Année de l'emprunt", y= "Nombre de mois pour le remboursement du crédit")
  labs(caption = "CC - Base de données CASSMIR, indicateur construit à partir des données PTZ, 2020")

dev.off()

Production cartographique et graphique sur les types de biens immobiliers

# Cartographie -  Carte en cercles proportionnels

## Sélection des variables dans la couche commune de la base. 
LogMapBox<- as.data.frame(cassmir_Communes) %>%
  filter(B_IMMO_APP_ANC2 > 0) %>%
  select(INSEE_COM, B_IMMO_APP_ANC2, -geom) %>%
  group_by(INSEE_COM) %>%
  summarise(Total_transactions = sum(B_IMMO_APP_ANC2, na.rm = TRUE)) 

## Discrétisation
Breaks <- getBreaks(LogMapBox$Total_transactions, nclass = 5, method = "quantile")

##  Couleurs
Col <-  carto.pal(pal1 = "red.pal", n1 = 4)

## Jointure avec les données géométriques "vierges"
LogMapBox <- full_join(ComMgp_box, LogMapBox, by = "INSEE_COM")   

## Carte
MapStockOldH <- tm_shape(LogMapBox, projection = 2154) +
  tm_polygons(col = "gray95", border.col = "gray80") +
  tm_bubbles(size = "Total_transactions",
             col = "gray40", alpha = 0.7,
             title.size = "Nombre de transactions") +
  tm_shape(DepMgp_box %>% filter(CODE_DEPT %in% c("75", "77", "78", "91", "92", "93", "94", "95"))) +
  tm_borders(col = "gray20",
             lwd = 1.5,
             alpha=0.8) +
  tm_scale_bar(color.dark = "grey60",
               position = c(0.8, 0.1),
               breaks = c(0,10,20,30),
               lwd= 0.15, 
               text.size = 0.6) + 
  tm_layout(frame = FALSE) +
  tm_credits(text = "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN\nANR WISDHOM - UMR Géographie-Cités, UMS RIATE, 2020", position = c(0.3, 0.02), size = 0.5)

## sauvegarde
png(file = "fig/CartoOlddHouses_rv.png", width = 250, height = 150, units = "mm", res = 300)
MapStockOldH
dev.off()

Production cartographique multi-scalaire

# Cartographie -  Cartes à trois échelles différentes sur le même indicateur: Mettre sur une planche commune les trois cartes avec représentations lissées pour carroyage 1km et 200m.

## Exemple sur le prix moyen pour l'année 2012
year <- 2012
var <- "B_PX_APP_M"
cassmir_Communes$B_PX_APP_M <-as.numeric(cassmir_Communes$B_PX_APP_M)
cassmir_Grid1km$B_PX_APP_M <-as.numeric(cassmir_Grid1km$B_PX_APP_M)
cassmir_Grid200m$B_PX_APP_M <-as.numeric(cassmir_Grid200m$B_PX_APP_M)

titleLeg <- "Prix moyen en 2012\n(méthode déciles)"
mapSource <- "CC - Base de données CASSMIR, indicateur construit à partir des données BIEN\nANR WISDHOM - UMR Géographie-Cités, UMS RIATE, 2020"

## Discrétisation sur les 3 couches pour bornes communes
Bv <- rbind (cassmir_Communes %>%
  filter(annee == year) %>%
  select(var),
  cassmir_Grid1km %>%
  filter(annee == year) %>%
  select(var),
   cassmir_Grid200m %>%
  filter(annee == year) %>%
  select(var))

Bv <- getBreaks(Bv$B_PX_APP_M, nclass = 10, method = "quantile")


## Couches d'habillage
seineBox <- st_read("CASSMIR_Geom/geom_WithBoxMgp.gpkg", quiet = TRUE, layer = "Seine")
marneBox <- st_read("CASSMIR_Geom/geom_WithBoxMgp.gpkg", quiet = TRUE, layer = "Marne")
VillesBox <- st_read("CASSMIR_Geom/geom_WithBoxMgp.gpkg", quiet = TRUE, layer = "Villes")

## Ne garder que les geom au seine de l'IDF
require(lwgeo)
idf <- DepIDF %>% filter(CODE_DEPT %in% c("75", "77", "78", "91", "92", "93", "94", "95"))

# site_snap <- st_snap(seineBox, idf, tol = 1e-9)
# seine <- st_intersection(site_snap, idf)
# site_snap <- st_snap(marneBox, idf, tol = 1e-9)
# marne <- st_intersection(site_snap, idf)
# 
# villes <- VillesBox %>% 
#   filter(!duplicated(name) & name != "Paris") %>% 
#   mutate(name = case_when(name == "Boulogne-Billancourt" ~ "Boulogne",
#                           name == "Mantes-la-Ville" ~ "Mantes",
#                           TRUE ~ name))

## Fonction : lissage
lissage <- function(sf, res, data){
  
  # Création d'un masque avec l'objet de type sf
  Mask <- st_union(sf)
  Mask <- st_sf(id = 1, geometry = Mask)
  
  # Simplification géométrique du masque
  Mask <- Mask %>%
    st_buffer(dist = res) %>%
    st_buffer(dist = -(res*2)) %>%
    st_simplify(preserveTopology = FALSE, dTolerance = res) %>%
    st_buffer(res)
  
  # rasterisation
  require(stars) 
  x <- data
  ## Résolution
  x_ras <- st_rasterize(x[var], dx = res, dy = res, 
                      xlim = st_bbox(x)[c(1,3)] + c(-res, res)/2, 
                      ylim = st_bbox(x)[c(2,4)] + c(-res, res)/2)
  st_raster_type(x_ras) ### [1] "regular"

  # construction du contour
  options(scipen = 999)
  x_contour <- st_contour(x_ras, breaks = Bv)
  #level = 1:length(row.names(x_contour))
  
  return(x_contour)
  
}

## Fonction : mapping
rastMap <- function(rast, titleMap){
  
  # Le choix de la couleur
    Col <-  carto.pal(pal1 = "green.pal", n1 = 5, pal2 = "red.pal", n2 = 5)
  #Si problème de couleurs :ajouter la derniere couleur à la lsite
  # Col <-  c("#17692C", "#287A22", "#5D9D52", "#8DBC80", "#B8D9A9", "#FCDACA", "#F4988D", "#EC4E49", "#DA0001", "#7C000C", "#7C000C")
  
  
  # Stock bbox
  bbx <- st_bbox(idf)
  
  # variable to plot
  level <-  1:length(row.names(rast))
  
  ## ggplot
  map <- ggplot() + 
  geom_sf(data = rast, aes(fill = as.factor(level)), 
                   color = Col, show.legend = FALSE)+
  scale_color_manual(values = Col) +
  scale_fill_manual(values = Col) +
  # geom_sf(data = marne, fill = NA, color = "#33f4f4", size = 0.2) +
  # geom_sf(data = seine, fill = NA, color = "#33f4f4", size = 0.2) +
  geom_sf(data = idf, fill = NA, color = "gray20", size = 0.3) +
  ## Facultatif : libellés 
  # geom_sf_text(data = villes %>% filter(name != "Boulogne"), aes(label = name), 
  #              size = 1.4, color = "grey20", fontface = 2,
  #              check_overlap = TRUE) +
  # geom_sf_label(data = villes %>% filter(name != "Boulogne"), aes(label = name), size = 1, color = "grey20",
  #              check_overlap = FALSE, label.padding = unit(0.05, "lines"), label.size = 0.1) +
  coord_sf(xlim = c(bbx[1], bbx[3]),
           ylim =  c(bbx[2], bbx[4]),
           datum = sf::st_crs(2154), expand = TRUE) +
  labs(title = titleMap) +
  theme_void()
  
  return(map)

}


## Carte grille 1km 
rast1 <- lissage(sf = grid1kIDF, 
                 res = 1000, 
                 data = cassmir_Grid1km %>% filter(annee == year) %>% select(Id_carr1km, var))

map_rast1 <- rastMap(rast = rast1, titleMap = "Carroyage 1km lissé")

# # sauvegarde
# ggsave(map_rast1, filename ="fig/map_rast1.png", 
#        device = "png", width = 8, height = 8, units = "in",  dpi = 300)
 

## Carte grille 200m
rast2 <- lissage(sf = grid200IDF, 
                 res = 200, 
                 data = cassmir_Grid200m %>% filter(annee == year) %>% select(IdINSPIRE, var))

map_rast2 <- rastMap(rast = rast2, titleMap = "Carroyage 200m lissé")

# # sauvegarde
# ggsave(map_rast2, filename ="fig/map_rast2.png", 
#        device = "png", width = 8, height = 8, units = "in",  dpi = 300)


## Carte communale
### Le choix de la couleur
Col <-  carto.pal(pal1 = "green.pal", n1 = 5, pal2 = "red.pal", n2 = 5)
  
### Stock bbox
bbx <- st_bbox(idf)

### Discrétisation 
dataComForMap <- cassmir_Communes %>% 
  filter(annee == year) %>% 
  select(ID_GEOFLA, var = var) %>%
  mutate(typo = cut(var, breaks = Bv, dig.lab = 2, 
                     include.lowest = TRUE, right = FALSE))

### Facultatif 
# dataComForMap <- dataComForMap %>% 
#   mutate(typo2 = recode(typo,
#                        "NA" = "",
#                        "[2.7e+04,1.2e+05)" = "27500 - 124603",
#                        "[1.2e+05,1.4e+05)" = "124603 - 142156",
#                        "[1.4e+05,1.6e+05)" = "142156 - 155830",
#                        "[1.6e+05,1.7e+05)" = "155830 - 167877",
#                        "[1.7e+05,1.8e+05)" = "167877 - 178695",
#                        "[1.8e+05,1.9e+05)" = "178695 - 191644",
#                        "[1.9e+05,2.1e+05)" = "191644 - 210668",
#                        "[2.1e+05,2.5e+05)" = "210668 - 245461",
#                        "[2.5e+05,3.2e+05)" = "245461 - 315658",
#                        "[3.2e+05,1.3e+06]" = "315658 - 1309232"))

dataComForMap <- dataComForMap %>% 
  mutate(typo2 = recode(typo,
                       "NA" = "",
                       "[2.7e+04,1.2e+05)" = "27500 - 124603",
                       "[1.2e+05,1.4e+05)" = "124603 - 142156",
                       "[1.4e+05,1.5e+05)" = "142156 - 155830",
                       "[1.5e+05,1.7e+05)" = "155830 - 167877",
                       "[1.7e+05,1.8e+05)" = "167877 - 178695",
                       "[1.8e+05,1.9e+05)" = "178695 - 191644",
                       "[1.9e+05,2.1e+05)" = "191644 - 210668",
                       "[2.1e+05,2.4e+05)" = "210668 - 245461",
                       "[2.4e+05,3.1e+05)" = "245461 - 315658",
                       "[3.1e+05,1.3e+06]" = "315658 - 1309232"))



### ggplot
require(ggsn)
map_com <- ggplot() + 
  geom_sf(data = dataComForMap, aes(fill = typo2), 
          colour = NA, lwd = 0.1,  
          show.legend = TRUE) +
  scale_fill_manual(name = titleLeg, 
                    values = Col, na.value = NA) +
  # geom_sf(data = marne, fill = NA, color = "#33f4f4", size = 0.2) +
  # geom_sf(data = seine, fill = NA, color = "#33f4f4", size = 0.2) +
  geom_sf(data = idf, fill = NA, color = "gray20", size = 0.3) +
  # geom_sf(data = villes %>% st_centroid(), colour = "grey20", size = 0.5, show.legend = NA) +
  # geom_sf_label(data = villes %>% filter(name != "Boulogne"), aes(label = name), size = 1, color = "grey20",
  #              check_overlap = TRUE, label.padding = unit(0.05, "lines"), label.size = 0.1) +
  coord_sf(xlim = c(bbx[1], bbx[3]),
           ylim =  c(bbx[2], bbx[4]),
           datum = sf::st_crs(2154), expand = TRUE) +
  labs(title = "Communes") +
  labs(caption = mapSource) +
  scalebar(dataComForMap, dist = 15, dist_unit = "km", 
           st.bottom = FALSE, st.color = "grey60", st.dist = 0.03,
           st.size = 1.5, location = "bottomleft", height = 0.015,
           box.color = "grey60", box.fill = c("grey60", "white"), border.size = 0.3,
           transform = FALSE) +
  theme_void() +
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 7))


# # sauvegarde
# ggsave(map_com, filename ="fig/map_com.png", 
#        device = "png", width = 8, height = 8, units = "in",  dpi = 300)


## Planche
require(patchwork)
threePlots <- map_rast1 + map_rast2 + map_com + plot_layout(ncol = 3, nrow = 1) 

### resize legend
# https://stackoverflow.com/questions/52297978/decrease-overal-legend-size-elements-and-text
addSmallLegend <- function(myPlot, pointSize = 1, textSizeT = 8, textSize = 6, spaceLegend = 0.7) {
    myPlot &
        guides(shape = guide_legend(override.aes = list(size = pointSize)),
               color = guide_legend(override.aes = list(size = pointSize))) &
        theme(plot.title = element_text(size = 10),
              plot.caption = element_text(size = 4),
              legend.title = element_text(size = textSizeT), 
              legend.text  = element_text(size = textSize),
              legend.key.size = unit(spaceLegend, "lines"))
}

### Apply on original plot
p <- addSmallLegend(threePlots)


## sauvegarde 
ggsave(p, filename ="fig/threePlots.png", 
       device = "png", width = 9, height = 3, units = "in",  dpi = 300)












BD CASSMIR - licence CC-BY-NC