#############################################
###### Script R corrigé des exercices #######
###### TD du 13/09/2024               #######
###### R. Le Goix                     #######
#############################################

## Créer un projet et son environnement de travail ----

# Contenant le chemin d'accès aux fichiers et les librairies / packages utilisés 
setwd("/Users/renaud/SIG/cours_m2_traitementsstatscarto/")
#getwd()

library(sf) # Pour manipuler des objets géographiques vectoriels
library(readr)


## Partie 1. ----

### Travail sur un tableau de données (dataframe) ----
#### Importer les données à partir d'une fichier texte (format CSV) ----

# Méthode simple
# Revenus_IRIS <- read.csv("public/data/Revenus_IRIS_PC.csv")
## Cette méthode nécessite de contrôler à postériori la structure des variables

# Méthode en utilisant l'interface RStudio qui permet 1. de contrôler la structure des variables que l'on importe ; 2. de copier / coller le code utilisé.
Revenus_IRIS_PC <- read_csv("public/data/Revenus_IRIS_PC.csv", 
                            col_types = cols(IRIS = col_character(), 
                                             DEC_MED21 = col_number(), DEC_GI21 = col_number()))


# Visualiser les colonnes et variables disponibles
head(Revenus_IRIS_PC)

# Sélectionner les colonnes qui nous intéressent (nettoyage des fichiers)
Revenus_IRIS_PC <- Revenus_IRIS_PC[,c("IRIS","DEP","DEC_MED21","DEC_GI21","DEC_PAUT21")]

# Nettoyer les fichiers. On a vu que des lignes comportaient des données numériques manquantes, codées NA. On va les supprimer. Il y a deux solutions

# On supprime les enregistrements qui contiennent des données manquantes
# Supprimer les lignes avec des valeurs NA
Revenus_IRIS_PC <- na.omit(Revenus_IRIS_PC)
# Il doit rester 2529 lignes

# On conserve les enregistrements qui sont complets
Revenus_IRIS_PC <- Revenus_IRIS_PC[complete.cases(Revenus_IRIS_PC), ]
# Il doit rester 2529 lignes (le monde est bien fait)

# On contrôle
head(Revenus_IRIS_PC)

# Ooops, il semblerait que la variable `DEC_PAUT21` soit toujours une chaîne de caractère. On a besoin de la transformer en variable numérique
# Ce n'est pas très grave
# Essayer puis voir ce que ça donne. On a un souci car la décimale est codée avec une virgule.
# On va manipuler la chaîne de caractère pour mettre un point à la place


Revenus_IRIS_PC$RevPatrimoine <- as.numeric(Revenus_IRIS_PC$DEC_PAUT21)

# On utilise la fonction `gsub` qui permet de remplacer des expressions régulières dans une chaîne de caractère
Revenus_IRIS_PC$RevPatrimoine <- gsub(",", ".", Revenus_IRIS_PC$DEC_PAUT21)

# On peut maintenant converir la variable dans un format numérique
Revenus_IRIS_PC$RevPatrimoine <- as.numeric(Revenus_IRIS_PC$RevPatrimoine)

# Et supprimer une variable devenue inutile
Revenus_IRIS_PC$DEC_PAUT21 <- NULL

## Partie 2. Créer une carte choroplèthe ----
### Installer le package sf et mapsf ----

# Ouvrir les fichiers .SHP

COM_PC = st_read("public/geom/COM_PC.shp")
IRIS_PC = st_read("public/geom/IRIS_PC.shp")

### Visualiser cartes communes et PC ----
# Méthode avec plot
plot(IRIS_PC)
plot(COM_PC)

# Méthode avec Mapview. (Utilise OpenStreetMap, nécessite une connexion internet)


mapview() + mapview (IRIS_PC) + mapview(COM_PC)


## Effectuer une jointure attributaire ----
IRIS_PC <- merge(
  x = IRIS_PC,             # l'objet sf
  y = Revenus_IRIS_PC,     # le data.frame contenant l'information statistique
  by.x = "CODE_IRIS",  # identifiant dans x
  by.y = "IRIS",  # identifiant dans y
  all.x = TRUE         # conserver toutes les lignes du fichier x (ici l'objet sf, ou fond de carte)
)

# Les deux objets ont bien été joints
head(IRIS_PC, 3)

## Une carte de base
library(mapsf)
mf_map(IRIS_PC, var = "DEC_MED21", type = "choro")

## Une carte plus travaillée
mf_init(x = COM_PC) # On utilise un fond de carte (ici la carte communale) pour baliser le domaine de travail
# plot population
mf_map(
  x = IRIS_PC,
  var = "DEC_MED21",
  type = "choro",
  breaks = "quantile",
  nbreaks = 5,
  pal = "Reds",
  border = "white",
  col_na = "lightgrey",
  lwd = 0.01,
  leg_pos = "left",
  leg_title = "Revenus médians 2021", leg_title_cex	= 0.5, leg_val_cex = 0.4, add = TRUE
)
mf_layout(title = "Les revenus des ménages à Paris et petite couronne en 2021",
          credits = "Sources: INSEE, 2021 ; IGN 2022. Réalisation R. Le Goix 2024",   scale = TRUE,
          arrow = TRUE, frame = TRUE)
# Et on rajoute la trame communale pour la lisibilité de la carte
mf_map(COM_PC, col=NA, add=TRUE)