---
title: "Préparation DVF"
author: "Boris Mericskay"
date: "02/07/2021"
output: html_document
---

> Ce script est adapaté pour le cours d'initiation à R du M2 Géoprisme (R. Le Goix, Université Paris Cité, Sept. 2024)

Source : https://github.com/bmericskay/Script1-DVF/blob/main/Script1.Rmd

Boris Mericskay et Florent Demoraes, « Préparer et analyser les données de "Demandes de valeurs foncières" en open data : proposition d’une méthodologie reproductible », Cybergeo: European Journal of Geography [En ligne], Cartographie, Imagerie, SIG, document 1031, mis en ligne le 28 septembre 2022, consulté le 25 septembre 2024. URL : http://journals.openedition.org/cybergeo/39583 ; DOI : https://doi.org/10.4000/cybergeo.39583

---

#  1- PRÉPARATION DES DONNÉES DVF OPENDATA
---
Cette première partie du script a comme objectif de préparer les données DVF en opendata avant de commencer les analyses (nettoyage, filtrage, agrégation,...).

Seul le package `tidyverse` dédié à la manipulation de données est nécessaire.

```{r}
library(tidyverse)
library(readxl) # Import de feuilles XLS / XLSX
library(sf)
```


Définition de l'environnement de travail

```{r setup}
knitr::opts_chunk$set(cache = FALSE)
data <- setwd("/Users/renaud/SIG/cours_m2_traitementsstatscarto")
```

---
## IMPORT DATASET
---

Lecture et assemblage des fichiers DVF
```{r }

DVF75 <- read.csv("data/DVF/dvf75.csv", encoding="UTF-8", stringsAsFactors=FALSE)
DVF92 <- read.csv("data/DVF/dvf92.csv", encoding="UTF-8", stringsAsFactors=FALSE)
DVF93 <- read.csv("data/DVF/dvf93.csv", encoding="UTF-8", stringsAsFactors=FALSE)
DVF94 <- read.csv("data/DVF/dvf94.csv", encoding="UTF-8", stringsAsFactors=FALSE)

DVF <- rbind (DVF75,DVF92,DVF93,DVF94)
rm(DVF75,DVF92,DVF93,DVF94)


```

```{r DimTableau}
# Les dimensions du tableau
dim(DVF)

# La fonction head() fournit les premières lignes`
# Visualiser les colonnes et variables disponibles
head(DVF)

```

---
## Filtre et nettoyage des données

![Les étapes de la préparation des données DVF. Source :Mericskay and Demoraes (2022) ](resources/WorkflowDataprep.png)


### Etape1> Sélection des mutations de type "Ventes" de "Maison" et "Appartement'

```{r }
etape1 <- DVF %>% filter(nature_mutation == "Vente")

etape1bis <- etape1 %>% filter(type_local == "Maison" | type_local == "Appartement")
```

### Etape2> Sélection et renommage des variables

```{r selection_var}
etape2 <- etape1bis %>% select(
    id = id_mutation, 
    disposition = numero_disposition, 
    parcelle= id_parcelle, 
    date = date_mutation, 
    nature = nature_mutation, 
    codecommune = code_commune, 
    departement = code_departement, 
    type = type_local, 
    surface = surface_reelle_bati, 
    piece = nombre_pieces_principales, 
    prix = valeur_fonciere, 
    latitude, 
    longitude)
```

### Etape3> Remplacement des cellules vides par des NA et suppression des NA

```{r NettoyageNA}
etape2[etape2==""] <- NA
etape3 <- etape2 %>% na.omit()
```

### Vérification geographique
```{r eval=FALSE, include=FALSE}
DVFtest <- st_as_sf(etape3, coords = c("longitude","latitude"), crs = 4326)
plot(st_geometry(DVFtest), pch = 3, cex = 0.01)
rm(DVFtest)
```

---
## Suppression des doublons et des mutations multiventes
### Etape4> Sélections des mutations simples

#### Regrouper les transactions selon l'ID, la surface et la vente
```{r}
unique <- etape3 %>% distinct(id, prix, surface)
nbunique <- unique %>% group_by(id) %>% summarise(nb = n())
```

#### Suppression des doublons et des mutations multiventes
```{r}
etape4 <- nbunique %>% filter(nb==1)
```


### Etape5> Jointure attributaire pour récupérer les informations de la mutation
```{r }
merge <- merge(etape4,etape3, by="id")
etape5 <- merge %>% 
  distinct(id, .keep_all = TRUE) %>% 
  select(id, date, type, nature, codecommune, prix, surface, piece, latitude, longitude)
```



#### Modification des formats des colonnes
```{r cars}
etape5$prix <- as.numeric(etape5$prix)
etape5$surface <- as.numeric(etape5$surface)
etape5$piece <- as.numeric(etape5$piece)
```

---
## Suppression des valeurs aberrantes

#### Fixer un seuil minimal des prix (premier et dernier percentile)
```{r cars}

quantile(etape5$prix,  probs = c(0.01, 0.999))

```


#### Fixer un seuil maximal des prix (histogramme)
```{r cars}
hist(etape5$prix, breaks = 200)
```


#### Créer deux jeux de données (maisons / appartements)

```{r Maisons_Apt}

Maisons <- etape5 %>% filter(type == 'Maison')
Appartement <- etape5 %>% filter (type == 'Appartement')    
```
 
 
#### Fixer un seuil maximal des surfaces (histogramme)
```{r }

hist(Maisons$surface, nclass = 500, xlim = c(0,600))
hist(Appartement$surface, nclass = 500, xlim = c(0,300))

```

### Etape 6> Sélection des bornes de prix (1er et dernier percentile) et de surface

```{r cars}
etape6 <- etape5 %>% 
  filter(between(prix, 10000, 7198250)) %>% 
  filter(case_when(type=='Appartement' ~  between(surface, 10, 200)) | 
                            case_when(type=='Maison' ~  between(surface, 10, 350)))
```

---
## Création du prix au m² et filtre des valeurs extrêmes et aberrantes 

### Etape7> Création de la variable prix/m²
```{r cars}
etape7 <- etape6 %>% mutate(prixm2 = prix/surface)
```


### Etape8 > Sélection des bornes de prix au m²

#### Fixer un seuil minimal des prix au m² (percentile)
```{r cars}
quantile(etape7$prixm2, probs = c(0.01, 0.999))
```

#### Fixer un seuil maximal des prix au m² (histogramme)
```{r cars}
hist(etape7$prixm2, breaks = 1000, xlim = c(0,20000))
```

#### Filtres des valeurs
```{r cars}
etape8 <- etape7 %>% filter(between(prixm2,330,20000))
```

#### Transformation de la date en année
```{r cars}
etape8$date <- as.character(etape8$date)
etape8 <- etape8 %>% mutate(ANNEE = substr(etape8$date, 1, 4))
```

#### Arrondir les variables numériques
```{r cars}
etape8$prix <- round(etape8$prix)
etape8$prixm2 <- round(etape8$prixm2)
```

---

## Enrichissement des mutations avec de nouvelles variables 
--------------------------------------------

### Etape9 >Ajout de variables relatives à la commune et à l'EPCI (Liste des communes du Code Officiel Géographique)

Voir Code Officiel Géographique : https://www.insee.fr/fr/information/6800675

```{r}

COG <- read.csv("data/v_commune_2023.csv", sep=",", stringsAsFactors=FALSE)
COG <- COG[, c("COM", "LIBELLE")]
COG <- COG %>% mutate(codecommune = COM)
etape9 <- merge(etape8, COG, "codecommune")


```




### Etape10 > Structuration du jeu de données final
```{r }

DVFOK <- etape9 %>% select(id, date, annee = ANNEE, type, prix, surface, prixm2, codecommune, commune = LIBELLE, latitude, longitude)

```

### Vérification geographique
```{r}
DVFtest <- st_as_sf(DVFOK, coords = c("longitude","latitude"), crs = 4326)
plot(st_geometry(DVFtest), pch = 3, cex = 0.01)
rm(DVFtest)
```


### Ecrire le jeu de données final en csv
```{r cars}
write.csv(DVFOK, 'data/DVFOK.csv')
```

---
### STOP ICI 

---
## Diagnostic des types de mutations

### Nombre de mutations

```{r cars}
length(unique(etape3$id))
```

### Nombre de mutations multiventes

```{r cars}
Mutationscomplexes <- nbunique %>% filter(nb>1)
length(unique(Mutationscomplexes$id))
```

### Nombre de mutations monoventes/monolignes

```{r cars}
uniquesimple <- etape3 %>% group_by(id) %>% summarise(prix = max(prix), surface = max(surface), nb = n())
nbtransactionstimplemonoligne <-uniquesimple %>% filter(nb==1)
length(unique(nbtransactionstimplemonoligne$id))
```
### Nombre de mutations monoventes/multilignes

```{r cars}
(length(unique(etape3$id)))-(length(unique(nbtransactionstimplemonoligne$id)))
```
### Diagnostic du nombre de biens par mutation

```{r cars}
diaggroup <- nbunique %>% mutate(as.character(nb)) %>% group_by(nb) %>% summarise(tot = n())
```

### Nombre de biens et nombre moyen de biens pour les mutations multiventes

```{r cars}
length(Mutationscomplexes$id)
sum(Mutationscomplexes$nb)
mean(Mutationscomplexes$nb)
```

---

# 2- Indicateur génériques et visualisation de données

---


## Tableau récapitulatif des indicateurs génériques

```{r cars}
recap <- DVFOK %>% group_by(type) %>% summarise(tot = n (), prixmed = median(prix), prixmoy = mean(prix), surfmed = median(surface), surfmoy = mean(surface), prixm2med = median(prixm2), prixm2moy = mean(prixm2))

recap <- recap %>% mutate(part = (tot/sum(tot)*100))
recap
```



## Visualisation de données avec ggplot2
---

### Récapitulatif des ventes par département

```{r cars}
recapdep <- DVFOK %>% group_by(DEP, typo) %>% summarise(nb= n())

ggplot(recapdep, aes(x=DEP, y=nb, fill=typo)) +
  geom_bar(position = 'stack', stat='identity') +
  xlab("") +  ylab("Nombre de mutations") +
  geom_text(aes(label=nb), color="black", size=4, position = position_stack(vjust = 0.5))+
  theme_bw()
```

### Recapitulatif des ventes par TypoINSEE

```{r cars}
recapinsee <- DVFOK %>% group_by(Typo_INSEE, type) %>% summarise(nb= n())

ggplot(recapinsee, aes(x=Typo_INSEE, y=nb, fill=type)) +
  geom_bar(position = 'stack', stat='identity') +
  geom_text(aes(label=nb), color="black", size=4, position = position_stack(vjust = 0.5))+
  xlab("") +  ylab("Nombre de mutations") +
  theme_bw()
```

###Histogramme des prix au m2 par Departement

```{r cars}
RecapPrixDep <- DVFOK %>% group_by(DEP, typo) %>% mutate(moydeptype = mean(prixm2))
RecapPrixDep$ moydeptype <- round(RecapPrixDep$moydeptype)

ggplot(RecapPrixDep, aes(x=prixm2, fill= typo)) +
  geom_histogram(position = "identity", color = "white") +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
  geom_vline(aes(xintercept=moydeptype),color="black", linetype="longdash", size=0.75) +
  xlab("Prix au m2") +  ylab("Effectifs") +
  geom_text(y = 5500, aes(x = moydeptype, label = moydeptype), size = 3, hjust = -.3) +
  facet_grid(typo~DEP)
```

### Histogramme des prix m2 par TypoINSEE
```{r cars}
RecapPrixTypoINSEE <- DVFOK %>% group_by(Typo_INSEE, type) %>% mutate(moyeinseetype = mean(prixm2))
RecapPrixTypoINSEE$ moyeinseetype <- round(RecapPrixTypoINSEE$moyeinseetype)

ggplot(RecapPrixTypoINSEE, aes(x=prixm2, fill= type)) +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
    geom_histogram(position = "identity", color = "white") +
  geom_vline(aes(xintercept=moyeinseetype),color="black", linetype="longdash", size=0.75) +
  xlab("Prix au m2") +  ylab("Effectifs") +
  geom_text(y = 5000, aes(x = moyeinseetype, label = moyeinseetype), size = 3, hjust = -.1) +
  facet_grid(type~Typo_INSEE)
```

### Evolution du nombre de mutations par années et par département

```{r cars}

evoldvfdep <- DVFOK %>% group_by(annee, Dep, type) %>% summarise(nb = n())
View(evoldvfdep)

ggplot(data=evoldvfdep, aes(x=annee, y=nb, fill=type)) +
  geom_bar(stat="identity") +
  xlab("Année") +  ylab("Nb de mutations") +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
  facet_wrap(~Dep, nrow = 1)
```

### Evolutiondu nombre de mutations par années et par type de commune

```{r cars}

evoldvfINSEE <- DVFOK %>% group_by(annee, Typo_INSEE, type) %>% summarise(nb = n())

ggplot(data=evoldvfINSEE, aes(x=annee, y=nb, fill=type)) +
  geom_bar(stat="identity") +
  xlab("Année") +  ylab("Nb de mutations") +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
  facet_wrap(~Typo_INSEE, nrow = 1)

```

### Evolution des prix au m2 par département


```{r cars}
DVFOK$annee <- as.numeric(DVFOK$annee)

EvolPrixDep <- DVFOK %>% group_by(annee, Dep, type) %>% summarise(prix_m2 = mean(prixm2))

ggplot(data=EvolPrixDep, aes(x=annee, y=prix_m2, color=type)) +
  geom_line(stat="identity", size= 1)+
  geom_point(stat="identity", size= 2)+
  scale_y_continuous(breaks=c(1500, 1700,1900, 2100, 2300, 2500, 2700)) +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
  xlab("") +  ylab("Prix moyen au m²") + 
  facet_wrap(~Dep, nrow = 1) 
```


### Evolution des prix au m2 par type de comumune INSEE

```{r cars}
EvolPrixINSEE <- DVFOK %>% group_by(annee, Typo_INSEE, type) %>% summarise(prix_m2 = mean(prixm2))

ggplot(data=EvolPrixINSEE, aes(x=annee, y=prix_m2, color=type)) +
  geom_line(stat="identity", size= 1)+
  geom_point(stat="identity", size= 2)+
  scale_y_continuous(breaks=c(1500, 1700,1900, 2100, 2300, 2500, 2700)) +
  theme_bw() +
  theme(strip.text = element_text(face = "bold")) +
  xlab("") +  ylab("Prix moyen au m²") + 
  facet_wrap(~Typo_INSEE, nrow = 1)
```

### BoxPlot des Prix par type de bien, type de commune et département

```{r cars}
DVFOK$Typo_INSEE <- factor(DVFOK$Typo_INSEE ,levels = c("Espace rural", "Couronne périurbaine", "Pôle urbain"))


ggplot(data = DVFOK, aes(x = Typo_INSEE, y = prix, color = type)) +
  geom_boxplot(notchwidth = 0.5) +
  ylim(0,200000)+
  facet_grid(type~Dep) +
  labs(x= "Type de commune", y= "Prix") +
  theme_bw() +
  theme(strip.text = element_text(face = "bold") ,axis.text.x = element_text(angle = 45, hjust = 1))
```