1.1 Packages nécessaires

# Librairies utilisées
library(sf)
library(cartography)
library(mapsf)
library(corrplot)
library(cowplot)
library(MTA)
library(readxl)
library(ggplot2)
library(FactoMineR) 
library(factoextra)
library(cluster)
library(reshape)
library(reshape2)
library(flows)
# NB : Pour le package flows, la version la plus récente est disponible ici :
# remotes::install_github("rCarto/flows") # ou # install.packages("mapsf")
# Pour obtenir une version plus ancienne (celle utilisée ici) : https://cran.r-project.org/src/contrib/Archive/flows/
# install.packages("packages/flows_1.1.1.tar.gz", repos=NULL, type="source")
library(sp)
library(knitr)
library(condformat) # https://cran.r-project.org/web/packages/condformat/vignettes/introduction.html
library(units)
library(stringr)
# library(dplyr)
library(questionr)
library(spdep) # Pour les matrices de contiguïté
library(rgeoda) # Pour les matrices de contiguïté
library(rmapshaper) # Pour calcul distance aux limites des départements

# Liste pour installer les packages si besoin :
# sf cartography mapsf readxl foreign dplyr flextable knitr stringr units condformat forcats ggplot2 rstatix questionr corrplot gtsummary broom GGally effects forestmodel ggeffects labelled cowplot spdep rgeoda

1.2 Import des couches géographiques consolidées et des données

Les données utilisées sont une compilation de données principalement socio-économiques, localisées. La source principale est l’INSEE. Pour plus de détails, cf. le Markdown “Communes nouvelles : préparation des données”. NB : Ont été exclues des analyses l’Outre-mer et la Corse, espaces non concernés par les communes nouvelles.

Les géométries sont ici importées.

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)  
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)

# Métadonnées
# Liste toutes les variables disponiles
variables_dispo <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ind_target"))
# Liste les variables marquées dans le fichier meta_budgets.xlsx comme nous intéressant
target <- subset(variables_dispo, variable_selec == "X")

Les données socio-économiques qui décrivent les communes en géographies 2011 et 2025 sont ici importées. On commence par extraire les communes ayant participé à la création d’une commune nouvelle, appelées ici communes fusionnantes (datafus2011), les communes nouvelles, avec les géométries au 1er janvier 2025 et caractérisées par les données à la géométrie 2011 agrégées (dataCN_new), ainsi que les communes, à la géométrie 2011, qui n’ont pas participé à la création d’une commune nouvelle (dataNfus2011)

load("data/refdata.Rdata")
datafus2011 <- subset(df2011, COM_NOUV == "OUI")
dataCN_new <- subset(df_new, COM_NOUV == "OUI")
dataNfus2011 <- subset(df2011, COM_NOUV == "NON") 

1.3 Jointure données / géométries

Dans un certain nombre de cas, il sera utile d’avoir, dans un même objet, les données et les géométries. Les données sont ici jointes aux couches géographiques d’intérêt.

geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geom_new <- merge(geom_new, df_new, by = "CODGEO_new")
geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")

1.4 Définition de sous-ensembles pour les tests

Un certain nombre d’espaces peuvent être étudiés à part : certaines communes nouvelles car elles ont été l’objet d’une enquête qualitative (@bideau2017, @bideau2019), d’autres car elles sont sur des territoires très touchés par les fusions communales. Ces sous-ensembles sont cartographiés dans la section 1.7.

testEdC <- subset (geom2011, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")
testEdCSavoie <- subset (geom2011, CODGEO_new =="73150" | CODGEO_new == "73006")
testNormandie <- subset (geom2011, REG == "23" | REG == "25") # Haute et Basse Normandie
test49 <- subset (geom2011, CODE_DEPT == "49" ) # Département du Maine-et-Loire
testOuest <- subset (geom2011, REG == "23" | REG == "25"| REG == "53"| REG == "52") # Normandies, Bretagne, Pays-de-la-Loire


# couleurs <- c("darkblue", "red4")
couleurs <- c("#2b83ba", "red4")
# couleurs <- c("darkblue", "lightblue", "red4")
# couleurs <- c("#2b83ba",  "#fdae61", "#e31a1c")
# couleurs <- c("#2b83ba",  "#fdae61", "red4")

1.5 Analyse des corrélations

RY : pas sûr qu’il faille intégrer cette analyse à la fin. J’ai surtout créé ce bout de code pour réfléchir aux indicateurs à retenir pour les analyses ultérieures (et éviter les variables redondantes dans les CAH par exemple). Les résultats de ces graphiques montrent par exemple et sans surprise l’importante auto-corrélation des variables démographiques

tmp <- na.omit(datafus2011[,50:length(datafus2011)])

# Juste les communes fusionnées
tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, conf.level = .99)


col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")


corrplot(tmp2, method="color",cl.lim = c(-1,1), order = "hclust",
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))


# Communes n'ayant pas fusionné
tmp <- na.omit(dataNfus2011[,50:length(dataNfus2011)])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, method = "pearson")

corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot(tmp2, method="color",cl.lim = c(-1,1), order = "hclust",
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

1.6 Cartographie de localisation

par(mar=c(0,0,1.2,0))

# Localisation des communes nouvelles
layoutLayer(title = "Communes fusionnantes (2011-2025(01))", theme = "red.pal",
            author = "G. Bideau, 2025.",
            sources = "Sources : INSEE, IGN, 2026.", extent = dep)
plot(st_geometry(geomfus2011), col = "red", add = TRUE)
plot(st_geometry(dep), col = NA, add = TRUE)

# Nombre de fusions par départements
tmp <- aggregate(list(NUMBER_FUS = datafus2011[,"CODE_DEPT"]),
                 by = list(DEPT = datafus2011$CODE_DEPT),
                 FUN = length)
dep <- merge (dep, tmp, by.x = "CODE_DEPT", by.y = "DEPT", all.x = TRUE)


dep$NUMBER_FUS_rt <- (dep$NUMBER_FUS / sum(dep$NUMBER_FUS, na.rm = TRUE)) * 100 

layoutLayer(title = "Communes fusionnantes (2011-2025(01))", theme = "red.pal",
            author = "G. Bideau, 2025",
            sources = "Sources : INSEE, IGN, 2025", extent = dep)

plot(st_geometry(dep), col = NA, add = TRUE)

propSymbolsChoroLayer(x = dep, var = "NUMBER_FUS", var2 = "NUMBER_FUS_rt",
                      col = carto.pal(pal1 = "blue.pal", n1 = 3,
                                      pal2 = "red.pal", n2 = 3),
                      inches = 0.2, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.var.pos = "topright", 
                      legend.var2.pos = "topleft",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Pourcentage des communes du\ndépartement ayant fusionné (%)",
                      legend.var.title.txt = "Nombre de fusions\ncommunales",
                      legend.var.style = "e",
                      add = TRUE)

# En fonction du nombre de communes fusionnantes
layoutLayer(title = "Communes fusionnantes (2011-2025)", #theme = "red.pal",
            author = "G. Bideau, 2025",
            sources = "Sources : INSEE, IGN, 2025", extent = dep)
choroLayer(x = dep , var = "NUMBER_FUS_rt",
           method = "quantile", nclass = 6,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           border = NA,
           legend.pos = "topright", legend.values.rnd = 2,
           legend.title.txt = "Pourcentage\nde communes\nfusionnantes",
           add = TRUE)
plot(st_geometry(dep), add = TRUE)
propSymbolsLayer(x = geomCN_new,
                 var = "NbrComFus", symbols ="circle",
                 col =  "red",
                 legend.pos = "topleft",
                 legend.title.txt = "Nombre de communes fusionnantes",
                 legend.style = "e",
                 inches = 0.1,
                 add = TRUE)

# En fonction du nombre de communes fusionnantes sans la trame départementale
layoutLayer(title = "Communes fusionnantes (2011-2025)", #theme = "red.pal",
            author = "G. Bideau, 2025",
            sources = "Sources : INSEE, IGN, 2025", extent = dep)
plot(st_geometry(dep), add = TRUE)
propSymbolsLayer(x = geomCN_new,
                 var = "NbrComFus", symbols ="circle",
                 col =  "red",
                 legend.pos = "topleft",
                 legend.title.txt = "Nombre de communes fusionnantes",
                 legend.style = "e",
                 inches = 0.1,
                 add = TRUE)

# Cartographie des zonages d'études choisis
# Savoie (dont communes enquêtées)
layoutLayer(title = "Communes fusionnantes (2011-2025) : les communes nouvelles (rouge) dont enquêtées (bleu) en Savoie", theme = "red.pal",
            author = "G. Bideau, R. Ysebaert, 2026.",
            sources = "Sources : INSEE, IGN, 2026.", extent = subset(geom_new, CODE_DEPT == "73"))
plot(st_geometry(subset(geom_new, CODE_DEPT == "73")), add = TRUE)
plot(st_geometry(subset(geomCN_new, CODE_DEPT == "73")), col = "red", add = TRUE)
plot(st_geometry(testEdCSavoie), col = "blue", add = TRUE)

# Normandie (dont commune enquêtée)
layoutLayer(title = "Communes fusionnantes (2011-2025) : les communes nouvelles (rouge) dont enquêtée (bleu) en Normandie", theme = "red.pal",
            author = "G. Bideau, R. Ysebaert, 2026.",
            sources = "Sources : INSEE, IGN, 2026.", extent = testNormandie)
plot(st_geometry(testNormandie), add = TRUE)
plot(st_geometry(subset(testNormandie, testNormandie$COM_NOUV == "OUI")), col = "red", add = TRUE)
plot(st_geometry(subset(geomCN_new, CODGEO_new == "61324")), col = "blue", border = NA, add = TRUE)

# Maine-et-Loire
layoutLayer(title = "Communes fusionnantes (2011-2025) : les communes nouvelles (rouge) en Maine-et-Loire", theme = "red.pal",
            author = "G. Bideau, R. Ysebaert, 2026.",
            sources = "Sources : INSEE, IGN, 2026.", extent = test49)
plot(st_geometry(test49), add = TRUE)
plot(st_geometry(subset(test49, COM_NOUV == "OUI")), col = "red", add = TRUE)

1.6bis Cartographie en fonction de la taille démographique

Cf. l’idée donnée ici : https://www.bnsp.insee.fr/ark:/12148/bc6p08tp83z/f1.pdf#page=2

# par(mfrow = c(1,1), mar=c(0,0,1.5,0))

# En fonction de la population de la commune fusionnante

choroLayer(x = geomfus2011 , var = "P09_POP",
           method = "quantile", nclass = 4,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)


layoutLayer(title = " ",# "Communes fusionnantes (2011-2024)\nen fonction de leur nombre d'habitants",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")

# On crée un tableau donnant, pour chaque département, la population moyenne des communes fusionnantes
tableau_dep <- as.data.frame(tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, mean))
colnames(tableau_dep) <- c("moyenne_Cfus")
tableau_dep$median_Cfus <- tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP)
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_dep_tt <- as.data.frame(tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_Cfus),]

kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_Cfus", "moyenne", "median_Cfus", "median")],
      col.names = c("Département", "Nombre de communes fusionnantes", "Population moyenne des communes fusionnantes", "Population communale moyenne", "Population médiane des communes fusionnantes", "Population communale médiane"), digits = 0)

# En fonction de la population de la commune nouvelle
choroLayer(x = geomCN_new , var = "P09_POP",
           method = "quantile", nclass = 4,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)

layoutLayer(title = "",# "Communes nouvelles (2011-2024)\nen fonction de leur nombre d'habitants",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")

# On crée un tableau donnant, pour chaque département, la population moyenne des communes nouvelles
tableau_dep <- as.data.frame(tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, mean))
colnames(tableau_dep) <- c("moyenne_CN")
tableau_dep$median_CN <- tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP) # On garde le nombre de communes fusionnantes pour faciliter les comparaisons entre les deux tableaux
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_dep_tt <- as.data.frame(tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_CN),]

kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_CN", "moyenne", "median_CN", "median")],
      col.names = c("Département", "Nombre de communes nouvelles", "Population moyenne des communes nouvelles", "Population communale moyenne", "Population médiane des communes nouvelles", "Population communale médiane"), digits = 0)

1.6ter Cartographie en fonction de la taille démographique, rapportée aux autres communes

On regarde la position des communes fusionnantes puis des communes nouvelles vis-à-vis des communes françaises

par(mfrow = c(1,2), mar=c(0,0,1.2,0))

# Création des données concernant les quantiles, mais en s'appuyant sur les données de l'ensemble des communes françaises
geomfus2011$P09_POP_quantiles <- as.factor(cut(geomfus2011$P09_POP,
                                     breaks=c(quantile(geom2011$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))
geomCN_new$P09_POP_quantiles <- as.factor(cut(geomCN_new$P09_POP,
                                     breaks=c(quantile(geom_new$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))

# En fonction de la population de la commune fusionnante

typoLayer(x = geomfus2011 , var = "P09_POP_quantiles",
          border = NA,
          legend.values.order = levels(geomfus2011$P09_POP_quantiles),
          col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
          # col = carto.pal(pal1 = "red.pal", n1 = 5),
          legend.pos = "topleft",
          legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)


layoutLayer(title = "Communes fusionnantes (2011-2025)\nen fonction de leur nombre d'habitant",
            author = "G. Bideau, 2025",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2025")


# En fonction de la population de la commune nouvelle
typoLayer(x = geomCN_new , var = "P09_POP_quantiles",
          col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
          border = NA,
          legend.values.order = levels(geomCN_new$P09_POP_quantiles),
          legend.pos = "topleft",
          legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)

layoutLayer(title = "Communes nouvelles (2011-2025)\nen fonction de leur nombre d'habitant",
            author = "G. Bideau, 2025",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2025")

On voit bien que les communes nouvelles sont bien davantage dans la strate supérieure. Mais il ne faut pas croire non plus que seules de petites communes fusionnent, ce n’est pas le cas.

1.7 Analyse multiscalaire de la position des communes nouvelles

Pour les xx ratios sélectionnés, nous souhaitons ici évaluer la position des communes nouvelles dans différents contextes spatiaux.

Pour chaque commune française nous calculons grâce au package MTA (@ysebaert2019a) et pour l’ensemble des indicateurs les déviations suivantes. La valeur 100 correspondant à la moyenne du contexte (définis par un rapport entre le numérateur et le dénominateur) :

• Déviation à la moyenne nationale (devgen).
• Déviation à la moyenne du département d’appartenance de la commune (devdep).
• Déviation à la moyenne de la catégorie d’aire urbaine d’appartenance de la commune (devcatau).
• Déviation à la moyenne des communes situées à moins de 10km des communes

Dans chaque cas, nous différencions les communes fusionnantes des autres (respectivement suffixes fus ou Nfus).

# Transformation en NA des valeurs qui tendent vers l'infini 
geom2011 <- st_set_geometry(geom2011, NULL)
geom2011 <- do.call(data.frame,lapply(geom2011, function(x) replace(x, is.infinite(x),NA)))

# Calculs déviations MTA pour caractériser la position des indicateurs dans différents contextes
ratio <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ratios"))

# Compilation des ratios pour les communes 2011
# Déviation à la moyenne française
for (i in 1:nrow(ratio)){
  geom2011[paste0(ratio[i,"CODE"],"_DEV_GEN")] <-gdev(geom2011, 
                                                      var1 = ratio[i, "Numerator_Code"],
                                                      var2 = ratio[i, "Denominator_Code"],
                                                      type = "rel")
}

# Déviation au département d'appartenance
for (i in 1:nrow(ratio)){
  geom2011[paste0(ratio[i,"CODE"],"_DEV_DEP")] <-tdev(geom2011, 
                                                      var1 = ratio[i, "Numerator_Code"],
                                                      var2 = ratio[i, "Denominator_Code"],
                                                      type = "rel",
                                                      key = "CODE_DEPT")
}

# Déviation à la catégorie d'aire urbaine d'appartenance
for (i in 1:nrow(ratio)){
  geom2011[paste0(ratio[i,"CODE"],"_DEV_CATAU")] <-tdev(geom2011, 
                                                      var1 = ratio[i, "Numerator_Code"],
                                                      var2 = ratio[i, "Denominator_Code"],
                                                      type = "rel",
                                                      key = "CATAEU2010")
}

# Déviation aux unités voisines (moins de 10 km) > Excessivement coûteux en calcul.
# Nécessite davantage d'investigations...

# Calcul matrice de distance > peut être à intégrer à dataprep si on veut 
# se servir de cette matrice de distance à l'avenir car très long à calculer
# st_agr(geom2011) = "constant"
# x.pt <- sf::st_centroid(geom2011, of_largest_polygon = TRUE)
# distmat <- sf::st_distance(x.pt, by_element = FALSE)
# colnames(distmat) <- x.pt[["CODGEO"]]
# rownames(distmat) <- x.pt[["CODGEO"]]
# attr(distmat, "units") <- NULL
# class(distmat) <- setdiff(class(distmat),"units")
# distmat[distmat <= 10000] <- 1
# distmat[distmat > 10000] <- 0
# 
# try(save(distmat, file = "data/distMat.Rdata"))

# for (i in 1:nrow(ratio)){
#   geom2011[paste0(ratio[i,"CODE"],"_DEV_10K")] <-sdev(geom2011, 
#                                                       var1 = ratio[i, "Numerator_Code"],
#                                                       var2 = ratio[i, "Denominator_Code"],
#                                                       type = "rel",
#                                                       dist = 1,
#                                                       mat = distmat)
# }



# Synthèse déviation générale
devgen <- geom2011[,endsWith(colnames(geom2011), "DEV_GEN")]
devgen <- cbind(devgen, geom2011$COM_NOUV)
colnames(devgen)[length(devgen)] <- "COM_NOUV"
# Médiane des communes non fusionnantes
devgenNfus <- devgen[devgen$COM_NOUV == "NON",]
devgenfus <- devgen[devgen$COM_NOUV == "OUI",]
devgenNfus <- as.data.frame(apply(devgenNfus[,1:length(devgenNfus)-1], 2, median, na.rm = TRUE))
devgenfus <- as.data.frame(apply(devgenfus[,1:length(devgenfus)-1], 2, median, na.rm = TRUE))

# Synthèse déviation départements
devdep <- geom2011[,endsWith(colnames(geom2011), "DEV_DEP")]
devdep <- cbind(devdep, geom2011$COM_NOUV)
colnames(devdep)[length(devdep)] <- "COM_NOUV"
# Médiane des communes non fusionnantes
devdepNfus <- devdep[devdep$COM_NOUV == "NON",]
devdepfus <- devdep[devdep$COM_NOUV == "OUI",]
devdepNfus <- as.data.frame(apply(devdepNfus[,1:length(devdepNfus)-1], 2, median, na.rm = TRUE))
devdepfus <- as.data.frame(apply(devdepfus[,1:length(devdepfus)-1], 2, median, na.rm = TRUE))


# Synthèse déviation catégories aires urbaines
devcatau <- geom2011[,endsWith(colnames(geom2011), "DEV_CATAU")]
devcatau <- cbind(devcatau, geom2011$COM_NOUV)
colnames(devcatau)[length(devcatau)] <- "COM_NOUV"
# Médiane des communes non fusionnantes
devcatauNfus <- devcatau[devcatau$COM_NOUV == "NON",]
devcataufus <- devcatau[devcatau$COM_NOUV == "OUI",]
devcatauNfus <- as.data.frame(apply(devcatauNfus[,1:length(devcatauNfus)-1], 2, median, na.rm = TRUE))
devcataufus <- as.data.frame(apply(devcataufus[,1:length(devcataufus)-1], 2, median, na.rm = TRUE))

Ces éléments nous permettent de comparer la déviation à la moyenne des communes fusionnantes et des communes inchangées. On peut, d’ailleurs, calculer la différence entre ces éléments.

Le tableau suivant propose une synthèse des données ainsi construites. Les couleurs sont construites en fonction des minima et maxima des valeurs de chaque partie du tableau (déviations et différences).

           Différence (%) entre communes nouvelles et communes inchangées concernant la déviation à la moyenne française 
                                                                                                                   29.64 
      Différence (%) entre communes nouvelles et communes inchangées concernant la déviation à la moyenne départementale 
                                                                                                                    5.32 

Différence (%) entre communes nouvelles et communes inchangées concernant la déviation aux mêmes catégories d’aires urbaines 6.95 Différence (%) entre communes nouvelles et communes inchangées concernant la déviation à la moyenne française -13.57 Différence (%) entre communes nouvelles et communes inchangées concernant la déviation à la moyenne départementale -10.07 Différence (%) entre communes nouvelles et communes inchangées concernant la déviation aux mêmes catégories d’aires urbaines -9.84

Le premier élément remarquable dans ce tableau est la proximité entre les communes fusionnantes et les communes inchangées concernant la déviation à la moyenne, qu’elle soit française et encore plus départementale ou liée à la catégorie d’aire urbaine. Cependant, certaines variables sortent du lot et pointent une certaine spécificité des communes nouvelles, que ce soit en comparant à la moyenne française, départementale ou de ZAU : la part des agriculteurs (ou de l’agriculture, sur-représentation), la part des ouvriers (ou de l’industrie, sur-représentation) et la part des actifs travaillant dans leur commune (sous-représentation). Sont également sous-représentés des éléments comme la part de la construction, de l’administration, des professions intermédaires et des cadres. Cela dessine le profil de communes fusionnantes globalement proche du profil moyen des communes françaises mais ayant comme spécificités des populations plus agricoles, plus ouvrières et, souvent, travaillant davantage hors de leur commune d’habitation.

1.8 Graphique des fusions en fonction du temps

Plusieurs vagues de fusions peuvent être notées depuis la création du statut de communes nouvelles.

Le graphique ci-dessous représente le nombre de communes participant à une fusion.

Dans le premier graphique, les fusions du 1er janvier de l’année n sont rattachées à l’année n-1 (donc, pour l’année 2019 par exemple, du 2 janvier 2019 au 1er janvier 2020 inclus). Cela s’explique par le faite que de très nombreuses fusions sont décidées en amont pour prendre effet au 1er janvier de l’année suivante. Dans ce cas, elle sont considérées comme existantes pour toute l’année suivante.

Le second graphique, pour être complet, présente les données en fonction de l’année d’effectivité (1er janvier 2019 rattaché à l’année 2019).

Cette distinction est particulièrement importante pour l’année 2019 justement, puisque le nombre important de fusions associé à cette année est lié à des fusions au 1er janvier 2019, donc avant le vote de la loi du 1er août.

datafus2011$FusDate <- as.Date(datafus2011$FusDate, tryFormats = "%Y-%m-%d")

datafus2011$FusDateAnnee <- cut(datafus2011$FusDate,
                              breaks = as.Date(c("2011-1-1", "2012-1-2", "2013-1-2", "2014-1-2", "2015-1-2", "2016-1-2", "2017-1-2","2018-1-2", "2019-1-2", "2020-1-2", "2021-1-2", "2022-1-2", "2023-1-2", "2024-1-2")), 
                        labels = c("2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022", "2023")) # NB : On intègre le 1er janvier d'une année n à l'année précédente n-1

table(datafus2011$FusDateAnnee)

## 
## 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 
##    2   29    2   37 1082  653   91  600    6    5   19   18   22

time_serie <- data.frame(table(datafus2011$FusDateAnnee))

par(bg = "darkgrey", fg = "white", lwd = 2)

barplot(time_serie$Freq,  xlab = " ", ylab = "Nombre de communes fusionnantes", main = "Nombre de communes fusionnantes en fonction de l'année de décision\n(1er janvier de l'année n rattaché à l'année n-1)", ylim = c(0, 2000), names.arg = time_serie$Var1)
grid(nx=0, ny = 10, col="gray", lwd=2)

# Indication d'éléments dans le graphique.
arrows(0.1, 1500, 0.1, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( -0.4, 1700 , labels = "Loi du 16 décembre 2010\ncréant le statut de\n« commune nouvelle »", col = "black", cex = 0.9, pos = 4, adj = c(1, 1))

arrows(3.7, 600, 3.7, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( 3.7, 1000 , labels = "Premières incitations\nfiscales (art. 133,\nloi de finances pour 2014)", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

arrows(5.3, 1300, 5.3, 1150, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text(5.3, 1800 , labels = "Dotation globale de\nfonctionnement\nmajorée de 5%\n(loi du 16 mars 2015)", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

arrows(10.3, 1200, 10.3, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( 10, 1500 , labels = "Loi du 1er août 2019\nmodifiant l'incitation\nfiscale (6€ par habitant)\net créant le statut de\ncommune-communauté", col = "black", cex = 0.9, pos = 4, adj = c(1, 1))

datafus2011$FusDateAnnee_effectivite <- cut(datafus2011$FusDate,
                        breaks = as.Date(c("2011-1-1", "2012-1-1", "2013-1-1", "2014-1-1", "2015-1-1", "2016-1-1", "2017-1-1","2018-1-1", "2019-1-1", "2020-1-1", "2021-1-1", "2022-1-1", "2023-1-1", "2024-1-1", "2025-1-1")),
                        labels = c("2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022", "2023", "2025")) # Si on utilise l'année (1er janvier laissé dans l'année concernée)

table(datafus2011$FusDateAnnee_effectivite)

## 
## 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2025 
##    0    2   31    0   46 1165  564   90  604    0    5   19   18   22

time_serie <- data.frame(table(datafus2011$FusDateAnnee_effectivite))

par(bg = "darkgrey",fg = "white", lwd = 2)
par(bg = "white",fg = "gray", lwd = 2)

barplot(time_serie$Freq,  xlab = " ", ylab = "Nombre de communes fusionnantes", main = "Nombre de communes fusionnantes en fonction de l'année d'effectivité", ylim = c(0, 2000), names.arg = time_serie$Var1)
grid(nx=0, ny = 10, col="gray", lwd=2)

# Indication d'éléments dans le graphique.
arrows(0.1, 1500, 0.1, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( -0.4, 1700 , labels = "Loi du 16 décembre 2010\ncréant le statut de\n« commune nouvelle »", col = "black", cex = 0.9, pos = 4, adj = c(1, 1))

arrows(3.5, 600, 3.5, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( 3.5, 1000 , labels = "Premières incitations\nfiscales (art. 133,\nloi de finances\npour 2014)", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

arrows(5.3, 1300, 5.3, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text(5.3, 1800 , labels = "Dotation globale de\nfonctionnement\nmajorée de 5%\n(loi du 16 mars 2015)", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

arrows(10.3, 1200, 10.3, 700, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( 10, 1700 , labels = "Loi du 1er août 2019\nmodifiant l'incitation\nfiscale (6€ par habitant)\net créant le statut de\ncommune-communauté", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

arrows(15.6, 700, 15.6, 100, length = 0.10, angle = 30,
       code = 2, col = "black", lty = par("lty"),
       lwd = par("lwd"))
text( 15.6, 1200 , labels = "Dotation d'amorçage\nrevalorisée et dotation\nde garantie sans limite\nde temps (loi de\nfinances pour 2024)", col = "black", cex = 0.9, pos = 1, adj = c(1, 1))

Cette chronologie peut être complétée par des éléments concernant les évolutions législatives :

• Loi du 16 décembre 2010 créant le statut de commune nouvelle

• Premières incitations fiscales (art. 133 de la loi de Finances pour 2014)

• Dotation globale de fonctionnement majorée de 5 % (loi du 16 mars 2015)

• Loi du 1er août 2019 créant le statut de commune-communauté

• Loi de Finances pour 2020 créant une incitation fiscale spécifique (6€ par habitant)

Plusieurs vagues de fusions peuvent être notées depuis la création du statut de communes nouvelles.

La première va du 1er janvier 2011 au 1er janvier 2015 et concerne 70 communes fusionnantes.

La seconde va du 2 janvier 2015 au 1er janvier 2016 et concerne 1085 communes.

La troisième va du 2 janvier 2016 au 1er janvier 2017 et concerne 658 communes.

La quatrième va du 2 janvier 2017 au 1er janvier 2017 et concerne 91 communes.

La cinquième va du 2 janvier 2018 au 1er janvier 2020 et concerne 608 communes.

La sixième va du 2 janvier 2020 au 1er janvier 2021 et concerne 5 communes. [NB : À compléter avec les communes créées depuis.]

datafus2011$FusDate <- as.Date(datafus2011$FusDate, tryFormats = "%Y-%m-%d")
dev.off()

## null device 
##           1

time_serie_phas <- data.frame(table(datafus2011$FusPhas))
barplot(time_serie_phas$Freq, names.arg = time_serie_phas$Var1)

On peut également représenter les fusions en fonction de leur date exacte

time_serie <- data.frame(table(datafus2011$FusDate)) # Si on utilise la date exacte
time_serie$Var1 <- as.Date(time_serie$Var1, tryFormats = "%Y-%m-%d")


par(bg = "darkgrey", fg = "white", lwd = 2)
plot(time_serie$Var1, time_serie$Freq, pch = 20, cex = 1.5, col = "black", xlab = "Date exacte (jour)",
ylab = "Nombre de communes fusionnantes", main = "Nombre de communes fusionnantes en fonction de la date exacte d'effectivité", ylim = c(0, 2000))
grid(lwd=2)
lines(time_serie$Var1, time_serie$Freq, col = 2, lwd = 2)

1.9 Cartographie par phases

col <- c("#ffffbf", "#fee090", "#fdae61", "#f46d43", "#d73027", "#a50026")  # Échelle rouge vers le plus foncé
col <- c("#ffffbf", "#e0f3f8","#abd9e9","#74add1","#4575b4","#313695") # Échelle bleue vers le plus foncé
col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090", "#ffffbf") # Échelle rouge vers le plus clair
plot(st_geometry(dep), border = "#1A1A19", col = "white", lwd = 1)
typoLayer(x = geomfus2011, var = "FusPhas",  
          col=col,
          border = NA, 
          legend.title.cex = 0.7,
          legend.values.cex = 0.6,
          # 
          legend.pos = "left", 
          legend.values.order = c("Phase 1", "Phase 2", "Phase 3", "Phase 4", "Phase 5", "Phase 6"),
          add = T)

layoutLayer(title = "Communes fusionnantes (2012-2025) par phases",
            author = "G. Bideau",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2025")

1.10 La population des communes fusionnantes

La taille est souvent considérée comme le premier facteur menant à une rationalisation de l’échelon communal par fusion. Les municipalités sont effectivement une maille héritée des paroisses, territoires du quotidien religieux, ayant perdu de leur sens dans le cadre des grandes dynamiques démographiques du XIXe et XXe siècle et n’étant plus nécessairement adaptées aux compétences communales. On pourrait alors s’attendre à une proportion de communes fusionnantes bien plus importante parmi les petites communes. On compare ici la proportion des communes en fonction de leur taille, en mettant en parallèle l’ensemble des communes françaises et les communes fusionnantes.

# Moyenne et médiane de l'ensemble des communes-françaises
mean(df2011$P09_POP)

## [1] 1716.749

median(df2011$P09_POP)

## [1] 426

# Moyenne et médiane des communes fusionnantes
mean(datafus2011$P09_POP)

## [1] 1012.266

median(datafus2011$P09_POP)

## [1] 404

PopCom <- df2011$P09_POP
PopCom <- cut(PopCom, breaks = c(-1,200,500,1000, 5000, 10000, 2230000))
levels(PopCom)<-c("moins de 200 habitants", "200-500 habitants", "500-1000 habitants", "1000-5000 habitants", "5000-10000 habitants", "plus de 10000 habitants")
levels(PopCom)<-c("moins de 200", "200-500", "500-1000", "1000-5000", "5000-10000", "plus de 10000")
table (PopCom)

## PopCom
##  moins de 200       200-500      500-1000     1000-5000    5000-10000 
##          9377         10469          6948          7480          1046 
## plus de 10000 
##           887

# Pourcentage de communes de moins de 200 habitants avant la fusion
table (PopCom)[1] / length (PopCom)

## moins de 200 
##    0.2589759

# Pourcentage de communes de moins de 200 habitants après la fusion
PopCom_new <- df_new$P09_POP
PopCom_new <- cut(PopCom_new, breaks = c(-1,200,500,1000, 5000, 10000, 2230000))
levels(PopCom_new)<-c("moins de 200 habitants", "200-500 habitants", "500-1000 habitants", "1000-5000 habitants", "5000-10000 habitants", "plus de 10000 habitants")
levels(PopCom_new)<-c("moins de 200", "200-500", "500-1000", "1000-5000", "5000-10000", "plus de 10000")
table (PopCom_new)

## PopCom_new
##  moins de 200       200-500      500-1000     1000-5000    5000-10000 
##          8720          9669          6606          7403          1082 
## plus de 10000 
##           898

table (PopCom_new)[1] / length (PopCom_new)

## moins de 200 
##    0.2536432

tabcont<-table(df2011$FUSION, PopCom)

# Test chi²
tab.chi2 <- chisq.test(tabcont)
tab.chi2

## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 57.647, df = 5, p-value = 3.719e-11

tab <- round(100*prop.table(tabcont,margin=1),2) # Pourcentages, le total se fait par lignes

# Pour avoir un graphique
prbarplot <- as.data.frame(tab)
ggplot(data = prbarplot, aes(x = PopCom, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  theme(axis.text.x = element_text(angle = 45, hjust=1), legend.position = c(0.8, 0.85)) +
  scale_x_discrete("Nombre d'habitants") + 
  scale_y_continuous("Part du groupe de communes (%)") +
  # scale_fill_manual("Communes\nfusionnantes", values = couleurs) +
  scale_fill_manual(labels = c("Communes inchangées", "Communes fusionnantes"), values = couleurs) +
  labs(fill = NULL)

# Proportion sur l'ensemble des communes françaises
Ensemble <- round(100* table(PopCom) / length(PopCom),2)
tab <- rbind(tab, Ensemble)

# Pour transposer (lignes deviennent colonnes et vice versa)
tab <- t(tab)
Total <- apply(tab, 2, sum)
tab <- rbind (tab, Total)

kable(tab, row.names = T, digits = 2, caption = "Taille des communes (nombre d'habitants) et fusions")

# col.names = c("Moins de 200", "200-500", "500-1 000", "1 000-5 000", "5000-10 000", "Plus de 10 000", "Total"

rm(tab, table, tabcont, Ensemble)

1.11 Quelques statistiques basiques

1.12 Communes nouvelles et périmètres de fusions

NB : Pour l’appartenance au 1er janvier 2025, pas d’information autres que les appartenances administratives (région, arrondissement, etc.). J’espère que cela viendra rapidement (note au 14/03/2025).

# Création d'un tableau pour l'appartenance de toutes les communes ayant créé une commune nouvelle
variables_appart <- c("CODGEO", "LIBGEO", "CODGEO_new", "LIBGEO_new", "CODE_DEPT", "CATAEU2010", "REG", "ARR", "CV", "UU2010", "AU2010", "ZE2010", "EPCI", "ChefLieu", "superficie")

AppartComFus <- datafus2011[, variables_appart]
AppartCom <- df2011[, c(variables_appart)]

# NB : Différence due à L'Oudon

# Nombre de communes ayant fusionné dans chaque EPCI, rapporté au nombre de communes total
# En rapide :
EPCIComNvll <- merge(data.frame(table(AppartComFus$EPCI)),
                     data.frame(table(AppartCom$EPCI)), by = "Var1", all.x = TRUE)
# En développé :
# EPCINbrFus <- data.frame(table(AppartComFus$EPCI))
# EPCINbrCom <- data.frame(table(AppartCom$EPCI))
# EPCIComNvll <- merge(EPCINbrFus, EPCINbrCom, by = "Var1", all.y = TRUE)

# Opérations pour créer la nouvelle donnée :
EPCIComNvll <- rename.variable(EPCIComNvll, "Freq.x", "NombreCommunesFusionnant")
EPCIComNvll <- rename.variable(EPCIComNvll, "Freq.y", "NombreCommunesTotal")
EPCIComNvll <- rename.variable(EPCIComNvll, "Var1", "EPCI")
EPCIComNvll$NombreCommunesFusionnant [is.na(EPCIComNvll$NombreCommunesFusionnant)] <- 0
EPCIComNvll$RatioComFus = EPCIComNvll$NombreCommunesFusionnant / EPCIComNvll$NombreCommunesTotal * 100
# On retire les communes n'appartenant pas à un EPCI
EPCIComNvll <- subset (EPCIComNvll, EPCIComNvll$EPCI != "ZZZZZZZZZ")
plot(sort(EPCIComNvll$RatioComFus, decreasing = TRUE),
     xlab= "Nombre d'EPCI dans lequel au moins une commune a fusionné",
     ylab= "Pourcentage de communes ayant fusionné au sein de l'EPCI",
     main = "Les communes fusionnantes au sein des EPCI")

# mean(EPCIComNvll$RatioComFus)
# median(EPCIComNvll$RatioComFus)
# quantile (EPCIComNvll$RatioComFus)
# sd(EPCIComNvll$RatioComFus)
summary (EPCIComNvll$RatioComFus)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.222  12.500  22.222  31.341  40.000 100.000

# Choix des périmètres d'appartenance pris en compte
Appart_etud <- c("REG", "CODE_DEPT", "ARR", "CV", "ZE2010", "AU2010", "EPCI")

# Création d'un tableau pour l'appartenance de la nouvelle commune à partir du Chef-Lieu
AppartComNvlles <- subset (AppartComFus, ChefLieu == "O")

# Création d'un tableau avec les caractéristiques : de l'ancienne commune (les ".x") et du chef-lieu de la commune nouvelle (les ".y").
# À noter qu'on n'a pris que les communes n'étant pas chefs-lieu, pour éviter les double-comptes
AppartComFus2 <- merge(
  subset(AppartComFus, ChefLieu == "N"),
  AppartComNvlles, by = "CODGEO_new", all.x = TRUE)

appart <- "AU2010"

tableau <- data.frame()

for (appart in Appart_etud) {
# Création d'un champ désignant, pour chaque zonage, si la commune fusionnante a la même appartenance que la commune chef-lieu.
  AppartComFus2[, paste0(appart, "_simil")] <- ifelse(AppartComFus2[, paste0(appart, ".x")] == AppartComFus2[, paste0(appart, ".y")], TRUE, FALSE)
  
  # Modifications particulières pour ZAU car l'aire urbaine notée "000" regroupe toutes les communes isolées hors influence des pôles. Une appartenance identique à cette dernière n'a donc pas beaucoup de sens.
  
   if (appart != "AU2010") {AppartComFus2 <- AppartComFus2} else {
     AppartComFus2$AU2010_simil[AppartComFus2$AU2010.x == "000"] <- "Non applicable"
     AppartComFus2$AU2010_simil[AppartComFus2$AU2010.y == "000"] <- "Non applicable"

   }
  
    # Modifications particulières pour EPCI car, de même, l'absence d'EPCI est notée "ZZZZZZZZZ" (neux fois la lettre Z)
   if (appart != "EPCI") {AppartComFus2 <- AppartComFus2} else {
     AppartComFus2$EPCI_simil[AppartComFus2$EPCI.x == "ZZZZZZZZZ"] <- "Non applicable"
     AppartComFus2$EPCI_simil[AppartComFus2$EPCI.y == "ZZZZZZZZZ"] <- "Non applicable"
     }
  
nbSimil <- sum(AppartComFus2[, paste0(appart, "_simil")] == TRUE, na.rm = TRUE)
nbDiff <- sum(AppartComFus2[, paste0(appart, "_simil")] == FALSE, na.rm = TRUE)
nbCom <- nbSimil + nbDiff


num_variable <- which(target$recoding == appart)
nom_variable <- target[num_variable, "DESCRIPTION"]

# On rajoute une mention de la superficie moyenne de l'échelon
superficie_appart <- tapply(df2011$superficie, INDEX = df2011[, appart], sum)
superf_moy <- round(mean(superficie_appart),0)

ligne <- c(nom_variable, nbSimil, 100*nbSimil/nbCom, nbDiff, 100*nbDiff/nbCom, nbCom, 100*nbCom/nbCom, superf_moy)

tableau <- rbind(tableau, ligne)

}


# On passe les variables en numérique
tableau[2:ncol(tableau)] <- apply(tableau[2:ncol(tableau)], 2, as.numeric)

colnames(tableau) <- c("", "Nombre", "%", "Nombre", "%", "Nombre", "%", "(km²)")

# Si on veut ordonner le tableau en fonction de la superficie moyenne décroissante
# tableau <- tableau[order(-tableau$`(km²)`),]

kable(tableau, row.names = F, digits = 2, align = "c", caption = "L'appartenance des communes fusionnantes à différents zonages") %>%
  kableExtra::add_header_above(c("Appartenance identique" = 1, "Oui" = 2, "Non" = 2, "Total" = 2, "Superficie moyenne de l'échelon" = 1))

# Réalisation d'un diagramme

temp1 <- tableau[, c(1, 2)]
temp1$Appart <- "Appartenance identique"

temp2 <- tableau[, c(1, 4)]
temp2$Appart <- "Appartenance différente"

tableau_plot <- rbind(temp1, temp2)
colnames(tableau_plot) <- c("Échelon", "Nombre", "Appartenance")


# Réalisation d'un diagramme
ggplot(tableau_plot, aes(x = Échelon,
                              y = Nombre,
                              fill = Appartenance)) + geom_bar(stat="identity") +
  # Changement de l'ordre dans la légende
  scale_x_discrete(limits=c("Région", "Département", "Arrondissement", "Canton de ville", "Zone d'emploi 2010", "Aire urbaine 2010", "EPCI")) +
  labs(title="Nombre de communes fusionnantes\nayant une appartenance identique à celle\nde la commune chef-lieu de la commune nouvelle")

rm(appart, nbSimil, nbDiff, nbCom, tableau, num_variable, nom_variable, superficie_appart, superf_moy, tableau_plot, temp1, temp2)

1.13 Communes nouvelles et superficie

EtudSuperf <- df2011[, c("CODGEO", "FUSION", "superficie")]

seuils <- quantile(EtudSuperf$superficie, probs=seq(0, 1, 0.1)) # Seuils à partir des déciles de l'ensemble des communes fusionnantes

EtudSuperf$deciles <- EtudSuperf$superficie
EtudSuperf$deciles <- cut(EtudSuperf$deciles, breaks = seuils)
# levels(EtudSuperf$deciles)<-c("1er décile", "2e décile", "3e décile", "4e décile", "5e décile", "6e décile", "7e décile", "8e décile", "9e décile", "10e décile")
table (EtudSuperf$deciles)

## 
## (0.0165,4.08]   (4.08,5.65]   (5.65,7.19]   (7.19,8.84]   (8.84,10.8] 
##          3620          3621          3621          3620          3621 
##   (10.8,13.1]   (13.1,16.3]   (16.3,21.1]   (21.1,29.9]    (29.9,758] 
##          3621          3620          3621          3621          3621

tabcont<-table(EtudSuperf$FUSION, EtudSuperf$deciles)
# Test chi²
tab.chi2 <- chisq.test(tabcont)
tab.chi2

## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 17.647, df = 9, p-value = 0.0395

# Pour avoir un graphique avec les barres juxtaposées en pourcentage des lignes
tab <- round(100*prop.table(tabcont,margin=1),2) # Pourcentages, le total se fait par lignes
prbarplot <- as.data.frame(tab)
ggplot(data = prbarplot, aes(x = Var2, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  scale_x_discrete("Superficie communale (km²)") + 
  scale_y_continuous("Part du groupe de communes (%)") +
  scale_fill_manual("Communes\nfusionnantes", values = couleurs) 

  # ggtitle("La répartition des communes fusionnantes\nen fonction des superficies\n(par déciles sur l'ensemble des communes françaises)")

# Pour avoir un graphique empilé avec le pourcentage total des communes françaises
tab <- round(100*prop.table(tabcont),2) # Pourcentages, le total se fait par lignes
prbarplot <- as.data.frame(tab)
ggplot(data = prbarplot, aes(x = Var2, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  scale_x_discrete("Superficie communale (km²)") + 
  scale_y_continuous("Part du groupe de communes (%)") +
  scale_fill_manual("Communes\nnouvelles", values = couleurs) +
  ggtitle("La répartition des communes fusionnantes\nen fonction des superficies\n(par déciles sur l'ensemble des communes françaises)")  +
  geom_hline(yintercept=100 * nrow(datafus2011)/nrow(df2011) /10)

rm(tab, tabcont, seuils, EtudSuperf)

Si on observe la superficie communale, la répartition des communes fusionnantes ou non laisse entrevoir un phénomène potentiellement contre-intuitif : la participation à une commune nouvelle ne paraît pas du tout être guidée par la taille des communes. Alors que les communes avec une très petite superficie (premier décile sur l’ensemble des communes françaises) ont effectivement un peu plus fusionné en moyenne (mais très légèrement), celles avec des superficies au-dessus de la moyenne (pour les classes entre 13,1 et 29,9 km²) ont elles aussi davantage fusionné. Dans tous les cas, les variations sont généralement faibles et la question de la superficie ne paraît pas être un élément explicatif, ce qui est confirmé par un test de Chi² ne permettant pas de rejeter l’hypothèse d’indépendance des deux variables (X-squared = 13.84, df = 9, p-value = 0.1281).

EtudSuperf <- df_new[, c("CODGEO", "FUSION", "superficie")]

seuils <- quantile(EtudSuperf$superficie, probs=seq(0, 1, 0.1)) # Seuils à partir des déciles de l'ensemble des communes fusionnantes

EtudSuperf$deciles <- EtudSuperf$superficie
EtudSuperf$deciles <- cut(EtudSuperf$deciles, breaks = seuils)
# levels(EtudSuperf$deciles)<-c("1er décile", "2e décile", "3e décile", "4e décile", "5e décile", "6e décile", "7e décile", "8e décile", "9e décile", "10e décile")
table (EtudSuperf$deciles)

## 
## (0.0165,4.13]   (4.13,5.73]   (5.73,7.28]   (7.28,8.99]     (8.99,11] 
##          3437          3438          3438          3438          3438 
##     (11,13.4]   (13.4,16.7]   (16.7,21.8]   (21.8,31.7]    (31.7,758] 
##          3437          3438          3438          3438          3438

tabcont<-table(EtudSuperf$FUSION, EtudSuperf$deciles)


# Pour avoir un graphique avec les barres juxtaposées en pourcentage des lignes
tab <- round(100*prop.table(tabcont,margin=1),2) # Pourcentages, le total se fait par lignes
prbarplot <- as.data.frame(tab)
ggplot(data = prbarplot, aes(x = Var2, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  scale_x_discrete("Superficie communale (km²)") + 
  scale_y_continuous("Pourcentage de communes") +
  scale_fill_manual("Communes\nfusionnantes", values = c("blue", "red")) +
  ggtitle("La répartition des communes fusionnantes\nen fonction des superficies\n(par déciles sur l'ensemble des communes françaises)")

# Pour avoir un graphique empilé avec le pourcentage total des communes françaises
tab <- round(100*prop.table(tabcont),2) # Pourcentages, le total se fait par lignes
prbarplot <- as.data.frame(tab)
ggplot(data = prbarplot, aes(x = Var2, y = Freq, fill = Var1)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  scale_x_discrete("Superficie communale (km²)") + 
  scale_y_continuous("Pourcentage de communes") +
  scale_fill_manual("Communes\nfusionnantes", values = c("blue", "red")) +
  ggtitle("La répartition des communes fusionnantes\nen fonction des superficies\n(par déciles sur l'ensemble des communes françaises)")  +
  geom_hline(yintercept=100 * nrow(dataCN_new)/nrow(df_new) /10) # Ligne représentant la moyenne des communes nouvelles par décile

rm(tab, tabcont, seuils, EtudSuperf)

En revanche, le passage en commune nouvelle a eu un effet majeur puisque les communes nouvelles sont bien plus présentes dans le décile avec la superficie les plus élevée que dans les autres.

2.0 Préparation des données

load("data/refdata.Rdata")
datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes
dataNfus2011 <- subset(df2011, COM_NOUV == "NON") # Les communes, à la géométrie 2011, qui n'ont pas participé à la création d'une commune nouvelle
datafus2011 <- datafus2011[, c("CODGEO", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "CATAEU2010", "FusPhas")]
dataNfus2011 <- dataNfus2011 [, c("CODGEO", "CODGEO_new", "LIBGEO", "CODGEO", "CATAEU2010")]
c("CATAEU2010", "ZAU_POL", "ZAU_RUR", "ZAU_MAR_SP", "ZAU_MAR", "ZAU_PERI", "ZAU_AU")

## [1] "CATAEU2010" "ZAU_POL"    "ZAU_RUR"    "ZAU_MAR_SP" "ZAU_MAR"   
## [6] "ZAU_PERI"   "ZAU_AU"

couleurs <- c("#2b83ba", "red4")

2.1 Les communes nouvelles en fonction des ZAU

À partir des données INSEE précisant le zonage en aire urbaine de chaque commune (CATAEU2010), nous pouvons observer le profil des communes fusionnantes en le comparant à celui des communes inchangées.

ZAU_non <- data.frame(dataNfus2011$CATAEU2010)
ZAU_oui <- data.frame(datafus2011$CATAEU2010)

ZAU_non$Fusion <- "Communes inchangées"
ZAU_oui$Fusion <- "Communes fusionnantes"

colnames(ZAU_non)[1] <- "CATAEU2010"
colnames(ZAU_oui)[1] <- "CATAEU2010"

ZAU <- rbind(ZAU_non, ZAU_oui)

ZAU_nb <- table(ZAU$Fusion, ZAU$CATAEU2010)

prop <- prop.table(ZAU_nb, 1) * 100

couleurs_ici <- rev(couleurs)
barplot(prop,
        xlab = "Catégories d'aire urbaine (ZAU INSEE)",
        ylab = "Part du total des communes (%)",
        # names.arg = c("Commune appartenant à un grand pôle", "Commune appartenant à la couronne d'un grand pôle", " Commune multipolarisée des grandes aires urbaines", "Commune appartenant à un moyen pôle", " Commune appartenant à la couronne d'un moyen pôle", "Commune appartenant à un petit pôle", "Commune appartenant à la couronne d'un petit pôle", "Autre commune multipolarisée", "Commune isolée hors influence des pôles"),
       # names.arg = c("Grand pôle", "Couronne d'un grand pôle", " Multipolarisée des\ngrandes aires urbaines", "Moyen pôle", "Couronne d'un moyen pôle", "Petit pôle", "Couronne d'un petit pôle", "Autre commune\nmultipolarisée", "Commune isolée hors\ninfluence des pôles"),

        main = " ",# "Graphique 2 :\nRépartition des communes fusionnantes ou non\nen fonction des catégories d'aire urbaine",
        las = 2,
        border = NA,
        col=couleurs_ici, # "#ff87a9","#f7d358"
        beside = TRUE)

legend(x="topright", legend = rownames(prop) , cex=0.8,
       fill=couleurs_ici,bty="n")     

prbarplot <- melt(prop)
prbarplot$Var1 <- factor(prbarplot$Var1, levels = rev(levels(prbarplot$Var1)))

prbarplot$Var2 <- as.character(prbarplot$Var2)
ggplot(data = prbarplot, aes(x = Var2, y = value, fill = Var1)) +
  geom_bar(stat = "identity", position=position_dodge()) +
  theme(axis.text.x = element_text(angle = 45, hjust=1), legend.position = c(0.8, 0.85)) +
  scale_x_discrete("Catégorie de Zonage en Aire Urbaine (ZAU)",
                   labels = c("Grand pôle", "Couronne\nd'un grand pôle", " Commune\nmultipolarisée\n par des grandes\naires urbaines", "Moyen pôle", "Couronne d'un\nmoyen pôle", "Petit pôle", "Couronne\nd'un petit pôle", "Autre commune\nmultipolarisée", "Commune isolée hors\ninfluence des pôles")) + 
  scale_y_continuous("Part du groupe de communes (%)") +
  scale_fill_manual(values = couleurs) + labs(fill = NULL)

# Chi2
tab.chi2 <- chisq.test(ZAU_nb)
tab.chi2 

## 
##  Pearson's Chi-squared test
## 
## data:  ZAU_nb
## X-squared = 181.24, df = 8, p-value < 2.2e-16

tab.chi2$observed

##                        
##                           111   112   120   211   212   221   222   300   400
##   Communes fusionnantes    81   827   310    24    74   100    32   590   633
##   Communes inchangées    3119 11353  3658   415   727   754   526  6385  6600

tab.chi2$expected

##                        
##                               111       112       120       211       212
##   Communes fusionnantes  236.0583   898.497  292.7123  32.38425  59.08835
##   Communes inchangées   2963.9417 11281.503 3675.2877 406.61575 741.91165
##                        
##                               221       222       300       400
##   Communes fusionnantes  62.99807  41.16267  514.5334  533.5656
##   Communes inchangées   791.00193 516.83733 6460.4666 6699.4344

tab.chi2$residuals

##                        
##                                 111         112         120         211
##   Communes fusionnantes -10.0921904  -2.3852264   1.0104527  -1.4733210
##   Communes inchangées     2.8481331   0.6731385  -0.2851615   0.4157883
##                        
##                                 212         221         222         300
##   Communes fusionnantes   1.9398795   4.6618769  -1.4281385   3.3269637
##   Communes inchangées    -0.5474565  -1.3156357   0.4030372  -0.9389077
##                        
##                                 400
##   Communes fusionnantes   4.3046989
##   Communes inchangées    -1.2148359

mosaicplot(ZAU_nb, shade = TRUE,  cex=0.6, main = NULL, xlab = "Communes fusionnantes",
           ylab = "Communes inchangées")

# Proportion sur l'ensemble des communes françaises
Ensemble <- apply(ZAU_nb, 2, sum)
ZAU_nb <- rbind(ZAU_nb, Ensemble)

ZAU_prop <- t(prop.table(ZAU_nb, 1) * 100)
Total <- apply(ZAU_prop, 2, sum)
ZAU_prop <- rbind (ZAU_prop, Total)

kable(ZAU_prop, row.names = T, digits = 2, caption = "Communes fusionnantes en fonction du ZAU")

Ces éléments permettent de se rendre compte que les profils des communes fusionnantes et des communes inchangées sont relativement proches. Les écarts les plus importants se situent aux extrêmes : ainsi, les communes fusionnantes sont marquées par une plus faible représentation des communes centres de grands ou moyens pôles dans les communes nouvelles et une sur-représentation des communes multipolarisées ou hors influence des pôles. Le test de Chi² réalisé permet de rejeter l’hypothèse d’indépendance des deux variables.

# Idée : faire la même chose avec variable "FusPhas" > est-ce que ces fusions ont significativement touché
# des communes rurales au début, à la fin... 
# A priori pas très concluant

test <- table(datafus2011$FusPhas, datafus2011$CATAEU2010)
test
prop <- prop.table(test, 1) * 100

barplot(prop, main= NULL,
        xlab = "Catégories d'aire urbaine",
        ylab = "Part du total des communes d'une phase (%)",
        las = 2,
        border = NA,
        col= rainbow(nrow(prop)),
        beside = TRUE)

legend(x="topright", legend = rownames(prop) , cex=0.8, fill= rainbow(nrow(prop)),bty="n")     


# Chi2
tab.chi2 <- chisq.test(test)
tab.chi2 

tab.chi2$observed
tab.chi2$expected
tab.chi2$residuals

mosaicplot(test, shade = TRUE,  cex=0.6, main = NULL, xlab = "Phase de fusion",
           ylab = "Catégorie d'aire urbaine")

# Comparaison en fonction de catégorisations de ZAU
# Analyse peu concluante

listeZAU <- c("CATAEU2010", "ZAU_POL", "ZAU_RUR", "ZAU_MAR_SP", "ZAU_MAR", "ZAU_PERI", "ZAU_AU")


results_chi2 <- data.frame()
for (i in listeZAU) { # Pour chaque identifiant de commune nouvelle, 
  a <- which(colnames(dataNfus2011) == i)
  ZAU_non <- data.frame(dataNfus2011[a])
  ZAU_oui <- data.frame(datafus2011[a])
  ZAU_non$Fusion <- "Pas de fusion"
  ZAU_oui$Fusion <- "Communes fusionnées"
  colnames(ZAU_non)[1] <- "CATAEU2010"
  colnames(ZAU_oui)[1] <- "CATAEU2010"
  ZAU <- rbind(ZAU_oui, ZAU_non)
  ZAU <- table(ZAU$Fusion, ZAU$CATAEU2010)

  tab.chi2 <- chisq.test(ZAU)
  Xsquared <- tab.chi2$statistic
  df <- tab.chi2$parameter
  p.value <- tab.chi2$p.value
  resultpartiel <- c(i, Xsquared, df, p.value)
  results_chi2 <- rbind (results_chi2, resultpartiel, stringsAsFactors = FALSE)
  rm(a, ZAU_non, ZAU_oui, ZAU, tab.chi2, p.value, df, Xsquared, resultpartiel)
  }

colnames(results_chi2) <- c("Catégorisation testée", "X-squared", "df", "p.value")

2.2. Chaque commune nouvelle est-elle composée d’un seul type de ZAU ?

# Définition d'un tableau comportant le nombre de communes fusionnantes par commune nouvelle
count_CN <- plyr::count(datafus2011, "CODGEO_new")

mesCommunes <- count_CN$CODGEO_new
#mesCommunes <- c("61324","73150","73006")

df <- datafus2011

results <- data.frame()
for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
  toto <- subset (df, CODGEO_new == i ) # Ne garder que les communes fusionnantes y ayant participé
  a <- table(toto$CATAEU2010) # Relever les ZAU des communes fusionnantes
  results <- rbind(results,a) # Combiner les résultats, par lignes
 # rm(a) # Supprimer "a"
}
x <- c("111","112","120","211","212","221","222","300","400") # On renomme les colonnes
colnames(results) <- x
count_CN <- cbind(count_CN, results)

# On identifie la variable la plus fréquente
count_CN$max <- apply(count_CN[, 3:11], 1, function(x) max(x, na.rm = TRUE))
# Quel est le code ZAU revenant le plus fréquemment dans une CN
count_CN$ZAUmaj2 <- colnames(count_CN[, 3:11])[apply(count_CN[, 3:11], 1, which.max)]
# On note si une CN a des communes fusionnantes avec un ZAU identique
count_CN$ZAUident <- ifelse(count_CN$max == count_CN$freq, TRUE, FALSE)
summary(count_CN$ZAUident)

##    Mode   FALSE    TRUE 
## logical     338     506

rm(results)
# On extrait les CN ayant des communes fusionnante avec un ZAU identique
CNIdent <- subset(count_CN, ZAUident==TRUE)
CNIdent <- merge(CNIdent, datafus2011[, c("CODGEO", "CATAEU2010")], by.x = "CODGEO_new", by.y = "CODGEO", all.x = TRUE)

tabcont<-summary(CNIdent$CATAEU2010)
round(100*prop.table(tabcont,margin=),1) # Pourcentages de chaque ZAU chez les communes nouvelles

##  111  112  120  211  212  221  222  300  400 NA's 
##  4.5 42.9  9.1  0.4  2.2  1.4  0.4 12.6 26.1  0.4

tabcont2<- summary(datafus2011$CATAEU2010)
round(100*prop.table(tabcont2,margin=),1) # Pourcentages de chaque ZAU chez les communes françaises

##  111  112  120  211  212  221  222  300  400 
##  3.0 31.0 11.6  0.9  2.8  3.7  1.2 22.1 23.7

# On peut donc en conclure qu'une large majorité des communes nouvelles sont, du point de vue du ZAU, composées de communes au profil homogène.

3.0 Préparation des données

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) # Les communes selon la géographie en vigueur au 1er janvier 2011.
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE) # Les  communes nouvelles, avec les géométries au 1er janvier 2025 et caractérisées par les données à la géométrie 2011 agrégées.
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE) # Départements


# Import et sélection des données
load("data/refdata.Rdata")

df2011 <- df2011 [, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "CATAEU2010", "FusPhas", "COM_NOUV", "FUSION", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT","C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "P09_POP75PY_RT", "P11_POT_FIN", "P09_POP", "REG")]
df_new <- df_new [, c("LIBGEO_new", "CODGEO_new", "CODGEO", "CATAEU2010", "FusPhas", "COM_NOUV", "FUSION", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT","C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "P09_POP75PY_RT", "P11_POT_FIN", "P09_POP", "REG")]

# Pour prendre les mêmes données que sur la CAH
# df_new <- df_new [, c("CODGEO","CODE_DEPT", "CATAEU2010", "LIBGEO_new", "P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT", "COM_NOUV", "CODGEO_new")]
# df2011 <- df2011 [, c("CODGEO","CODE_DEPT", "CATAEU2010", "LIBGEO_new", "P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT", "COM_NOUV", "CODGEO_new")]


datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes
dataCN_new <- subset(df_new, COM_NOUV == "OUI") # Les  communes nouvelles, avec les géométries au 1er janvier 2025 et caractérisées par les données à la géométrie 2011 agrégées.

# Appariement des données
geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")
geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")

3.1 Proximité statististique (distance euclidienne) entre communes fusionnantes

On calcule ici la proximité statistique qui existe, ou non, entre communes ayant fusionné au sein d’un même ensemble. À partir des données socioéconomiques, nous pouvons observer la proximité entre des communes (LIBGEO_2011) de même appartenance communale (LIBGEO_new).

tmp <- datafus2011 [, c("CODGEO_new", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT","C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "P09_POP75PY_RT")]
# tmp <- datafus2011 [, c("P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT", "COM_NOUV", "CODGEO_new")]
tmp <- na.omit(tmp)

d <- by(tmp, tmp$CODGEO_new, function(x)dist(x, method = "euclidean"))
df <- as.data.frame(sapply(d,mean))

df$CODGEO_new <- row.names(df)
colnames(df)[1] <- "dist_CN"

geomCN_new  <- merge(geomCN_new , df, by.x = "CODGEO_new", all.x = TRUE)

# Export des données sur la distance euclidienne
st_write(obj = df, dsn = "sorties/disteucl.gpkg", layer = "disteucl", delete_layer = TRUE, quiet = TRUE)

# Cartographie
par(mar=c(0,0,1.2,0))

plot(st_geometry(dep), col = NA)

choroLayer(x = geomCN_new , var = "dist_CN",
           method = "quantile", nclass = 6,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Distance statistique (euclidienne) moyenne par groupe de communes nouvelles",
           add = TRUE)

layoutLayer(title = "Distance euclidienne intra-groupe", theme = "red.pal",
            author = "G. Bideau, R. Ysebaert, 2026.",
            sources = "Sources : INSEE, IGN, 2026\nIndicateurs : P09_RETR1564_RT, C09_ACT1564_Agr_RT, C09_ACT1564_ArtCom_RT, C09_ACT1564_Cadr_RT, C09_ACT1564_ProfInt_RT, C09_ACT1564_Empl_RT,C09_ACT1564_Ouvr_RT, P09_POP0014Y_RT, P09_POP1529Y_RT, P09_POP3044Y_RT, P09_POP4559Y_RT, P09_POP6074Y_RT, P09_POP75PY_RT")

# Zoom sur la Normandie 
norm <- subset (geomCN_new , REG == "23" | REG == "25") # Si anciennes régions
# norm <- subset (geomCN_new , REG == "28")
bb <- st_bbox(norm)
bbnorm <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
norm <- norm[!is.na(norm$LIBGEO),]
plot(st_geometry(bbnorm), col = NA, border = NA)

plot(st_geometry(dep), add = TRUE)

propSymbolsChoroLayer(x = norm, var = "P09_POP", var2 = "dist_CN",
                      method = "quantile", nclass = 6,
                      col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
                      border = NA,
                      legend.var.pos = "left",
                      legend.title.cex = 0.7,
                      legend.var.title.txt = "Population totale",
                      legend.var2.pos = "topleft", legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Distance statistique (euclidienne)\nmoyenne par groupe\nde communes nouvelles",
                      add = TRUE)

layoutLayer(title = "Distance euclidienne intra-groupe",
            author = "G. Bideau, R. Ysebaert, 2026.",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2026.\nIndicateurs : ACT_AGR_RT, ACT_ART_RT,
ACT_INT_RT, ACT_EMP_RT, ACT_OUV_RT, ACT_RET_RT, 
ACT_NON_RT, POT_FIN, POP_0010Y_RT, POP_1117Y_RT, 
POP_2539Y_RT,  POP_4054Y_RT, POP_5564Y_RT,
POP_6579Y_RT, POP_80PY_RT")

# To do : Jouer peut être cette distance avec moins de variables pour mieux contrôler ce qui se passe : par variables démo, éco, emploi et voir si les natures des distances sont les mêmes

Une grande diversité de situations est visible concernant cette distance euclidienne, mais avec des tropismes régionaux intéressants. En particulier, en Bretagne et Pays-de-la-Loire, ce sont des communes nouvelles plutôt marquées par une faible distance euclidienne qui sont majoritaires. C’est le cas aussi, bien que moins nettement, en Savoie ou en Alsace. Certaines régions sont davantage marquées par des communes nouvelles davantage hétérogènes : les Pyrénées, le sud du Massif Central. Enfin, dans les régions très marquées par les fusions, on peut noter le cas de la Normandie, où les situations sont très variables.

GB : développement de l’analyse des distances euclidiennes à la fin de l’ACP/CAH.

4.1 Préparation des données et calcul des déviations

À partir du package MTA, on regarde la proximité statistique des communes nouvelles Cf. https://cran.r-project.org/web/packages/MTA/vignettes/MTA_Scenario.html (@ysebaert2019). On trouvera une vignette présentant le déroulé logique et la description précise des différentes fonctions du package ici.

Cette analyse porte ici sur le potentiel financier des communes.

On importe le ratio d’intérêt (potentiel financier par habitant par exemple), défini par un numérateur (potentiel financier) et un dénominateur (population municipale). On réalise cette extraction pour l’ensemble des communes françaises, afin de caractériser la position des communes nouvelles sur cet indicateur au regard de la moyenne nationale ; ainsi que juste pour les communes nouvelles, afin de pouvoir caractériser les inégalités existant au sein de celles-ci.

# Import des géométries
geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE)
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE)
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)

# Import des données d'intérêt : on ne garde que les variables utiles et
# celles qui peuvent être d'intérêt pour définir des appartenances
load("data/refdata.Rdata")

df2011 <- df2011 [, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "CATAEU2010", 
                      "P11_POT_FIN", "P11_Rev_Fisc", "P09_ACT1564", "P09_CHOM1564", "P09_POP", "REG", "CODE_DEPT","ChefLieu", "FUSION")]


# Jointures
geom2011 <- merge(geom2011, df2011, by = "CODGEO", all.x = TRUE)
geomfus2011 <- merge(geomfus2011, df2011, by = "CODGEO", all.x = TRUE)
geomCN_new <- merge(geomCN_new, df_new, by = "CODGEO_new", all.x = TRUE)

Grâce au package MTA, on calcule pour cet indicateur plusieurs déviations, dont la moyenne est exprimée en indice 100 : - Écart général à la moyenne française et à celle des communes nouvelles. - Écart aux communes de même catégorie de ZAU, pour l’ensemble des communes françaises et uniquement les communes nouvelles.

num <- "P11_POT_FIN"
nomnum <- "Potentiel financier"
denom <- "P09_POP"
# 
# num <-  "P11_Rev_Fisc"
# denom <- "P09_POP"
# geom2011 <- na.omit(geom2011)
# df2011 <- na.omit(df2011)
# 
num <-  "P09_CHOM1564"
nomnum <- "Taux de chômage des 15-64 ans\n"
denom <- "P09_ACT1564"

# Calcul des déviations générales
geom2011$gdevfr <- gdev(x = geom2011,  var1 = num, var2 = denom, type = "rel")

# Calcul des déviations à la ZAU de même appartenance
geom2011$mdevfr <- tdev(x = geom2011, var1 = num, var2 = denom, type = "rel",
                        key = "CATAEU2010")


# On ne s'intéresse qu'aux communes fusionnantes
df <- st_set_geometry(geom2011, NULL)
geomfus2011 <- merge(geomfus2011, df[,c("CODGEO", "gdevfr", "mdevfr")],
                     by = "CODGEO", all.x = TRUE)

4.2 Position des communes fusionnantes au regard de la moyenne française et de leur ZAU d’appartenance

On s’intéresse ici à caractériser globalement les communes fusionnantes au regard d’une variable (potentiel financier ou taux de chômage). Grâce au package MTA, nous commençons par réaliser deux représentations cartographiques. La première pour évaluer la position des communes nouvelles au regard de l’ensemble de ces communes.

La seconde pour observer uniquement la situation des communes fusionnantes entre elles. Ceci afin de mettre en évidence quelles sont celles qui disposent du plus de potentiel financier entre elles.

cette analyse permet d’identifier 3 cas d’étude qui semblent intéressants : Normandie et dep 49 (situations défavorables) ; communes nouvelles de Savoie (qui s’en tirent toutes bien) > ne serait-ce pas des communes abritant des stations de ski ?

# Visualisation des déviations générales et territoriales
cols <- carto.pal(pal1 = "blue.pal", n1 = 3, pal2 = "wine.pal", n2 = 3)

# plot a choropleth map of the relative global deviation
par(mfrow = c(1,2), mar = c(0,0,1.2,0))
choroLayer(x = geomfus2011, var = "gdevfr", legend.pos = "topleft",
           legend.title.txt = "Déviation relative\n(100 = moyenne française)",
           legend.title.cex = 0.7,
           breaks = c(min(geomfus2011$gdevfr, na.rm = TRUE),
                      75, 90, 100, 111, 133,
                      max(geomfus2011$gdevfr, na.rm = TRUE)), 
           border = NA, col = cols)

plot(st_geometry(dep), col = NA, add = TRUE)

layoutLayer(#title = "Potentiel financier\nécart à la moyenne française",
            title = paste0(nomnum, "\nécart à la moyenne française"),
            sources = "INSEE, 2009-2011.", author = "G. Bideau, R. Ysebaert, 2026.",
            scale = FALSE, tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black")
dev.print(device = svg, file = paste0("figures/Dev ", num, " ecart_Fr", ".svg"))

## png 
##   2

choroLayer(x = geomfus2011, var = "mdevfr", legend.pos = "topleft",
           legend.title.txt = "Déviation relative\n(100 = moyenne de la catégorie de ZAU d'appartenance)",
           legend.title.cex = 0.7,
           breaks = c(min(geomfus2011$mdevfr, na.rm = TRUE),
                      75, 90, 100, 111, 133,
                      max(geomfus2011$mdevfr, na.rm = TRUE)), 
           border = NA, col = cols)

plot(st_geometry(dep), col = NA, add = TRUE)

layoutLayer(#title = "Potentiel financier\nÉcart à la moyenne de la ZAU d'appartenance",
            title = paste0(nomnum, "\nécart à la moyenne française"),
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black")

dev.print(device = svg, file = paste0("figures/Dev ", num, " ecart_ZAU", ".svg"))

## png 
##   2

La fonction bidev présente la synthèse du positionnement des communes nouvelles sur ces deux déviations. Ce graphique reprend une des fonctionnalités d’HyperAtlas, le logiciel de référence historique qui permet l’analyse des inégalités territoriales. Cette fonction est désormais implémentée dans le package MTA.

Pour synthétiser la positionnement des communes nouvelles sur deux déviations, on commence par les situer dans un repère orthonormé : les communes représentées en rouge sont celles dont la moyenne du potentiel financier par habitant se situe au dessus de la moyenne française et de leur ZAU d’appartenance ; celles représentées en bleu en dessous de la moyenne pour les deux contextes ; et en jaune et vert dans des situations contradictoires en fonction des contextes spatiaux.

On s’intéresse aussi ici à l’éloignement statistique au regard des valeurs moyennes : plus la valeur des communes se situe loin des valeurs moyennes des contextes, plus la tonalité de couleur est saturée. Les communes qui se situent autour des valeurs moyennes (compris entre 75 et 125 sur les deux contextes, 100 représentant la moyenne) sont représentées en blanc.

Un graphique est associé à la représentation cartographique afin de visualiser le nuage de point sur lequel porte la représentation cartographique. Les communes de chaque même ZAU s’alignent alors selon un segment de droite dont la longueur exprime l’importance des disparités des communes qui ont fusionné et la pente la différence entre la moyenne de la ZAU et la moyenne nationale (@grasland2004) L’alignement des points selon différentes courbures de droites dans le repère permet ici d’apprécier l’effet de l’appartenance à une ZAU sur les inégalités induites par le potentiel financier. Basiquement, plus la pente du nuage de point de communes de même ZAU est importante, moins les régions qui la compose sont différenciées (les écarts intra-zone globaux étant faibles).

RY : référence de C.Grasland à rajouter en biblio sur ce point si on retient l’analyse. La lecture de cet article peut donner des billes pour interpréter correctement ces sorties d’ailleurs: Claude Grasland. Les inégalités régionales dans une Europe élargie. Les incertitudes du grand élargissement : L’Europe centrale et balte dans l’intégration européenne, L’Harmattan, pp.181-214, 2004, “ Pays de l’Est ”.

# Supprimer 0 (transformations log impossible)
geomfus2011 <- geomfus2011[geomfus2011$gdevfr != 0,]
geomfus2011 <- geomfus2011[geomfus2011$mdevfr != 0,]

# Calculer map_bidev
xx <- map_bidev(x = geomfus2011, dev1 = "gdevfr", dev2 = "mdevfr",
                breaks = c(75, 150, 300))

# Générer les géométries et le vecteur de couleur
com <- xx$geom
cols <- xx$cols

par(mfrow = c(1,2), mar = c(0,4,0,0))

typoLayer(x = com, var = "bidev", border = NA, col = cols, lwd = 0.2, 
          legend.pos = "n")
plot(dep$geom, col = NA, lwd = 1, add = TRUE)

layoutLayer(title = paste0("Synthèse des deux déviations.\nVariable : ", num),
            author = "Source : INSEE, 2026.", scale = 50, col = "white", coltitle = "black")

plot_bidev(x = com,  dev1 = "gdevfr",  dev2 = "mdevfr", 
           dev1.lab = "Déviation nationale",
           dev2.lab = "Déviation aux ZAU de même appartenance",
           cex.lab = 0.8, breaks = c(75, 150, 300), cex.pt = 0.2, cex.axis = 0.5)

4.3 Zoom

On prend ici 4 cas d’étude : communes de l’ancienne région Rhône-Alpes (qui présente des situations contrastées), avec un zoom sur sa partie Est (la plus touchée par les communes nouvelles), la région Normandie (dont le potentiel financier des communes est visiblement faible au regard contextes spatiaux identifiés en partie précédente), tout comme le département 49 (très dense en communes nouvelles).

# Extraction cas d'étude
norm <- subset (geom2011, REG == "23" | REG == "25") # Normandie
normCN <- subset (geomCN_new, REG == "28")

dep49 <- subset (geom2011, CODE_DEPT == "49" ) # Dep 49
dep49CN <- subset (geomCN_new, CODE_DEPT == "49")

ralp <- subset (geom2011, REG == "82") # Rhône-Alpes
ralpCN <- subset (geomCN_new, CODE_DEPT == "69" | CODE_DEPT == "73"
                  | CODE_DEPT == "74" | CODE_DEPT == "38"
                  | CODE_DEPT == "42" | CODE_DEPT == "01"
                  | CODE_DEPT == "26" | CODE_DEPT == "07")

ralppartiel <- subset (geom2011, CODE_DEPT == "69" | CODE_DEPT == "73"
                  | CODE_DEPT == "74" | CODE_DEPT == "38" | CODE_DEPT == "01") # Rhône-Alpes partiel
ralppartielCN <- subset (geomCN_new, CODE_DEPT == "69" | CODE_DEPT == "73"
                  | CODE_DEPT == "74" | CODE_DEPT == "38" | CODE_DEPT == "01")

# Gestion des couleurs et des labels
colsmst <- function(x){
  
  coldf <- data.frame(colvec = c("#f0f0f0", "#fdc785","#ffffab","#fba9b0",
                                 "#addea6","#ffa100","#fff226","#e30020"),
                      mst = seq(0,7,1),
                      leg_val = c("∅","R -    -   ","   - Z -   ","R - Z -   ","   -    - C", "R -    - C", "   - Z - C", "R - Z - C"), # ligne pour disposer les lettres de manière plus lisible
                      # leg_val = c("0","R","Z","R-Z","C", "R-C", "Z-C", "R-Z-C"), # ligne d'origine
                      stringsAsFactors = FALSE)
  
  xx <- st_set_geometry(x, NULL)
  xx <- coldf[coldf$mst %in% xx[,"mst"],]
  cols <- xx$colvec
  leg_val <- as.vector(xx$leg_val)
  
  return(list("cols" = cols,  "leg_val" = leg_val))
}

# Calculer les trois déviations
norm$gdev <- gdev(x = norm, var1 = num, var2 = denom)
norm$tdev <- tdev(x = norm, var1 = num, var2 = denom, key = "CATAEU2010")
norm$sdev <- sdev(x = norm, var1 = num, var2 = denom, order = 1)

# Situation supérieur à 125 (souvent favorable)
mst <- map_mst(x = norm, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 125, superior = TRUE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
fuz <- st_set_geometry(geomfus2011, NULL)
fuz$fuz <- 1
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]

# Aire d'étude pour faciliter la lisibilité
study_area <-subset (dep, CODE_DEPT == "27" | CODE_DEPT == "50"
                  | CODE_DEPT == "14" | CODE_DEPT == "61" | CODE_DEPT == "76" | CODE_DEPT == "35" )
plot(st_geometry(study_area), border = NA)

typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n", add = TRUE)

depnorm <- subset (dep, CODE_DEPT == "27" | CODE_DEPT == "50"
                  | CODE_DEPT == "14" | CODE_DEPT == "61" | CODE_DEPT == "76")
plot(st_geometry(depnorm), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nau-dessus de l'indice 125 :",
           nodata = FALSE, pos = c(3410000, 2930000), title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation favorable (Normandie)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "center",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "NORM ", num, " sit_fav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
norm <- merge (norm, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(norm)[length(norm)-1] <- "mst_sup" # Classification selon déviation supérieur à 125


# Situation inférieur à 75 (souvent défavorable)
mst <- map_mst(x = norm, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 75, superior = FALSE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
plot(st_geometry(study_area), border = NA)

typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n", add = TRUE)

plot(st_geometry(depnorm), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nen-dessous de l'indice 75 :",
           nodata = FALSE, pos =  c(3410000, 2940000), title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation défavorable (Normandie)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "center",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "NORM ", num, " sit_defav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
norm <- merge (norm, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(norm)[length(norm)-1] <- "mst_inf" # Classification selon déviation inférieur à 75

# Calculer les trois déviations
dep49$gdev <- gdev(x = dep49, var1 = num, var2 = denom)
dep49$tdev <- tdev(x = dep49, var1 = num, var2 = denom, key = "CATAEU2010")
dep49$sdev <- sdev(x = dep49, var1 = num, var2 = denom, order = 1)

# Situation supérieur à 125 (souvent favorable)
mst <- map_mst(x = dep49, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 125, superior = TRUE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
fuz <- st_set_geometry(geomfus2011, NULL)
fuz$fuz <- 1
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

depfig <- subset (dep, CODE_DEPT == "49")
plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nau-dessus de l'indice 125 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation favorable (Maine-et-Loire)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "49 ", num, " sit_fav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
dep49 <- merge (dep49, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(dep49)[length(dep49)-1] <- "mst_sup" # Classification selon déviation supérieur à 125


# Situation inférieur à 75 (souvent défavorable)
mst <- map_mst(x = dep49, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 75, superior = FALSE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nen-dessous de l'indice 75 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation défavorable (Maine-et-Loire)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "49 ", num, " sit_defav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
dep49 <- merge (dep49, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(dep49)[length(dep49)-1] <- "mst_inf" # Classification selon déviation inférieur à 75

# Calculer les trois déviations
ralp$gdev <- gdev(x = ralp, var1 = num, var2 = denom)
ralp$tdev <- tdev(x = ralp, var1 = num, var2 = denom, key = "CATAEU2010")
ralp$sdev <- sdev(x = ralp, var1 = num, var2 = denom, order = 1)

# Situation supérieur à 125 (souvent favorable)
mst <- map_mst(x = ralp, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 125, superior = TRUE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
fuz <- st_set_geometry(geomfus2011, NULL)
fuz$fuz <- 1
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

depfig <- subset (dep, CODE_DEPT == "69" | CODE_DEPT == "73"
                  | CODE_DEPT == "74" | CODE_DEPT == "38"
                  | CODE_DEPT == "42" | CODE_DEPT == "01"
                  | CODE_DEPT == "26" | CODE_DEPT == "07")
plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nau-dessus de l'indice 125 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation favorable (Rhône-Alpes)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "RALP ", num, " sit_fav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
ralp <- merge (ralp, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(ralp)[length(ralp)-1] <- "mst_sup" # Classification selon déviation supérieur à 125


# Situation inférieur à 75 (souvent défavorable)
mst <- map_mst(x = ralp, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 75, superior = FALSE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nen-dessous de l'indice 75 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation défavorable (Rhône-Alpes)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "RALP ", num, " sit_defav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
ralp <- merge (ralp, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(ralp)[length(ralp)-1] <- "mst_inf" # Classification selon déviation inférieur à 75

# Calculer les trois déviations
ralppartiel$gdev <- gdev(x = ralppartiel, var1 = num, var2 = denom)
ralppartiel$tdev <- tdev(x = ralppartiel, var1 = num, var2 = denom, key = "CATAEU2010")
ralppartiel$sdev <- sdev(x = ralppartiel, var1 = num, var2 = denom, order = 1)

# Situation supérieur à 125 (souvent favorable)
mst <- map_mst(x = ralppartiel, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 125, superior = TRUE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
fuz <- st_set_geometry(geomfus2011, NULL)
fuz$fuz <- 1
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

depfig <- subset (dep, CODE_DEPT == "69" | CODE_DEPT == "73"
                  | CODE_DEPT == "74" | CODE_DEPT == "38" | CODE_DEPT == "01")
plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nau-dessus de l'indice 125 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 4),]# Sélection, au hasard, de 4 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation favorable (Rhône-Ain-Isère-Savoies)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "RALP-Partiel ", num, " sit_fav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
ralppartiel <- merge (ralppartiel, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(ralppartiel)[length(ralppartiel)-1] <- "mst_sup" # Classification selon déviation supérieur à 125


# Situation inférieur à 75 (souvent défavorable)
mst <- map_mst(x = ralppartiel, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev",
               threshold = 75, superior = FALSE)

# Unlist outputs of the function
com <- mst$geom

# Ne garder que les communes fusionnantes
com <- merge(com, fuz[,c("CODGEO", "fuz")], by.x = "CODGEO",
             by.y = "CODGEO", all.x = TRUE)
com <- com[!is.na(com$fuz),]
com <- com[order(com$mst),]

# Légende et couleurs
xx <- colsmst(x = com)
cols <- xx$cols
leg_val <- xx$leg_val

# Cartographie
par(mar = c(0,0,0,0), mfrow = c(1,1))
bb <- st_bbox(com)
bb <- st_as_sfc(st_bbox(bb + c(-60000,0,0,0), crs = 3035))
com <- com[!is.na(com$LIBGEO),]
typoLayer(x = com, var = "mst", border = "grey50",
          col = cols, lwd = 0.2, legend.pos = "n")

plot(st_geometry(depfig), lwd = 1, col = NA, add = TRUE)

legendTypo(col = cols, categ = leg_val,
           title.txt = "Communes fusionnantes\nen-dessous de l'indice 75 :",
           nodata = FALSE, pos = "topleft", title.cex = .7, values.cex = .6)

toponymes <- com[com$mst == 7,] # Sélection des noms des communes correspondant à la classe 7
toponymes <- toponymes[sample(1:nrow(toponymes), 6),]# Sélection, au hasard, de 6 toponymes
labelLayer(toponymes, txt = "LIBGEO", halo = TRUE, overlap = FALSE, cex = 0.6)

layoutLayer(# title = paste0(nomnum, " : Communes nouvelles en situation défavorable (Rhône-Ain-Isère-Savoies)"),
            title = paste0("\n\n\n\n\n\n", nomnum), postitle = "right",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            source = "Source : INSEE, 2026.", author =  "G. Bideau, R. Ysebaert, 2026.
100: Moyenne des déviations
R: Situation au regard de la moyenne des communes de l'espace d'étude (région)
Z: Situation au regard de la moyenne des ZAU d'appartenance sur l'espace d'étude.
C: Situation au regard des communes contigues")

# dev.print(device = svg, file = paste0("figures/Dev ", "RALP-Partiel ", num, " sit_defav", ".svg"))

# Intégration des données de déviation pour sauvegarde
st_geometry(com) <- NULL # On supprime les géométries car sinon fusion impossible
ralppartiel <- merge (ralppartiel, com[, c("CODGEO", "mst")], by = "CODGEO", all.x = TRUE)
colnames(ralppartiel)[length(ralppartiel)-1] <- "mst_inf" # Classification selon déviation inférieur à 75

st_write(obj = norm [, c("CODGEO", "LIBGEO", "CODGEO_new", "LIBGEO_new", "mst_sup", "mst_inf")], dsn = paste0("sorties/norm_", num, ".gpkg"), layer = "norm", delete_layer = TRUE, quiet = TRUE)
st_write(obj = dep49 [, c("CODGEO", "LIBGEO", "CODGEO_new", "LIBGEO_new", "mst_sup", "mst_inf")], dsn = paste0("sorties/dep49_", num, ".gpkg"), layer = "dep49", delete_layer = TRUE, quiet = TRUE)
st_write(obj = ralp [, c("CODGEO", "LIBGEO", "CODGEO_new", "LIBGEO_new", "mst_sup", "mst_inf")], dsn = paste0("sorties/ralp_", num, ".gpkg"), layer = "ralp", delete_layer = TRUE, quiet = TRUE)



# Pour importer les données concernant le taux de chômage
norm_P09_CHOM1564 <- st_read("sorties/norm_P09_CHOM1564.gpkg", layer = "norm", quiet = TRUE) 
dep49_P09_CHOM1564 <- st_read("sorties/dep49_P09_CHOM1564.gpkg", layer = "dep49", quiet = TRUE) 
ralp_P09_CHOM1564 <- st_read("sorties/ralp_P09_CHOM1564.gpkg", layer = "ralp", quiet = TRUE) 

# Pour importer les données concernant le potentiel financier
norm_P11_POT_FIN <- st_read("sorties/norm_P11_POT_FIN.gpkg", layer = "norm", quiet = TRUE) 
dep49_P11_POT_FIN <- st_read("sorties/dep49_P11_POT_FIN.gpkg", layer = "dep49", quiet = TRUE) 
ralp_P11_POT_FIN <- st_read("sorties/ralp_P11_POT_FIN.gpkg", layer = "ralp", quiet = TRUE) 

dep49_P11_POT_FIN <- as.data.frame(dep49_P11_POT_FIN) # Pour que le sf devienne un simple data frame (sinon, on peut aussi utiliser st_geometry(data) <- NULL qui retire/supprime la géométrie
dep49Croisement <- merge (dep49_P09_CHOM1564, dep49_P11_POT_FIN[, c("CODGEO", "mst_sup", "mst_inf")], by = "CODGEO")

colnames(dep49Croisement) <- c("CODGEO", "LIBGEO", "CODGEO_new", "LIBGEO_new", "mst_sup_CHOM", "mst_inf_Chom", "mst_sup_PotFin", "mst_inf_PotFin","geometry")

summary(dep49Croisement$mst_inf_Chom >0)

# Définition de la variable à étudier
geom2011$variableproxim <- geom2011$P11_POT_FIN
# Définition de l'ensemble de référence choisi
# Choisir "Fusion" si on souhaite comparer à l'ensemble des communes fusionnantes, "REG" aux régions etc
geom2011$CATAEUFusion <- interaction (geom2011$CATAEU2010, geom2011$FUSION)
geom2011$key <- geom2011$CATAEUFusion
geom2011$Reference <- geom2011$P09_POP

# Déviation globale
geom2011$gdevrel <- gdev(x = geom2011, 
                         var1 = "variableproxim",
                         var2 = "Reference",
                         type = "rel")

# general absolute deviation 
geom2011$gdevabs <- gdev(x = geom2011, 
                         var1 = "variableproxim", 
                         var2 = "Reference", 
                         type = "abs")

# general deviation in million Euros
geom2011$gdevabsmil <- geom2011$gdevabs / 1000000

####  Cartography ----
prcarto<-subset(geom2011, FUSION == "OUI")
# margins
#par(mar = c(0, 0, 1.2, 0))

# Plot territories
plot(st_geometry(geom2011), col = "grey70", border = "#EDEDED", lwd = 0.25)
#plot(st_geometry(ept), border = "#1A1A19", lwd = 1, add = TRUE)

# Global deviation (relative and absolute) cartography
propSymbolsChoroLayer(x = prcarto,
                      var = "gdevabsmil", var2 = "gdevrel",
                      add = TRUE,
                      inches = 0.1,
                      col = carto.pal(pal1 = "blue.pal", n1 = 3,
                                      pal2 = "wine.pal", n2 = 3),
                      breaks = c(min(prcarto$gdevrel, na.rm = TRUE),
                                 75, 90, 100, 111, 133,
                                 max(prcarto$gdevrel, na.rm = TRUE)),
                      #border = "#f0f0f0",
                      #lwd = 0.25,
                      legend.var.pos = "left", legend.var2.pos = "topleft",
                      legend.var.title.txt = "Potentiel financier (millions d'euros)",
                      legend.var2.title.txt = "Déviation par rapport au contexte global (100 = Ensemble des communes françaises)",
                      legend.var.style = "e",
                      legend.var.values.rnd = 0,
                      legend.var2.values.rnd = 0)

# layout 
layoutLayer(title = "Global deviation - Potentiel financier",
            sources = "INSEE, 2009-2011.",
            north = TRUE, scale = 5, tabtitle = TRUE, frame = FALSE, theme = "red.pal",
            author = "!!!")

# Territorial relative deviation calculation
geom2011$mdevrel <- tdev(x = geom2011, 
                         var1 = "variableproxim", 
                         var2 = "Reference", 
                         type = "rel",
                         key = "key")

# Territorial absolute deviation calculation
geom2011$mdevabs <- tdev(x = geom2011, 
                         var1 = "variableproxim", 
                         var2 = "Reference", 
                         type = "abs",
                         key = "key")

# Territorial deviation in million Euros
geom2011$mdevabsmil <- geom2011$mdevabs / 1000000

####  Cartography ----
# Plot layout
par(mfrow = c(1, 1))

prcarto<-subset(geom2011, FUSION == "OUI")
prcarto<- na.omit(prcarto)
# Plot territories
plot(st_geometry((geom2011)), col = "grey70", border = "#EDEDED",lwd = 0.25)
plot(st_geometry((dep)), border = "#1A1A19", lwd = 1)

# Territorial deviation (relative and absolute) cartography
propSymbolsChoroLayer(x = prcarto,
                      var = "mdevabsmil", var2 = "mdevrel",
                      add = TRUE,
                      inches = 0.3,
                      col = carto.pal(pal1 = "blue.pal", n1 = 3,
                                      pal2 = "wine.pal", n2 = 3),
                      breaks = c(min(prcarto$mdevrel, na.rm = TRUE),
                                 75, 90, 100, 111, 133,
                                 max(prcarto$mdevrel, na.rm = TRUE)),
                      border = "#f0f0f0",
                      lwd = 0.25,
                      legend.var.pos = "left", legend.var2.pos = "topleft",
                      legend.var.title.txt = "Potentiel financier (millions d'euros)",
                      legend.var2.title.txt = "Déviation par rapport à l'ensemble régional (100 = moyenne des communes de la région)",
                      legend.var.style = "e",
                      legend.var.values.rnd = 0,
                      legend.var2.values.rnd = 0)


# layout 
layoutLayer(title = "!!!!",
            sources = "!!!!",
            north = TRUE, scale = 5, tabtitle = TRUE, frame = FALSE, theme = "red.pal",
            author = "!!!!")

####  Réalisation d'un box-plot pour faciliter la visualisation ----
par(cex.lab = 1)
par(cex.axis = 0.75)
par(mar = c(4, 4, 2, 2))

# Drop geometries
df <- st_set_geometry(geom2011, NULL)
# Reorder EPT according to gdev value
df$key <- with(df, reorder(key, gdevrel, mean, na.rm = TRUE))
df$key <- with(df, reorder(key, CATAEUFUSION, na.rm = TRUE))

# Colors management
col <- carto.pal(pal1 = "red.pal", n1 = (nlevels(df$key) / 2), 
                 pal2 = "green.pal", n2 = (nlevels(df$key) / 2),
                 middle = FALSE, transparency = TRUE)
library(ggplot2)
ggplot(df, aes(x=CATAEU2010, y=gdevrel, fill=FUSION)) +
  geom_violin() +
  geom_boxplot()
toto<- subset(df, gdevrel >= 600)

# Boxplot
boxplot(df$gdevrel ~ df$key,
        col = col,
        ylab = "Global deviation",
        xlab = "Territorial deviation",
        varwidth = TRUE,
        range = 1,
        outline = TRUE,
        las = 1) 

# Horizontal Ablines
abline (h = seq(40, 300, 10), col = "#00000060", lwd = 0.5, lty = 3)

# Plot mean values
xi<- tapply(df$gdevrel, df$key, mean, na.rm = TRUE)
points(xi, col = "navy", pch = 19)

head(df)

# Plot mean values (Communes nouvelles)
dfCfus <- subset(df, COM_NOUV == "OUI")
xi<- tapply(dfCfus$gdevrel, dfCfus$key, mean, na.rm = TRUE)
points(xi, col = "#7C0000", pch = 19)


# Legend for the boxplot
# df$EPTName<- as.factor(df$LIBEPT)
# df$EPTName <- with(df, reorder(EPTName, gdevrel, mean, na.rm = TRUE))
# legend("topleft",
#        legend = levels(df$REG),
#        pch = 15,
#        col = col,
#        cex = 0.6,
#        pt.cex = 1,
#        title = "Territorial contexts (ordered by mean value of global deviation)")

# Spatial relative deviation calculation
geom2011$ldevrel <- sdev(x = geom2011, xid = "CODGEO", var1 = "variableproxim", var2 = "P09_POP",
                         order = 1, type = "rel")


# Spatial absolute deviation calculation
geom2011$ldevabs <- sdev(x = geom2011, xid = "CODGEO", var1 = "variableproxim", var2 = "P09_POP",
                         order = 1, type = "abs")

# Spatial deviation in million Euros
geom2011$ldevabsmil <- geom2011$ldevabs / 1000000

# Cartography
# Plot layout
par(mfrow = c(1, 1), mar = c(0, 0, 1.2, 0))

# Plot territories
plot(st_geometry(geom2011), col = "grey70", border = "#EDEDED",lwd = 0.25)
plot(st_geometry(dep), border = "#1A1A19",lwd = 1, add = T)

# Territorial deviation (relative and absolute) cartography
propSymbolsChoroLayer(x = geom2011,
                      var = "ldevabsmil", var2 = "ldevrel",
                      add = TRUE,
                      inches = 0.3,
                      col = carto.pal(pal1 = "blue.pal", n1 = 3,
                                      pal2 = "wine.pal", n2 = 3),
                      breaks = c(min(geom2011$ldevrel, na.rm = TRUE),
                                 75, 90, 100, 111, 133,
                                 max(geom2011$ldevrel, na.rm = TRUE)),
                      border = "#f0f0f0",
                      lwd = 0.25,
                      legend.var.pos = "left", legend.var2.pos = "topleft",
                      legend.var.title.txt = "Redistribution (Million euros)",
                      legend.var2.title.txt = "Deviation to the spatial context 
                      (100 = average of the contiguous territorial units - order 1)",
                      legend.var.style = "e",
                      legend.var.values.rnd = 0,
                      legend.var2.values.rnd = 0)

# layout 
layoutLayer(title = "Spatial deviation - xxxxxxxxxxxxxxx",
            sources = "x.xxxxx",
            north = TRUE, scale = 5, tabtitle = TRUE, frame = FALSE, theme = "red.pal",
            author = "xxxx")

# Drop geometries
df <- st_set_geometry(geom2011, NULL)

# Spatial deviation - Top 10 of the potential contributors as regards to their total amount of income
df$ldevabsPerc<- df$ldevabs / df$variableproxim * 100
df<-df[order(df$ldevabsPerc, decreasing = TRUE), ]
df[1:10, c("ldevabsmil","ldevabsPerc")]

# Spatial deviation - Top 10 of the potential receivers as regards to their total amount of income
df<-df[order(df$ldevabsPerc, decreasing = FALSE), ]
df[1:10, c("ldevabsmil","ldevabsPerc")]

5.0 Préparation des données

geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE) # Départements

# Import des données
load("data/refdata.Rdata")
rm(df_new)


df2011 <- df2011[, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "CATAEU2010", "FusPhas", "COM_NOUV", "FUSION",  "P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT", "P11_POT_FIN", "P09_POP", "REG", "CODE_DEPT", "P11_POT_FIN_RT", "P09_POP")]
datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes

#Calcul de la surface si besoin
# datafus2011$surface <- st_area(geomfus2011) # Attention, unités : m²
# datafus2011$surface <- set_units(datafus2011$surface, km^2) # On passe en km²

# Appariement des données
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")

# Import données pour les noms de variables
ratio <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ratios"))
row.names(ratio) <- ratio$Numerator_Code

5.1 Sélectionner les variables d’intérêt

selecVar <- c("CODGEO","CODE_DEPT", "CATAEU2010", "LIBGEO_new", "P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT")

# test en intégrant la population et la surface
# selecVar <- c("CODGEO","CODE_DEPT", "CATAEU2010", "LIBGEO_new", "P09_CHOM1564_RT", "P09_ETUD1564_RT", "P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "C09_ACT1564_ArtCom_RT", "C09_ACT1564_Cadr_RT", "C09_ACT1564_ProfInt_RT", "C09_ACT1564_Empl_RT", "C09_ACT1564_Ouvr_RT", "P09_POP0014Y_RT", "P09_POP1529Y_RT", "P09_POP3044Y_RT", "P09_POP4559Y_RT", "P09_POP6074Y_RT", "C09_ACTOCC_OUT_RT", "P09_POP", "surface")

datafus2011 <- datafus2011[,selecVar]

PourCAH <- na.omit(datafus2011)
row.names(PourCAH) <- PourCAH$CODGEO
PourCAH$CODGEO <- NULL
PourCAH$CODE_DEPT <- NULL
PourCAH$CATAEU2010 <- NULL
PourCAH$LIBGEO_new  <- NULL

# PourCAH <- scale(PourCAH) # La standardisation peut se faire après

NbrVariables <- ncol(PourCAH)

selecVarCAH <- colnames(PourCAH)


# Si problèmes de doublons
PourCAH$CODGEO[duplicated(PourCAH$CODGEO)==TRUE]

## NULL

5.2 Réalisation de la typologie et valeurs

# On commence par réaliser une ACP (conseil donné dans R et espace, chapitre 7, p. 135)
res.pca <- PCA(PourCAH, graph = FALSE)
#res.pca <- PCA(PourCAH, ncp = 5, graph = FALSE) # 50 % of cumulative variance

# res.hcpc <- HCPC(res.pca, graph = TRUE, method = "ward", metric = "euclidian") #Ici coupure à 8 classes pas trop mal
#res.hcpc <- HCPC(res.pca, graph = FALSE, nb.clust = 7, method = "ward", metric = "euclidian") #Pour faire la coupure de manière déterministe

res.hcpc <- agnes(res.pca$ind$coord, metric = "euclidiean", method = "ward")


# Réalisation du dendrogramme en passant par plot
dendro.hcpc <- as.dendrogram(res.hcpc)
plot(dendro.hcpc, leaflab = "none", ylab = "Dissimilarité")

inertie <- sort(res.hcpc$height, decreasing = TRUE)
plot(inertie[1:20], type = "s", xlab = "Nombre de classes", ylab = "Inertie")

# Affichage de différents graphiques pour aider au découpage
sortedHeight <- sort(res.hcpc$height^2, decreasing = TRUE)
plot(sortedHeight,
     type = "h",
     xlab = "Noeuds",
     ylab = "Niveau d'agrégation")

relHeight <- sortedHeight / sum(sortedHeight) * 100
cumHeight <- cumsum(relHeight)
barplot(relHeight[1:30], names.arg = seq(1, 30, 1),
        col = "black", border = "white", xlab = "Noeuds",
        ylab = "Part de l'inertie totale (%)",
        main = "Diagramme de niveaux")

# Coupure de l'arbre (k = nombre de classes)
nclass <- 5
cluspop <- cutree(res.hcpc, k = nclass)

# cluspop

5.2.1 Valeurs absolues

# Définition de noms de communes
# Noms pour typo sans les navetteurs, avec 7 classes
# NomsGroupesCAH <- c("Communes d'urbains", "Communes familiales",
#                  "Communes au profil moyen, tendance laborieuse", "Communes rurales actives",
#                  "Communes jeunes et industrielles", "Communes agricoles vieillissantes",
#                  "Communes industrielles vieillissantes" )
# Noms pour typo avec les navetteurs, avec 5 classes
# NB : NOMS ADAPTÉS POUR LES DONNÉES 2011-2020 UNIQUEMENT
NomsGroupesCAH <- c("Groupe 1 PPCQ", # Proche Périphérie, Cadres, Quinquagénaires
                    "Groupe 2 PFO", # Périphéries Familiales et Ouvrières
                    "Groupe 3 MOARO", # Moyen, Ouvriers, Agriculteurs, Retraités, Out.
                    "Groupe 4 VAR", # Vieillissantes, Agricoles et Rurales
                    "Groupe 5 REAV" # Rurales, Enclavées, Agricoles et Vieillissantes
)
# Si on ne veut pas définir de noms
# NomsGroupesCAH <- paste ("Groupe",1 :nclass)

# On intègre ces données dans le tableau de départ
PourCAH <- as.data.frame(PourCAH, stringsAsFactors = FALSE)
PourCAH$Groupes <- factor(cluspop,
                          levels = 1:nclass,
                          labels = paste(NomsGroupesCAH))

# On identifie le nombre de chaque groupe
Groupes_nbr <- paste0(NomsGroupesCAH, " (n=", table(PourCAH$Groupes), ")")
names(Groupes_nbr) <- NomsGroupesCAH

PourCAH$CODGEO <- row.names(PourCAH)
# Calcul de la moyenne des variables
clusProfile <- aggregate(PourCAH [, 1:NbrVariables],
                         by = list(PourCAH$Groupes),
                         mean)
colnames(clusProfile)[1] <- "Groupes"
clusLong <- melt(clusProfile, id.vars = "Groupes")

ggplot(clusLong) +
  geom_bar(aes(x = variable, y = value, fill = Groupes),
           stat = "identity") +
  # scale_fill_grey() +
  # scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628","#f781bf")) + # Pour 8
   scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628")) + # Pour 7
  facet_wrap(~ Groupes, labeller = labeller(Groupes = Groupes_nbr)) +
  coord_flip() +
  theme_bw()

ggplot(clusLong) +
  geom_bar(aes(x = variable, y = value, fill = Groupes),
           stat = "identity") +
  # scale_fill_grey() +
  # scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628","#f781bf")) + # Pour 8
   scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628")) + # Pour 7
  facet_wrap(~ Groupes) +
  coord_flip() +
  theme_bw()

5.2.2 Valeurs moyennes

# Calculs des valeurs moyennes (pour comparaison)
# On identifie les variables en stock
VarCAHBrutes <- stringr::str_replace(selecVarCAH, "Y_RT", "")
VarCAHBrutes <- stringr::str_replace(VarCAHBrutes, "_RT", "")
# Import données totales
load("data/refdata.Rdata")
rm(df_new)
# Sélection données utiles
pourmoyennesCAH <- na.omit(subset(df2011, COM_NOUV == "OUI"))
moyennesCAH <- data.frame()
for (i in VarCAHBrutes){
  a <- ratio[i, "Denominator_Code"]
  b <- ratio[i, "Coeff"]
  c <- sum(df2011[, i], na.rm = TRUE)
  d <- sum(df2011[, a], na.rm = TRUE)
  moyfr <- round(b * c/d, 2)
  e <- sum(pourmoyennesCAH[, i], na.rm = TRUE)
  f <- sum(pourmoyennesCAH[, a], na.rm = TRUE)
  moyCFus <- round(b * e/f, 2)
  g <- ratio[i, "CODE"]
  h <- c(g, moyfr, moyCFus)
  moyennesCAH <- rbind(moyennesCAH, h, stringsAsFactors= FALSE)
  rm( a, b, c, d, e, f, g, h, moyfr, moyCFus)
}

colnames(moyennesCAH) <- c("Variable", "France", "CommunesFusionnantes")
moyennesCAH$Variable <- as.factor(moyennesCAH$Variable)
moyennesCAH$France <- as.numeric(moyennesCAH$France)
moyennesCAH$CommunesFusionnantes <- as.numeric(moyennesCAH$CommunesFusionnantes)
moyennesCAH$DiffComFusComFr <- moyennesCAH$CommunesFusionnantes - moyennesCAH$France

aa<- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = France), stat = "identity") + coord_flip()
bb <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = CommunesFusionnantes), stat = "identity") + coord_flip()
cc <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = DiffComFusComFr), stat = "identity") + coord_flip()
cowplot::plot_grid(aa, bb, cc, ncol = 1, nrow = 3)

rm(aa,bb,cc)

5.2.3 Valeurs standardisées

# On intègre ces données dans le tableau de départ

PourCAHz <- scale(PourCAH[,c(1:NbrVariables)])
PourCAHz <- as.data.frame(PourCAHz, stringsAsFactors = FALSE)

PourCAHz$Groupes <- PourCAH$Groupes

PourCAHz$CODGEO <- row.names(PourCAHz)
# Calcul de la moyenne des variables
clusProfileStd <- aggregate(PourCAHz [, 1:NbrVariables],
                         by = list(PourCAH$Groupes),
                         mean)
colnames(clusProfileStd)[1] <- "Groupes"
clusLongStd <- melt(clusProfileStd, id.vars = "Groupes")
clusLongStd <- merge(clusLongStd, ratio[, 1:2], by.x = "variable", by.y = "CODE") # Pour lisibilité du graphique
colnames(clusLongStd)[4] <- "Variable"

ggplot(clusLongStd) +
  geom_bar(aes(x = Variable, y = value, fill = Groupes),
           stat = "identity") +
  #  scale_fill_grey() +
  scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628","#f781bf")) + # Pour 8
  #  scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628")) + # Pour 7
  facet_wrap(~ Groupes, labeller = labeller(Groupes = Groupes_nbr)) +
    labs (y = "Valeurs standardisées", x = " ") +
  coord_flip() + theme_bw()

ggplot(clusLongStd) +
  geom_bar(aes(x = Variable, y = value, fill = Groupes),
           stat = "identity") +
  #  scale_fill_grey() +
  scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628","#f781bf")) + # Pour 8
  #  scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628")) + # Pour 7
  facet_wrap(~ Groupes) +
    labs (y = "Valeurs standardisées", x = " ") +
  coord_flip() + theme_bw()

5.2.4 Cartographie

# Remettre groupes dans geomfus2011 pour cartographie
typo <- merge(geomfus2011,PourCAH[ , c("CODGEO","Groupes")], by = "CODGEO")
## Carte Typoe nationale des communes fusionnées


# mf_export(x = typo, export = "svg", filename = "figures/TypoEspaceGeo.svg", 
#           width = 5, theme = "nevermind") # Si souhait d'export

par(mfrow = c(1, 1))
par(mar=c(0,0,1.2,0))
plot(st_geometry(dep), border = "#1A1A19",lwd = 1)

typoLayer(x = typo, var = "Groupes",  
          col = c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628","#f781bf"),
          border = NA, 
          #legend.title = "Typologie des communes fusionnantes/nà partir de critères socio-économiques",
          legend.title.cex = 0.7,
          legend.values.cex = 0.6,
          legend.pos = "left", add = T)

layoutLayer(title = "Typologie exploratoire des communes fusionnantes (2012-2026)",
            author = "G. Bideau, R. Ysebaert, 2026.",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2026.")

# dev.off() # Si souhait d'export

5.3 Comparaison et analyse des groupes

Dans le but de comparer un profil avec la moyenne française, on extrait les données de chaque groupe défini par la CAH (données issues de clusLong et clusProfile) et on les compile dans un tableau contenant déjà les données des moyennes françaises et des communes fusionnantes (moyennesCAH). En sortie, un tableau est créé pour chaque groupe.

# Tableau pour comparaison entre les groupes
clusProfile2 <- as.data.frame(t(clusProfile))
colnames(clusProfile2) <- clusProfile2[1,]
clusProfile2 <- clusProfile2[ c(2:nrow(clusProfile2)),]
compare <- cbind(moyennesCAH, clusProfile2)

for (i in NomsGroupesCAH){
  a <- which(colnames(compare)== i)
  b <- as.numeric(compare[,a])
  compare[,a] <- as.numeric(compare[,a])
  compare[paste0(i,"_DiffAvecFrance")] <- b - compare$France
  compare[paste0(i,"_DiffAvecComFus")] <- b - compare$CommunesFusionnantes
  rm(a, b)
}

kable(compare, row.names = F, digits = 2, caption = "Tableau de comparaison entre les groupes")

comparaisons <- data.frame()
for (i in NomsGroupesCAH){
  a <- subset(typo, typo$Groupes == i)
  Nb <- nrow(a)
  b <- round(mean(a$P11_POT_FIN), 2)
  bb <- round(mean(a$P11_POT_FIN_RT), 2)
  c <- table(a$CATAEU2010)
  h <- c(i, Nb, b, bb, c)
  comparaisons <- rbind(comparaisons, h, stringsAsFactors= FALSE)
  rm( a, b, bb, c, h, Nb)
} 
colnames(comparaisons) <- c("Groupes", "Nombre", "Moy_P11_POT_FIN", "Moy_P11_POT_FIN_RT", "111","112","120","211","212","221","222","300","400")
mean(typo$P11_POT_FIN)

mean(typo$P11_POT_FIN_RT)

tmp <- comparaisons[, c("Groupes", "Nombre", "111","112","120","211","212","221","222","300","400")]
kable(tmp, row.names = F, digits = 1, caption = "Code de Zonage en aire urbaine dans chaque groupe",
      col.names = c("Groupes", "Nombre de communes", "Unité urbaine", "Couronne périurbaine", "Communes multipolarisées", "Unité urbaine", "Couronne périurbaine", "Unité urbaine", "Couronne périurbaine", "Autres communes multipolarisées", "Communes hors influence des pôles"), format = "html") %>%
  kableExtra::add_header_above(c(" " = 2, "Grandes aires urbaines" = 3, "Pôles moyens" = 2, "Petits pôles" = 2," " = 2)) 

x <- length(NomsGroupesCAH)

table <- compare[, 1:(4+x)]
table[, 5:(4+x)] <- round (table[, 5:(4+x)], 2)

min <- min(table[, c(2, 3, 5:ncol(table))])
max <- max(table[, c(2, 3, 5:ncol(table))])

# Pour donner la description complète et propre de la variable
table <- merge(table, ratio[, 1:2], by.x = "Variable", by.y = "CODE")

# condformat(table) %>%
condformat(table[, c(length(table), 2:(length(table)-1))]) %>% # Pour avoir la description complète des variables
          # rule_fill_bar("Valeur pour la France entière (A)", limits = c(0, 100), low = "darkgreen", high = "darkgreen") %>%
          
          rule_fill_gradient2("France", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("CommunesFusionnantes", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 1 PPCQ", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 2 PFO", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 3 MOARO", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 4 VAR", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 5 REAV", low = "darkblue", high = "red", limits = c(min, max)) %>%
          
  theme_caption(caption = "Tableau de comparaison entre les groupes") %>%
  theme_htmlWidget(number_of_entries = nrow(compare))

table2 <- table[, c(length(table), 2:4)]
min <- min(table2[4])
max <- max(table2[4])
colnames(table2) <- c ("Variable", 
                            "Valeur pour la France entière (A)",
                            "Valeur pour les communes fusionnantes (B)",
                            "Différence ComFus/France (B-A)")
print(condformat(table2) %>%
  rule_fill_gradient2("Différence ComFus/France (B-A)", low = "blue", high = "red", limits = c(min, max)) %>%
  theme_caption(caption = "Situation des communes fusionnantes vis-à-vis des communes françaises") %>%
  theme_htmlWidget(number_of_entries = nrow(compare)))


test <- melt(compare[, c(1, 2, 5:(4+x)) ], id.vars = "Variable")
colnames(test) <- c("Variable", "Groupes", "Valeur")
ggplot(test) +
  geom_bar(aes(x = Variable, y = Valeur, fill = Groupes),
           stat = "identity") +
  #  scale_fill_grey() +
  scale_fill_manual(values=c("#a65628", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00")) + 
  #  scale_fill_manual(values=c("#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33","#a65628")) + # Pour 7
  facet_wrap(~ Groupes) +
  coord_flip() + theme_bw()


test <- melt(compare[, c(1, 4, 11, 13, 15, 17, 19) ], id.vars = "Variable")
colnames(test) <- c("Variable", "Différence avec la moyenne des communes fusionnantes", "Valeur")
levels(test$`Différence avec la moyenne des communes fusionnantes`) <- c("France",
                         "Groupe 1", "Groupe 2", "Groupe 3", "Groupe 4", "Groupe 5")
ggplot(test) +
  geom_bar(aes(x = Variable, y = Valeur, fill = `Différence avec la moyenne des communes fusionnantes`),
           stat = "identity") +
  #  scale_fill_grey() +
  scale_fill_manual(values=c("#a65628", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00")) + 
  facet_wrap(~ `Différence avec la moyenne des communes fusionnantes`) +
  coord_flip() + theme_bw() +
  labs(title = "Différence avec la moyenne des communes fusionnantes")+ theme(legend.position='none')

x <- length(NomsGroupesCAH)
min <- min(compare[, (5+x):ncol(compare)])
max <- max(compare[, (5+x):ncol(compare)])
rm(x)

# Boucle donnant les valeurs pour chaque groupe, avec éléments de comparaison d'un profil avec la moyenne française
for (i in NomsGroupesCAH){
  clusLongGr <- subset(clusLong, clusLong$Groupes == i)
  compare <- merge(moyennesCAH, clusLongGr[, c("variable", "value")], by.x = "Variable", by.y = "variable")
  compare$DiffAvecFr <- compare$value - compare$France
  compare$DiffAvecComFus <- compare$value - compare$CommunesFusionnantes
  compare[, 2:length(compare)] <- round(compare[, 2:length(compare)], 1)

  n <- stringr::str_sub(i, 8, 8)
  assign(paste0("compare_Gr", n), compare)
  # print(kable(compare, row.names = F, caption = i))
  colnames(compare) <- c ("Variable", 
                            "Valeur pour la France entière (A)",
                            "Valeur pour les communes fusionnantes (B)",
                            "Différence ComFus/France (B-A)",
                            "Valeur pour le groupe étudié (C)",
                            "Différence Groupe/France (C-A)",
                            "Différence Groupe/ComFus (C-B)")
#  min <- min(compare[, c("Différence ComFus/France (B-A)", "Différence Groupe/France (C-A)", "Différence Groupe/ComFus (C-B)")])
#  max <- max(compare[, c("Différence ComFus/France (B-A)", "Différence Groupe/France (C-A)", "Différence Groupe/ComFus (C-B)")])

  table <- condformat(compare) %>%
          rule_fill_bar("Valeur pour la France entière (A)", limits = c(0, 100), low = "darkgreen", high = "darkgreen") %>%
          rule_fill_bar("Valeur pour les communes fusionnantes (B)", limits = c(0, 100), low = "darkgreen", high = "darkgreen") %>%
          rule_fill_bar("Valeur pour le groupe étudié (C)", limits = c(0, 100), low = "darkgreen", high = "darkgreen") %>%
          rule_fill_gradient2("Différence ComFus/France (B-A)", low = "darkblue", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Différence Groupe/France (C-A)", low = "darkblue", high = "red",  limits = c(min, max)) %>%
          rule_fill_gradient2("Différence Groupe/ComFus (C-B)", low = "darkblue", high = "red",  limits = c(min, max))%>%
          theme_caption(caption = i) %>%
          theme_htmlWidget(number_of_entries = nrow(compare))
  print(knit_print(table))
  rm(compare, n, table)
} 

# Autre possibilité, qui m'a moins convaincue : https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html#celltext_specification
# library(formattable)
# library(kableExtra)
# library(dplyr)
# compare_Gr2 %>%
#   mutate (
#     DiffComFusComFr = color_tile("white", "orange") (DiffComFusComFr),
#     DiffAvecFr = color_tile("white", "orange") (DiffAvecFr),
#     DiffAvecComFus = color_tile("white", "orange") (DiffAvecComFus),
#     France = color_bar("lightgreen")(France),
#     CommunesFusionnantes = color_bar("lightgreen")(CommunesFusionnantes),
#     value = color_bar("lightgreen")(value)
#     )%>%
#   kable(escape = F, align = "r", digits = 1) %>%
#   kable_styling(c("striped", "condensed"), full_width = F)

# NB : On reprend les données déjà élaborées dans le 5.2.2
# Import données totales
moyennesComCAH <- data.frame()
for (i in selecVarCAH){
  a <- round(mean(df2011[, i], na.rm = TRUE), 2)
  b <- round(median(df2011[, i], na.rm = TRUE), 2)
  c <- round(mean(datafus2011[, i], na.rm = TRUE), 2)
  d <- round(median(datafus2011[, i], na.rm = TRUE), 2)
  
  h <- c(i, a, b, c, d)
  moyennesComCAH <- rbind(moyennesComCAH, h, stringsAsFactors= FALSE)
  rm( a, b, c, d, h)
}
# Pour calcul valeur moyenne retraités ("Inf" vient perturber) : 
# a <- subset(df2011$P09_RETR1564_RT, df2011$P09_RETR1564_RT != Inf)
# mean(a, na.rm =TRUE)
# rm(a)

colnames(moyennesComCAH) <- c("Variable", "Comfr_moy", "Comfr_med", "Comfus_moy", "Comfus_med")
moyennesComCAH <- merge(moyennesComCAH, ratio[, 1:2], by.x = "Variable", by.y = "CODE")
row.names(moyennesComCAH) <- moyennesComCAH$DESCRIPTION
moyennesComCAH[1] <- moyennesComCAH$DESCRIPTION
kable(moyennesComCAH[, 1:5], row.names = FALSE, col.names = c("", "Moyenne", "Médiane", "Moyenne", "Médiane"), digits = 2, format = "html")  %>%
  kableExtra::add_header_above(c("Variable" = 1, "Communes françaises" = 2, "Communes fusionnantes" = 2))

5.4 Chaque CN est-elle composée d’un seul type en fonction de la typologie réalisée ?

À partir de la typologie réalisée, permettant de catégoriser les communes fusionnantes sur la base de données socio-économiques, nous pouvons étudier la composition des communes nouvelles sur ce plan.

# Définition d'un tableau comportant le nombre de communes fusionnantes par commune nouvelle
count_CN_typo <- plyr::count(typo, "CODGEO_new")

mesCommunes <- count_CN_typo$CODGEO_new
# mesCommunes <- c("61324","73150","73006")
#pourinfo <- subset (typo, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")

df <- typo
results <- data.frame(matrix(ncol=nclass, nrow=0))

for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
  toto <- subset (df, CODGEO_new == i ) # Ne garder que les communes fusionnantes y ayant participé
  a <- table(toto$Groupes) # Relever les groupes des communes fusionnantes
  results <- rbind(results,a) # Combiner les résultats, par lignes
  rm(a, toto) # Supprimer "a"
}
# On renomme les colonnes
colnames(results) <- NomsGroupesCAH
count_CN_typo <- cbind(count_CN_typo, results)

# On identifie la variable la plus fréquente
count_CN_typo$max2 <- apply(count_CN_typo[, 3:(nclass+2)], 1, function(x) max(x, na.rm = TRUE))
# Quel est le groupe de typo revenant le plus fréquemment dans une CN
count_CN_typo$Typomaj <- colnames(count_CN_typo[, 3:(nclass+2)])[apply(count_CN_typo[, 3:(nclass+2)], 1, which.max)]

# On note si une CN a des communes fusionnantes avec un groupe de typo identique
count_CN_typo$TypoIdent <- ifelse(count_CN_typo$max == count_CN_typo$freq, TRUE, FALSE)
summary(count_CN_typo$TypoIdent)

##    Mode   FALSE    TRUE 
## logical     584     260

# On note si une CN a des communes fusionnantes avec un groupe de typo presque identique (max 25% de communes ayant un type différent)
count_CN_typo$TypopresqIdent <- ifelse((count_CN_typo$freq - count_CN_typo$max2)/count_CN_typo$freq <= 0.25, TRUE, FALSE)
summary(count_CN_typo$TypopresqIdent)

##    Mode   FALSE    TRUE 
## logical     535     309

# On extrait les CN ayant des communes fusionnante avec un groupe de typo identique
CNTypoIdent <- subset(count_CN_typo, TypoIdent==TRUE)
CNTypoIdent <- merge(CNTypoIdent, typo[, c("CODGEO", "Groupes")], by.x = "CODGEO_new", by.y = "CODGEO", all.x = TRUE)

  
# Quel pourcentage de chaque groupe dans les communes nouvelles homogènes ?
tabcont<-summary(CNTypoIdent$Groupes)
# Quel pourcentage de chaque groupe est majoritaire ?
count_CN_typo$Typomaj <- as.factor(count_CN_typo$Typomaj)
tabcont2<-summary(count_CN_typo$Typomaj)
# Pour comparaison : pourcentages de chaque groupe dans les communes fusionnantes
tabcont3<-summary(typo$Groupes)
round(100*prop.table(tabcont,margin=),1) # Pourcentages de chaque groupe en fonction de la totalité des communes nouvelles homogènes

##  Groupe 1 PPCQ   Groupe 2 PFO Groupe 3 MOARO   Groupe 4 VAR  Groupe 5 REAV 
##           33.8           20.8           28.8            2.3           13.1 
##           NA's 
##            1.2

round(100*prop.table(tabcont2,margin=),1) # Pourcentages de chaque groupe majoritaire en fonction de la totalité des communes fusionnantes

##  Groupe 1 PPCQ   Groupe 2 PFO Groupe 3 MOARO   Groupe 4 VAR  Groupe 5 REAV 
##           38.9           22.4           27.3            3.9            7.6

round(100*prop.table(tabcont3,margin=),1) # Pourcentages de chaque groupe en fonction de la totalité des communes fusionnantes

##  Groupe 1 PPCQ   Groupe 2 PFO Groupe 3 MOARO   Groupe 4 VAR  Groupe 5 REAV 
##           24.7           22.6           33.1            4.8           14.9

tableau <- rbind(
  summary(CNTypoIdent$Groupes),
  round(100*prop.table(tabcont,margin=),1) * 100 / nrow(typo),
  round(100*prop.table(tabcont,margin=),1),
  summary(count_CN_typo$Typomaj),
  round(100*prop.table(tabcont2,margin=),1),
  summary(typo$Groupes),
  round(100*prop.table(tabcont3,margin=),1)
)
row.names(tableau) <- c(
  "Nombre de communes nouvelles regroupant des communes homogènes du point de vue de la typologie",
  "Pourcentages en fonction de la totalité des communes nouvelles",
  "Pourcentages en fonction de la totalité des communes nouvelles homogènes",
  "Nombre de communes nouvelles ayant un type majoritaire",
  "Pourcentages de chaque groupe majoritaire en fonction de la totalité des communes fusionnantes",
  "Nombre de communes fusionnantes par groupe",
  "Pourcentages de chaque groupe en fonction de la totalité des communes fusionnantes")

kable(tableau, align = "c", digits = c(0, 1, 1, 0, 1, 0, 1))

Conclusion : Une large majorité des communes nouvelles sont, du point de vue de notre typologie, composées de communes au profil hétérogène. Lorsque les communes nouvelles sont homogène du point de vue du profil des communes fusionnantes, il s’agit, plus fréquemment, des groupes 1, 2 et 3. Mais si on compare ces pourcentages à la place de ces groupes dans la population totale, on observe surtout qu’un groupe est bien plus présent dans les communes homogènes que dans la population totale : le groupe 1. Inversement, le groupe 5 est bien moins présent dans les communes nouvelles homogènes que dans l’ensemble des communes fusionnantes. On pourrait donc mettre en avant la création de communes nouvelles plus fréquemment homogènes en contexte urbain.

5.5 Typologie et distance euclidienne

# Import des données concernant la typologie
# typo <- st_read("sorties/typo.gpkg", layer = "typo", quiet = TRUE)
# typo <- subset(typo, ChefLieu == "O")
# count_CN_typo <- plyr::count(typo, "CODGEO_new")

disteucl <- st_read("sorties/disteucl.gpkg", layer = "disteucl", quiet = TRUE)


disteucl <- merge (disteucl, count_CN_typo, by = "CODGEO_new", all.y = FALSE)


GroupesDisteucl <- by(disteucl$dist_CN, disteucl$Typomaj, mean, na.rm = TRUE)
GroupesDisteucl <- data.frame(as.table (GroupesDisteucl))
colnames(GroupesDisteucl) <- c("Typomaj", "mean")
GroupesDisteucl$sd <- as.numeric(by(disteucl$dist_CN, disteucl$Typomaj, sd, na.rm = TRUE))
GroupesDisteucl$median <- as.numeric(by(disteucl$dist_CN, disteucl$Typomaj, median, na.rm = TRUE))
# On rajoute les valeurs moyennes pour faciliter les comparaisons
levels(GroupesDisteucl$Typomaj) <- c(levels (GroupesDisteucl$Typomaj), "Ensemble")
GroupesDisteucl[ nrow(GroupesDisteucl)+1, ] <- c("Ensemble",
                                                 mean(disteucl$dist_CN, na.rm = TRUE),
                                                 sd(disteucl$dist_CN, na.rm = TRUE),
                                                 median(disteucl$dist_CN, na.rm = TRUE))
# Esthétique pour faciliter la lecture
GroupesDisteucl[, 2] <- as.numeric(GroupesDisteucl[, 2])
GroupesDisteucl[, 3] <- as.numeric(GroupesDisteucl[, 3])
GroupesDisteucl[, 4] <- as.numeric(GroupesDisteucl[, 4])
GroupesDisteucl[, 2:4] <- round(GroupesDisteucl[, 2:4], 1)

tmp <- GroupesDisteucl
colnames(tmp) <- c("Groupe majoritaire dans la commune nouvelle", "Moyenne", "Écart-type", "Médiane")
kable(tmp, row.names = F, digits = 1, caption = "Typologie et distance euclidienne intra-communes nouvelles")

condformat(tmp[, c(1,2, 4)]) %>%
          
          rule_fill_gradient2("Moyenne", low = "darkblue", high = "red", limits = c(min (tmp$Moyenne), max (tmp$Moyenne))) %>%
          rule_fill_gradient2("Médiane", low = "darkblue", high = "red", limits = c(min(tmp$Médiane), max(tmp$Médiane))) %>%
         
  theme_caption(caption = "Typologie et distance euclidienne intra-communes nouvelles") %>%
  theme_htmlWidget(number_of_entries = nrow(tmp))

L’observation de la distance euclidienne en fonction des catégories de communes précédemment élaborés permet de visualiser les différences entre les groupes. Les communes dont le type majoritaire est le 1 et le 2, par exemple, sont en moyenne plus homogènes que les communes nouvelles où les groupes 3, 5 et surtout 4 dominent.

5.6 Typologie vis-à-vis de quelques autres données

# Création d'une variable surface
# Si on part de rien, jouer les deux lignes ci-dessous
# typo <- st_read("sorties/typo.gpkg", quiet = TRUE)
comparaisons <- data.frame( table (typo$Groupes))

# Création d'une variable de surface
typo$surface <- st_area(typo) # Attention, unités : m²
typo$surface <- set_units(typo$surface, km^2) # On passe en km²
typo$densite <- typo$P09_POP/typo$surface

comparaisons$surface_moy <- round(tapply (typo$surface, typo$Groupes, mean), 2) # On fait la moyenne de la surface
comparaisons$surface_med <- round(tapply (typo$surface, typo$Groupes, median), 2) # On fait la médiane de la surface

comparaisons$pop_moy <- round(tapply (typo$P09_POP, typo$Groupes, mean), 2) # On fait la moyenne de la population
comparaisons$pop_med <- round(tapply (typo$P09_POP, typo$Groupes, median), 2) # On fait la médiane de la population

comparaisons$dens_moy <- round(tapply (typo$densite, typo$Groupes, mean), 2) # On fait la moyenne de la population
comparaisons$dens_med <- round(tapply (typo$densite, typo$Groupes, median), 2) # On fait la médiane de la population


tmp <- comparaisons[, c("Var1", "Freq", "surface_moy", "surface_med", "pop_moy", "pop_med", "dens_moy", "dens_med")]

colnames(tmp) <- c("Groupe", "Fréquence", "Surface moyenne (km²)", "Surface médiane (km²)", "Population moyenne", "Population médiane", "Densité moyenne", "Densité médiane")

table <- condformat(tmp)%>%
          rule_fill_gradient2("Surface moyenne (km²)", low = "blue", high = "red") %>%
          rule_fill_gradient2("Surface médiane (km²)", low = "blue", high = "red") %>%
          rule_fill_gradient2("Population moyenne", low = "blue", high = "red") %>%
          rule_fill_gradient2("Population médiane", low = "blue", high = "red") %>%
          rule_fill_gradient2("Densité moyenne", low = "blue", high = "red") %>%
          rule_fill_gradient2("Densité médiane", low = "blue", high = "red") %>%
  theme_caption(caption = "Groupes de la typologie vis-à-vis de la surface et de la population des communes fusionnantes") %>%
  theme_htmlWidget(number_of_entries = nrow(tmp))
print(table)

# On liste, pour chaque région, le nombre de communes appartenant à chaque groupe
typoregions <- data.frame(table(typo$REG, typo$Groupes))


# Import d'un fichier donnant le nom des régions en fonction de leur code : https://www.insee.fr/fr/information/2560625#titre-bloc-29
names <- data.frame(read_excel("data-raw/stats_insee/table-appartenance-geo-communes-11.xls", sheet = "Niv_supracom", skip = 5))
names <- subset(names , names$NIVGEO == "REG")
typoregions <- merge (typoregions, names[, c("CODGEO", "LIBGEO")], by.x = "Var1", by.y = "CODGEO", all.x = TRUE, all.y = FALSE)

# Pour avoir les groupes en colonnes, plus pertinent de le faire avant l'ajout des noms
typoregions1 <- dcast(typoregions, LIBGEO ~ Var2, value.var = "Freq")

# Pour avoir les régions en colonnes
typoregions2 <- dcast(typoregions, Var2 ~ LIBGEO, value.var = "Freq")


min <- min(typoregions1[, 2:ncol(typoregions1)])
max <- max(typoregions1[, 2:ncol(typoregions1)])

# Si le format choisi est avec les groupes en colonnes
table <- condformat(typoregions1) %>% # ici les couleurs sont définies par les limites min/max
          rule_fill_gradient2("Groupe 1 PPCQ", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 2 PFO", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 3 MOARO", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 4 VAR", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Groupe 5 REAV", low = "white", high = "red", limits = c(min, max))
print(table)

table2 <- condformat(typoregions1) %>% # ici les couleurs sont les min/max automatiques, par colonnes
          rule_fill_gradient2("Groupe 1 PPCQ", low = "blue", high = "red") %>%
          rule_fill_gradient2("Groupe 2 PFO", low = "blue", high = "red") %>%
          rule_fill_gradient2("Groupe 3 MOARO", low = "blue", high = "red") %>%
          rule_fill_gradient2("Groupe 4 VAR", low = "blue", high = "red") %>%
          rule_fill_gradient2("Groupe 5 REAV", low = "blue", high = "red")
print(table2)

# Si le format choisi est avec les régions en colonnes
table3 <- condformat(typoregions2) %>% # ici les couleurs sont les min/max automatiques, par colonnes
          rule_fill_gradient2("Alsace", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Aquitaine", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Auvergne", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Basse-Normandie", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Bourgogne", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Bretagne", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Centre", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Champagne-Ardenne", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Franche-Comté", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Haute-Normandie", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Île-de-France", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Languedoc-Roussillon", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Limousin", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Lorraine", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Midi-Pyrénées", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Nord-Pas-de-Calais", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Pays de la Loire", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Picardie", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Poitou-Charentes", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Provence-Alpes-Côte d'Azur", low = "white", high = "red", limits = c(min, max)) %>%
          rule_fill_gradient2("Rhône-Alpes", low = "white", high = "red", limits = c(min, max))
print(table3)

table4 <- condformat(typoregions2) %>% # ici les couleurs sont les min/max automatiques, par colonnes
          rule_fill_gradient2("Alsace", low = "white", high = "red") %>%
          rule_fill_gradient2("Aquitaine", low = "white", high = "red") %>%
          rule_fill_gradient2("Auvergne", low = "white", high = "red") %>%
          rule_fill_gradient2("Basse-Normandie", low = "white", high = "red") %>%
          rule_fill_gradient2("Bourgogne", low = "white", high = "red") %>%
          rule_fill_gradient2("Bretagne", low = "white", high = "red") %>%
          rule_fill_gradient2("Centre", low = "white", high = "red") %>%
          rule_fill_gradient2("Champagne-Ardenne", low = "white", high = "red") %>%
          rule_fill_gradient2("Franche-Comté", low = "white", high = "red") %>%
          rule_fill_gradient2("Haute-Normandie", low = "white", high = "red") %>%
          rule_fill_gradient2("Île-de-France", low = "white", high = "red") %>%
          rule_fill_gradient2("Languedoc-Roussillon", low = "white", high = "red") %>%
          rule_fill_gradient2("Limousin", low = "white", high = "red") %>%
          rule_fill_gradient2("Lorraine", low = "white", high = "red") %>%
          rule_fill_gradient2("Midi-Pyrénées", low = "white", high = "red") %>%
          rule_fill_gradient2("Nord-Pas-de-Calais", low = "white", high = "red") %>%
          rule_fill_gradient2("Pays de la Loire", low = "white", high = "red") %>%
          rule_fill_gradient2("Picardie", low = "white", high = "red") %>%
          rule_fill_gradient2("Poitou-Charentes", low = "white", high = "red") %>%
          rule_fill_gradient2("Provence-Alpes-Côte d'Azur", low = "white", high = "red") %>%
          rule_fill_gradient2("Rhône-Alpes", low = "white", high = "red")
print(table4)

# On liste, pour chaque département, le nombre de communes appartenant à chaque groupe
# méthode #
dep_typo <- unique(typo$CODE_DEPT)
typodep <- data.frame()
for (i in dep_typo) {
  numero <- i
  subset <- subset (typo, CODE_DEPT == i)
  b <- table(subset$Groupes)
 resultdep <- c(numero, b)
 typodep <- rbind (typodep, resultdep, stringsAsFactors = FALSE)
 rm(numero, subset, b, resultdep)
}

# méthode 2 (revient au môme)
typodep <- data.frame(table(typo$CODE_DEPT, typo$Groupes))
typodep <- dcast(typodep, Var1 ~ Var2)


# Cartographie
dep <- merge(dep, typodep, by.x = "CODE_DEPT", by.y = "Var1", all.x = TRUE)
dep[is.na(dep)] <- 0

# Manipulation des noms de groupes car sinon variables non reconnues
NomsGroupesCAH2 <- str_replace_all(NomsGroupesCAH, " ", "_")
colnames(dep)[4:8] <- NomsGroupesCAH2

par(mfrow=c(2,3)) 
i <- NomsGroupesCAH2[5]
for (i in NomsGroupesCAH2) {
plot(st_geometry(dep))
propSymbolsLayer(x = dep,
                 var = i, symbols ="circle",
                 col =  "#B00EF0",
                 legend.pos = "left",
                 legend.title.txt = "Nombre de communes fusionnantes", inches = 0.2, fixmax = 150)
}

6.0 Préparation des données

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) # Les communes selon la géographie en vigueur au 1er janvier 2011.
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE) # L'ensemble des communes au 01/01/2025, selon la géographie en vigueur.
# NB : dans geom_new, les communes nouvelles ont comme caractéristiques l'addition des stock des communes qui les composent.
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE) # Les  communes nouvelles, avec les géométries au 1er janvier 2025 et caractérisées par les données à la géométrie 2011 agrégées.
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE) # Départements

# Import des données
load("data/refdata.Rdata")


datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes
dataCN_new <- subset(df_new, COM_NOUV == "OUI") # Les  communes nouvelles, avec les géométries au 1er janvier 2025 et caractérisées par les données à la géométrie 2011 agrégées.
dataNfus2011 <- subset(df2011, COM_NOUV == "NON") # Les communes, à la géométrie 2011, qui n'ont pas participé à la création d'une commune nouvelle

# Appariement des données
geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geom_new <- merge(geom_new, df_new, by = "CODGEO_new")
geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")

# Définition de sous-ensembles
testEdC <- subset (geom2011, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")
testEdCSavoie <- subset (geom2011, CODGEO_new =="73150" | CODGEO_new == "73006")
testNormandie <- subset (geom2011, REG == "23" | REG == "25")
test49 <- subset (geom2011, CODE_DEPT == "49" )
testOuest <- subset (geom2011, REG == "23" | REG == "25"| REG == "53"| REG == "52") # Normandies, Bretagne, Pays-de-la-Loire

# Import des données de la typologie
typo <- st_read("sorties/typo.gpkg", quiet = TRUE)

6.1 Croisement de la typologie avec d’autres variables

Y <- typo$P11_POT_FIN / typo$P09_POP
Y [1:10]

##  [1]  587.2569  605.2329  975.9486  661.8765  796.7216  575.3696  702.9351
##  [8] 1028.8303 1018.0400  619.4632

# Pour variables qualitatives
summary(Y)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   296.8   490.6   592.0   740.3   765.4 48890.3

Y<-cut(Y,breaks=c(quantile(Y)))
levels(Y)<-c("Q1","Q2", "Q3", "Q4")
Y [1:10]

##  [1] Q2 Q3 Q4 Q3 Q4 Q2 Q3 Q4 Q4 Q3
## Levels: Q1 Q2 Q3 Q4

X<-typo$Groupes
X [1:10]

##  [1] "Groupe 1 PPCQ"  "Groupe 2 PFO"   "Groupe 3 MOARO" "Groupe 1 PPCQ" 
##  [5] "Groupe 3 MOARO" "Groupe 3 MOARO" "Groupe 3 MOARO" "Groupe 2 PFO"  
##  [9] "Groupe 2 PFO"   "Groupe 1 PPCQ"

tabcont<-table(X,Y)
tabcont # En valeur absolue

##                 Y
## X                 Q1  Q2  Q3  Q4
##   Groupe 1 PPCQ  166 173 158 162
##   Groupe 2 PFO   233 173 107  89
##   Groupe 3 MOARO 210 218 229 224
##   Groupe 4 VAR     4  23  47  53
##   Groupe 5 REAV   53  80 126 139

round(100*prop.table(tabcont,margin=1),1) # Pourcentages, le total se fait par lignes

##                 Y
## X                  Q1   Q2   Q3   Q4
##   Groupe 1 PPCQ  25.2 26.3 24.0 24.6
##   Groupe 2 PFO   38.7 28.7 17.8 14.8
##   Groupe 3 MOARO 23.8 24.7 26.0 25.4
##   Groupe 4 VAR    3.1 18.1 37.0 41.7
##   Groupe 5 REAV  13.3 20.1 31.7 34.9

round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population

##                 Y
## X                 Q1  Q2  Q3  Q4
##   Groupe 1 PPCQ  6.2 6.5 5.9 6.1
##   Groupe 2 PFO   8.7 6.5 4.0 3.3
##   Groupe 3 MOARO 7.9 8.2 8.6 8.4
##   Groupe 4 VAR   0.1 0.9 1.8 2.0
##   Groupe 5 REAV  2.0 3.0 4.7 5.2

round(100*prop.table(tabcont,margin=2),1) # Pourcentages, le total se fait par colonnes

##                 Y
## X                  Q1   Q2   Q3   Q4
##   Groupe 1 PPCQ  24.9 25.9 23.7 24.3
##   Groupe 2 PFO   35.0 25.9 16.0 13.3
##   Groupe 3 MOARO 31.5 32.7 34.3 33.6
##   Groupe 4 VAR    0.6  3.4  7.0  7.9
##   Groupe 5 REAV   8.0 12.0 18.9 20.8

6.2 Test du Chi-2

Y <- testOuest$FUSION
# Pour variables qualitatives
summary(Y)
# Y<-cut(Y,breaks=c(quantile(Y)))
# levels(Y)<-c("Q1","Q2", "Q3", "Q4")
X <- testOuest$C09_ACTOCC

head(testOuest)

summary(X)
tabcont<-table(X,Y)
tabcont # En valeur absolue
round(100*prop.table(tabcont,margin=1),1) # Pourcentages, le total se fait par lignes
round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population
round(100*prop.table(tabcont,margin=2),1) # Pourcentages, le total se fait par colonnes


test<-chisq.test(tabcont)
test$observed
round(test$expected,1)
round(test$residuals,2)
test

6.3 Test de Fisher

6.4 ANOVA

Cf. http://perso.ens-lyon.fr/lise.vaudor/grimoireStat/_book/modeliser-le-lien-entre-deux-variables.html#lanova-lien-entre-une-variable-categorielle-et-une-variable-quantitative (@vaudor2018) mais aussi @feuillet2019 qui conseille, page 49, de préférer l’ANOVA au test de Fischer si la variable quantitative a plus de deux modalités.

7.0 Préparation des données

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) # Les communes selon la géographie en vigueur au 1er janvier 2011.
load("data/refdata.Rdata")
rm(df_new)
df2011 <- df2011[, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "CATAEU2010", "FusPhas", "COM_NOUV", "FUSION", "P09_POP", "REG", "CODE_DEPT")]
geom2011 <- merge(geom2011, df2011, by = "CODGEO")

7.1 Données et méthodes utilisées pour l’étude des navettes

L’INSEE fournit de nombreux fichiers concernant les navettes, disponibles ici : https://www.insee.fr/fr/statistiques/2022117. L’ensemble des données ne sont présentes que dans le document BTT_FM_DTR_2009.txt (le fichier BTT_FM_DTR_2009.xls ne comprend que les flux supérieurs à 100 unités ou les données concernant le nombre de navetteurs pour une commune donnée). Dans ce fichier, chaque flux est défini par sa commune de départ, sa commune d’arrivée et son volume (nombre de navetteurs). Il s’agit des données 2009, à la géographie en vigueur au 1er janvier 2011. Les métadonnées sont précisées dans le fichier BTX_FM_DTR_2009.xls.

# Import de l'intégralité des flux
NavettesToutes <- read.csv2("data-raw/Flux_domicile_travail/BTT_FM_DTR_2009.txt", stringsAsFactors=TRUE)

Pour exploiter ces données, nous nous sommes appuyés sur le package flows, dont la démarche est présentée ici (https://cran.r-project.org/web/packages/flows/vignettes/flows.html) (@giraud2016) et là (https://journals.openedition.org/netcom/2134) (@beauguitte2015).

# On choisit les données qu'on veut étudier pour les CN enquêtées
# Si on veut utiliser les données de démonstration du package flows
# data(nav)
# Si on veut toutes les données
nav <- NavettesToutes
# Si on veut toutes les communes nouvelles enquêtées
nav <- subset(NavettesToutes,
              CODGEO == "61155" | DCLT == "61155" | CODGEO == "61324" | DCLT == "61324"
              | CODGEO == "61455" | DCLT == "61455" | CODGEO == "73006" | DCLT == "73006" 
              | CODGEO == "73038" | DCLT == "73038" | CODGEO == "73093" | DCLT == "73093"
              | CODGEO == "73126" | DCLT == "73126" | CODGEO == "73150" | DCLT == "73150"
              | CODGEO == "73169" | DCLT == "73169" | CODGEO == "73305" | DCLT == "73305") 
# Pour ne prendre que les communes enquêtées en Savoie
nav <- subset(NavettesToutes, CODGEO == "73006" | DCLT == "73006" 
              | CODGEO == "73038" | DCLT == "73038" | CODGEO == "73093" | DCLT == "73093"
              | CODGEO == "73126" | DCLT == "73126" | CODGEO == "73150" | DCLT == "73150"
              | CODGEO == "73169" | DCLT == "73169" | CODGEO == "73305" | DCLT == "73305")
# Pour ne prendre que les communes fusionnantes de Passais-Villages
nav <- subset(NavettesToutes,  CODGEO == "61155" | DCLT == "61155" | CODGEO == "61324" | DCLT == "61324"
              | CODGEO == "61455" | DCLT == "61455") 
# Pour ne prendre que Aime (entrées et sorties)
nav <- subset(NavettesToutes, CODGEO == "73006" | DCLT == "73006")

# Préparation des données pour qu'elles soient adaptées au package
nav$CODGEO <- as.character(nav$CODGEO)
nav$DCLT <- as.character(nav$DCLT)
nav$L_DCLT <- as.character(nav$L_DCLT)

myflows <- prepflows(mat = nav, i = "CODGEO", j = "DCLT", fij = "NBFLUX_C09_ACTOCC15P")
myflows[1:4,1:4]

##       04096 05046 10323 29135
## 04096     0     0     0     0
## 05046     0     0     0     0
## 10323     0     0     0     0
## 29135     0     0     0     0

diag(myflows) <- 0

7.2 Premiers aperçus des flux

str(nav)

## 'data.frame':    98 obs. of  5 variables:
##  $ CODGEO              : chr  "04096" "05046" "10323" "29135" ...
##  $ LIBGEO              : Factor w/ 34125 levels "Aast","Abainville",..: 12939 9283 24425 16952 19836 30191 6125 9991 24257 28325 ...
##  $ DCLT                : chr  "73006" "73006" "73006" "73006" ...
##  $ L_DCLT              : chr  "Aime" "Aime" "Aime" "Aime" ...
##  $ NBFLUX_C09_ACTOCC15P: num  4 4.05 3.18 4.09 4.72 ...

# Affichage de graphiques
statmat(mat = myflows, output = "all", verbose = FALSE)

# Affichage de statistiques
statmat(mat = myflows, output = "none", verbose = TRUE)

## matrix dimension: 77 X 77 
## nb. links: 97 
## density: 0.01657553 
## nb. of components (weak) 1 
## nb. of components (weak, size > 1) 1 
## sum of flows: 1602.16 
## min: 0.791505 
## Q1: 3.890673 
## median: 4.033352 
## Q3: 8.156913 
## max: 252.8937 
## mean: 16.51711 
## sd: 37.39685

7.3 Sélection de certains flux

Grâce au package flows, plusieurs fonctions permettent de sélectionner certains flux. Parmi ces fonctions, On peut utiliser : - nfirst: sélectionne les k premiers flux de toutes les origines; - xfirst: sélectionne tous les flux supérieurs à un seuil k; - xsumfirst: sélectionne autant de flux que nécessaire pour chaque origine afin que leur somme soit au moins égale à k.

# Sélection des premiers flux
flowSel1 <- firstflows(mat = myflows, method = "nfirst", ties.method = "first",
                       k = 1)
# Selection des flux > 500
flowSel1 <- firstflowsg(mat = myflows, method = "xfirst", k = 5)

# Sélection à partir de plusieurs conditions
# Sélection des flux qui représentent au moins 20% des départs de la commune
flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xfirst", k = 1)
# Applique le 2e critère de Nystuen et Dacey (1961) :
# On garde le flux fij si la somme des flux reçus par j est supérieure à la somme des flux reçus par i.
# Sinon, on garde le flux fji.
flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)

# On multiplie pour obtenir uniquement les flux souhaités (flowSel 1 et 2 sont binaires)
flowSel <- myflows * flowSel1 * flowSel2

# Si on souhaite prendre tous les flux
flowSel <- myflows 

# Comparaison de la matrice initiale et de la matrice avec les flux sélectionnés
compmat(mat1 = myflows, mat2 = myflows * flowSel, digits = 1)

##                mat1     mat2  absdiff reldiff
## nblinks        97.0     97.0      0.0     0.0
## sumflows     1602.2 160721.4 159119.2  9931.5
## connectcompx    1.0      1.0      0.0      NA
## min             0.8      0.6       NA      NA
## Q1              3.9     15.1       NA      NA
## median          4.0     16.3       NA      NA
## Q3              8.2     66.5       NA      NA
## max           252.9  63955.2       NA      NA
## mean           16.5   1656.9       NA      NA
## sd             37.4   8019.2       NA      NA

7.4 Calculs des départs/arrivées totaux et de quelques indicateurs

Sur les indicateurs calculés, cf. http://grasland.script.univ-paris-diderot.fr/agreg/module6/index.html

# 1e solution
# On réalise un dataframe à partir de la matrice des flux pour pouvoir ajouter des données
# Flux <- as.data.frame(myflows)
# Calcul de la somme de chaque ligne (pour avoir la somme des départs)
# Flux$DépartsTotaux <- apply(Flux, 1, sum)
#Flux <- rbind(Flux, apply(Flux, 2, sum))

# 2e solution, en créant des dataframes à part
inflows <- data.frame(id = colnames(flowSel), w = colSums(flowSel))
outflows <- data.frame(id = rownames(flowSel), w = apply(flowSel, 1, sum))

# Attention, ne prennent en compte que les flux entrant ou sortant des communes sélectionnées
FluxTotaux <- merge(inflows, outflows, by = "id", all.x = TRUE)
navNames <- unique(NavettesToutes[, c("DCLT", "L_DCLT")])
FluxTotaux <- merge(FluxTotaux, navNames, by.x = "id", by.y = "DCLT", all.x = TRUE)
colnames(FluxTotaux)[2] <- "Arrivées"
colnames(FluxTotaux)[3] <- "Départs"
# FluxTotaux[order(FluxTotaux$Arrivées, decreasing = TRUE),][1:10,] # Affichage des communes qui attirent le plus
# FluxTotaux[order(FluxTotaux$Départs.Sorties, decreasing = TRUE),][1:10,] # Affichage des communes qui ont le plus de sorties

# Attention, n'a de sens que pour les communes qui ont la totalité de leurs flux renseignés
FluxTotaux$Volume <- FluxTotaux$Arrivées + FluxTotaux$Départs
FluxTotaux$Solde <- FluxTotaux$Arrivées - FluxTotaux$Départs
FluxTotaux$Attractivité <- FluxTotaux$Solde / FluxTotaux$Volume

7.5 Introduction aux pôles dominants/intermédiaires/dominés

# Transformation des géométries en un spdf utilisable pour la cartographie
ShpCommunes2011 <- as(geom2011, "Spatial")

# Que Passais-Villages
ShpEdC <- subset(ShpCommunes2011, CODGEO_new == "61324")
ShpEdC <- subset(ShpCommunes2011, CODE_DEPT == "61" | CODE_DEPT == "53")

# Que communes nouvelles étudiées en Savoie
ShpEdC <- subset(ShpCommunes2011, CODGEO_new =="73150" | CODGEO_new == "73006")

# rm(ShpCommunes2011)

sp::plot(ShpEdC, col = "#cceae7")
plotMapDomFlows(mat = flowSel, spdf = ShpEdC, spdfid = "CODGEO", w = inflows, wid = "id",
                wvar = "w", wcex = 0.05, add = TRUE,
                legend.flows.pos = "topright",
                legend.flows.title = "Nb. of commuters")
title("Dominant Flows of Commuters")
mtext(text = "INSEE, 2011", side = 4, line = -1, adj = 0.01, cex = 0.8)

Le package flows propose de catégoriser les entités en pôles dominants/intermédiaires/dominés. Ils sont définis selon les règles suivantes : (1) le flux le plus important de i est émis en direction de j ; (2) la somme des flux reçus par j est supérieure à la somme des flux reçus par i. - i est dominé si son flux le plus important est émis vers j qui reçoit un volume de flux plus important que lui. - j est intermédiaire si son flux le plus important est émis vers k qui reçoit un volume de flux plus important que lui et que, par ailleurs, il domine i. - k est dominant si son flux le plus important est vers un nœud qui reçoit un volume de flux moins important que lui.

# Ce script donne la liste des unités appartenant à chacune des catégories.
# https://gist.github.com/rCarto/f06aabf2852e4a9b8ba51e93c17eb5f6

Dominated <- colnames(flowSel)[colSums(flowSel) == 0]
Dominant <- rownames(flowSel)[rowSums(flowSel) == 0]
Intermediary <- rownames(flowSel)[!rownames(flowSel) %in% c(Dominant, Dominated)]
head(Dominated)

## [1] "04096" "05046" "10323" "29135" "34172" "34301"

head(Intermediary)

## [1] "73003" "73006" "73011" "73015" "73032" "73038"

head(Dominant)

## [1] "ZZZZZ" "16015" "57221" "73034" "73057" "73132"

# Pour cartographie
# labelLayer(spdf = ShpEdC[ShpEdC@data$CODGEO %in% Dominated,], txt = "CODGEO", col= "red")
# labelLayer(spdf = ShpEdC[ShpEdC@data$CODGEO %in% Dominant,], txt = "CODGEO", col= "blue")
# labelLayer(spdf = ShpEdC[ShpEdC@data$CODGEO %in% Intermediary,], txt = "CODGEO", col= "green")

7.6 Catégorisation de toutes les communes fusionnantes (par itérations successives)

Nous souhaitons observer la composition des communes nouvelles du point de vue de cette catégorisation en pôles dominants/intermédiaires/dominés.

# Script permettant de lancer les analyses des flux pour différents paramètres
# Lancé sur le serveur R d'Humanum fin mars 2021
# Export des résultats globaux dans sorties/Navettes et des synthèses dans sortie/

# Import des données  =======
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
load("data/refdata.Rdata")
df2011 <- df2011[, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "COM_NOUV")]
datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")
# testEdC <- subset (geomfus2011, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")
# testEdCSavoie <- subset (geomfus2011, CODGEO_new =="73150" | CODGEO_new == "73006")
NavettesToutes <- read.csv2("data-raw/Flux_domicile_travail/BTT_FM_DTR_2009.txt", stringsAsFactors=TRUE)


# Définition d'un vecteur avec les identifiants des communes nouvelles
mesCommunes <- unique(datafus2011$CODGEO_new)
# mesCommunes <- c("61324","73150","73006")
# mesCommunes <- c("73150","73006")


df <- NavettesToutes
# Création d'un champ indiquant si la commune origine appartient à une commune nouvelle
df <- merge(df, datafus2011[, c("CODGEO", "CODGEO_new")], by = "CODGEO", all.x = TRUE)
# Création d'un champ indiquant si la commune destination appartient à une commune nouvelle
df <- merge(df, datafus2011[, c("CODGEO", "CODGEO_new")], by.x = "DCLT", by.y = "CODGEO", all.x = TRUE)
# Adaptation des données au package flows
df$CODGEO <- as.character(df$CODGEO)
df$DCLT <- as.character(df$DCLT)
df$L_DCLT <- as.character(df$L_DCLT)
rm(NavettesToutes, geomfus2011) # Supprimer les objets pour alléger R


# xfirst % =====
# Flux supérieurs au pourcentage défini par le seuil s
seuils <- c(5, 10, 20, 50)

for (s in seuils) { # Pour chaque seuil
  
  results <- data.frame()
  for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
    nav1 <- subset (df, CODGEO_new.x == i | CODGEO_new.y == i ) # Ne garder que navettes partant ou arrivant des communes fusionnantes y ayant participé
    listeCommunes <- c(unique(nav1$DCLT), unique(nav1$CODGEO)) # On prélève les identifiants des communes concernées
    listeCommunes <- unique(listeCommunes) # En supprimant les doubles comptes
    compil_Navettes <- data.frame(matrix(ncol=7, nrow=0))
    
    for (z in listeCommunes) { # Pour chaque commune qui apparaît en source ou destination d'une des communes fusionnantes
      nav2 <- subset (df, DCLT == z) # On intègre tous ses flux...
      compil_Navettes <- rbind(compil_Navettes, nav2) #... dans un tableau...
    }
    
    nav <- unique(compil_Navettes) #... dont on supprime les doubles comptes
    myflows <- prepflows(mat = nav, i = "CODGEO", j = "DCLT", fij = "NBFLUX_C09_ACTOCC15P")
    diag(myflows) <- 0
    flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xfirst", k = s) 
    flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)
    # On multiplie pour obtenir uniquement les flux souhaités (flowSel 1 et 2 sont binaires)
    flowSel <- myflows * flowSel1 * flowSel2
    
    
    Dominated <- colnames(flowSel)[colSums(flowSel) == 0]
    Dominant <- rownames(flowSel)[rowSums(flowSel) == 0]
    Intermediary <- rownames(flowSel)[!rownames(flowSel) %in% c(Dominant, Dominated)]
    
    resultspartiels <- subset(datafus2011[, c("CODGEO", "CODGEO_new")], CODGEO_new == i) # Pour toutes les communes fusionnantes appartenant à la commune nouvelle étudiée (i), on note si elle est dominante, dominée ou intermédiaire
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominated] <- "Dominated"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominant] <- "Dominant"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Intermediary] <- "Intermediary"
    results <- rbind(results,resultspartiels) # Combiner les résultats, par lignes
    
    rm(resultspartiels, nav, nav1, nav2, compil_Navettes, myflows, flowSel, flowSel1, flowSel2, listeCommunes, Dominant, Dominated, Intermediary)
    
  }
  write.table(results, paste0("sorties/Navettes/export_results_flows_xfirst", s, "pct_domflows1.txt"), sep="\t", row.names=FALSE) 
}

# xfirst =====
# Flux supérieurs au seuil s
seuils <- c(05, 10, 50, 100)

for (s in seuils) { # Pour chaque seuil
  
  results <- data.frame()
  for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
    nav1 <- subset (df, CODGEO_new.x == i | CODGEO_new.y == i ) # Ne garder que navettes partant ou arrivant des communes fusionnantes y ayant participé
    listeCommunes <- c(unique(nav1$DCLT), unique(nav1$CODGEO)) # On prélève les identifiants des communes concernées
    listeCommunes <- unique(listeCommunes) # En supprimant les doubles comptes
    compil_Navettes <- data.frame(matrix(ncol=7, nrow=0))
    
    for (z in listeCommunes) { # Pour chaque commune qui apparaît en source ou destination d'une des communes fusionnantes
      nav2 <- subset (df, DCLT == z) # On intègre tous ses flux...
      compil_Navettes <- rbind(compil_Navettes, nav2) #... dans un tableau...
    }
    
    nav <- unique(compil_Navettes) #... dont on supprime les doubles comptes
    myflows <- prepflows(mat = nav, i = "CODGEO", j = "DCLT", fij = "NBFLUX_C09_ACTOCC15P")
    diag(myflows) <- 0
    flowSel1 <- firstflowsg(mat = myflows, method = "xfirst", k = s) 
    flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)
    # On multiplie pour obtenir uniquement les flux souhaités (flowSel 1 et 2 sont binaires)
    flowSel <- myflows * flowSel1 * flowSel2
    
    Dominated <- colnames(flowSel)[colSums(flowSel) == 0]
    Dominant <- rownames(flowSel)[rowSums(flowSel) == 0]
    Intermediary <- rownames(flowSel)[!rownames(flowSel) %in% c(Dominant, Dominated)]
    
    resultspartiels <- subset(datafus2011[, c("CODGEO", "CODGEO_new")], CODGEO_new == i) # Pour toutes les communes fusionnantes appartenant à la commune nouvelle étudiée (i), on note si elle est dominante, dominée ou intermédiaire
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominated] <- "Dominated"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominant] <- "Dominant"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Intermediary] <- "Intermediary"
    results <- rbind(results,resultspartiels) # Combiner les résultats, par lignes
    
    rm(resultspartiels, nav, nav1, nav2, compil_Navettes, myflows, flowSel, flowSel1, flowSel2, listeCommunes, Dominant, Dominated, Intermediary)
    
  }
  write.table(results, paste0("sorties/Navettes/export_results_flows_xfirst", s, "_domflows1.txt"), sep="\t", row.names=FALSE) 
}



# nfirst =====
# Les "s" premieurs flux
seuils <- c(1, 2, 5, 10)

for (s in seuils) { # Pour chaque seuil
  
  results <- data.frame()
  for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
    nav1 <- subset (df, CODGEO_new.x == i | CODGEO_new.y == i ) # Ne garder que navettes partant ou arrivant des communes fusionnantes y ayant participé
    listeCommunes <- c(unique(nav1$DCLT), unique(nav1$CODGEO)) # On prélève les identifiants des communes concernées
    listeCommunes <- unique(listeCommunes) # En supprimant les doubles comptes
    compil_Navettes <- data.frame(matrix(ncol=7, nrow=0))
    
    for (z in listeCommunes) { # Pour chaque commune qui apparaît en source ou destination d'une des communes fusionnantes
      nav2 <- subset (df, DCLT == z) # On intègre tous ses flux...
      compil_Navettes <- rbind(compil_Navettes, nav2) #... dans un tableau...
    }
    
    nav <- unique(compil_Navettes) #... dont on supprime les doubles comptes
    myflows <- prepflows(mat = nav, i = "CODGEO", j = "DCLT", fij = "NBFLUX_C09_ACTOCC15P")
    diag(myflows) <- 0
    flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "nfirst", k = s)
    flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)
    # On multiplie pour obtenir uniquement les flux souhaités (flowSel 1 et 2 sont binaires)
    flowSel <- myflows * flowSel1 * flowSel2
    
    Dominated <- colnames(flowSel)[colSums(flowSel) == 0]
    Dominant <- rownames(flowSel)[rowSums(flowSel) == 0]
    Intermediary <- rownames(flowSel)[!rownames(flowSel) %in% c(Dominant, Dominated)]
    
    resultspartiels <- subset(datafus2011[, c("CODGEO", "CODGEO_new")], CODGEO_new == i) # Pour toutes les communes fusionnantes appartenant à la commune nouvelle étudiée (i), on note si elle est dominante, dominée ou intermédiaire
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominated] <- "Dominated"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominant] <- "Dominant"
    resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Intermediary] <- "Intermediary"
    results <- rbind(results,resultspartiels) # Combiner les résultats, par lignes
    
    rm(resultspartiels, nav, nav1, nav2, compil_Navettes, myflows, flowSel, flowSel1, flowSel2, listeCommunes, Dominant, Dominated, Intermediary)
    
  }
  write.table(results, paste0("sorties/Navettes/export_results_flows_nfirst", s, "_domflows1.txt"), sep="\t", row.names=FALSE) 
}

# Analyse des résultats =====
# Import pour lancer rapidement les tests
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
datafus2011 <- st_set_geometry(geomfus2011, NULL)

listetests <- c("xfirst5pct_domflows1", "xfirst10pct_domflows1", "xfirst20pct_domflows1", "xfirst50pct_domflows1", 
                "xfirst5_domflows1", "xfirst10_domflows1","xfirst50_domflows1","xfirst100_domflows1",
                "nfirst1_domflows1", "nfirst2_domflows1", "nfirst5_domflows1", "nfirst10_domflows1")

for (test in listetests) {
  
  ResultsFlux <- read.table(paste0("sorties/Navettes/export_results_flows_", test, ".txt"),
                            sep="\t", colClasses = "character", head = TRUE, stringsAsFactors = TRUE)
  
  mesCommunes <- unique(ResultsFlux$CODGEO_new)
  
  results <- data.frame() # On prépare le tableau avec les résultats
  
  for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
    Pranalyse <- subset(ResultsFlux, CODGEO_new == i) # On extrait les données concernant les communes fusionnantes d'une CN
    a <- sum(Pranalyse$StatutNavettes %in% "Dominant") # Nombre de pôles dominants dans la commune nouvelle
    b <- sum(Pranalyse$StatutNavettes %in% "Intermediary") # Nombre de pôles intermédiaires dans la commune nouvelle
    c <- sum(Pranalyse$StatutNavettes %in% "Dominated") # Nombre de pôles dominés dans la commune nouvelle
    d <- Pranalyse$StatutNavettes[Pranalyse$CODGEO == Pranalyse$CODGEO_new]
    
    resultspartiels <- c(i, a, b, c, d)
    results <- rbind(results,resultspartiels, stringsAsFactors= FALSE) # Combiner les résultats, par lignes
    
    rm(i, Pranalyse, a, b, c, d, resultspartiels)
    
  }
  colnames(results) <- c("CODGEO_new", "NbrDominant", "NbrIntermediary", "NbrDominated", "StatutChefLieu")
  
  
  write.table(results, paste0("sorties/Synthese_Dominants_", test, ".txt"), sep="\t", row.names=FALSE) 
}

# Import pour lancer rapidement les tests

geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
load("data/refdata.Rdata")

df2011 <- df2011[, c("LIBGEO_new", "CODGEO_new", "LIBGEO", "CODGEO", "ChefLieu", "COM_NOUV")]
datafus2011 <- subset(df2011, COM_NOUV == "OUI") # Désigne les données concernant les communes ayant participé à la création d'une commune nouvelle, appelées ici communes fusionnantes
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")
# testEdC <- subset (geom2011, CODGEO_new == "61324" | CODGEO_new =="73150" | CODGEO_new == "73006")
# testEdCSavoie <- subset (geom2011, CODGEO_new =="73150" | CODGEO_new == "73006")
# testNormandie <- subset (geom2011, REG == "23" | REG == "25")
NavettesToutes <- read.csv2("data-raw/Flux_domicile_travail/BTT_FM_DTR_2009.txt", stringsAsFactors=TRUE)


# Définition d'un vecteur avec les identifiants des communes nouvelles
mesCommunes <- unique(datafus2011$CODGEO_new)
# mesCommunes <- c("61324","73150","73006")
# mesCommunes <- unique(testEdCSavoie$CODGEO_new)
# mesCommunes <- c("73150","73006")
# mesCommunes <- unique(testNormandie$CODGEO_new)
# mesCommunes <- "73006"

# test <- datafus2011
# test$CODGEO <- as.numeric (test$CODGEO)
# test <- subset(test, CODGEO > "50400")
# mesCommunes <- unique(test$CODGEO_new)

df <- NavettesToutes
# Création d'un champ indiquant si la commune origine appartient à une commune nouvelle
df <- merge(df, datafus2011[, c("CODGEO", "CODGEO_new")], by = "CODGEO", all.x = TRUE)
# Création d'un champ indiquant si la commune destination appartient à une commune nouvelle
df <- merge(df, datafus2011[, c("CODGEO", "CODGEO_new")], by.x = "DCLT", by.y = "CODGEO", all.x = TRUE)
# Adaptation des données au package flows
df$CODGEO <- as.character(df$CODGEO)
df$DCLT <- as.character(df$DCLT)
df$L_DCLT <- as.character(df$L_DCLT)
rm(NavettesToutes, geomfus2011) # Supprimer les objets pour alléger R


results <- data.frame()

for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
  nav1 <- subset (df, CODGEO_new.x == i | CODGEO_new.y == i ) # Ne garder que navettes partant ou arrivant des communes fusionnantes y ayant participé
  listeCommunes <- c(unique(nav1$DCLT), unique(nav1$CODGEO)) # On prélève les identifiants des communes concernées
  listeCommunes <- unique(listeCommunes) # En supprimant les doubles comptes
  compil_Navettes <- data.frame(matrix(ncol=7, nrow=0))
  
  for (z in listeCommunes) { # Pour chaque commune qui apparaît en source ou destination d'une des communes fusionnantes
    nav2 <- subset (df, DCLT == z) # On intègre tous ses flux...
    compil_Navettes <- rbind(compil_Navettes, nav2) #... dans un tableau...
  }
  
  nav <- unique(compil_Navettes) #... dont on supprime les doubles comptes
  myflows <- prepflows(mat = nav, i = "CODGEO", j = "DCLT", fij = "NBFLUX_C09_ACTOCC15P")
  diag(myflows) <- 0
  # flowSel1 <- firstflowsg(mat = myflows, method = "xfirst", k = 5) # test 0
  # flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xfirst", k = 50) # test à faire 1
  flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xsumfirst", k = 60) # test à faire 2, avec autant de flux que nécessaires pour arriver à 60% des flux d'une commune
  # flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "nfirst", k = 2) # test à faire 3, avec, à chaque fois, les cinq premiers flux

  flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)
    # On multiplie pour obtenir uniquement les flux souhaités (flowSel 1 et 2 sont binaires)
  flowSel <- myflows * flowSel1 * flowSel2
  
  # flowSel <- myflows # Possibilité de sélectionner ?
  Dominated <- colnames(flowSel)[colSums(flowSel) == 0]
  Dominant <- rownames(flowSel)[rowSums(flowSel) == 0]
  Intermediary <- rownames(flowSel)[!rownames(flowSel) %in% c(Dominant, Dominated)]
  
  resultspartiels <- subset(datafus2011[, c("CODGEO", "CODGEO_new")], CODGEO_new == i) # Pour toutes les communes fusionnantes appartenant à la commune nouvelle étudiée (i), on note si elle est dominante, dominée ou intermédiaire
  resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominated] <- "Dominated"
  resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Dominant] <- "Dominant"
  resultspartiels$StatutNavettes[resultspartiels$CODGEO %in% Intermediary] <- "Intermediary"
  results <- rbind(results,resultspartiels) # Combiner les résultats, par lignes
  
  rm(resultspartiels, nav, nav1, nav2, compil_Navettes, myflows, flowSel, flowSel1, flowSel2, listeCommunes, Dominant, Dominated, Intermediary)
  
}

# write.table(results, "sorties/export_results_flows_xfirst05_domflows1.txt", sep="\t", row.names=FALSE) # test 0
# write.table(results, "sorties/export_results_flows_xfirst50pct_domflows1.txt", sep="\t", row.names=FALSE) # Fait
# write.table(results, "sorties/export_results_flows_xsumfirst60pct_domflows1.txt", sep="\t", row.names=FALSE) # test à faire 2
# write.table(results, "sorties/export_results_flows_nfirst5_domflows1.txt", sep="\t", row.names=FALSE) # fait sur serveur Humanum
# write.table(results, "sorties/export_results_flows_nfirst2_domflows1.txt", sep="\t", row.names=FALSE) # fait sur serveur Humanum

# Import pour lancer rapidement les tests
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) # Les  communes qui ont participé à la création de communes nouvelles (appelées communes fusionnantes).
datafus2011 <- st_set_geometry(geomfus2011, NULL)

test <- 
  # "xfirst05_domflows1" # Flux supérieurs à 5 navetteurs
  # "xfirst50pct_domflows1" # Flux représentant au moins 50% des navetteurs 
  # "xfirst20pct_domflows1" # Flux représentant au moins 20% des navetteurs
  # "xfirst10pct_domflows1" # Flux représentant au moins 10% des navetteurs
  # "xfirst05pct_domflows1" # Flux représentant au moins 05% des navetteurs
  # "nfirst2_domflows1" # Deux premiers flux
  # "nfirst5_domflows1" # Cinq premiers flux
  # "nfirst10_domflows1" # Dix premiers flux
  # "sumfirst60_domflows1" # Flux pour que leur addition représente au moins 60 % des flux.

ResultsFlux <- read.table(paste0("sorties/export_results_flows_", test, ".txt"), # Test 0
                           sep="\t", colClasses = "character", head = TRUE, stringsAsFactors = TRUE)

# Test 2 >>> Export absent
# ResultsFlux <- read.table("sorties/export_results_flows_xsumfirst60pct_domflows1.txt", sep="\t",
#                        colClasses = "character", head = TRUE, stringsAsFactors = TRUE)
# ResultsFlux <- read.table("sorties/sorties_20200710_Renaud/export_results_flows_xsumfirst60pct_domflows1.txt", sep="\t",
#                        colClasses = "character", head = TRUE, stringsAsFactors = TRUE)

# Seulement pour adapter les résultats sous anciennes dénominations
# summary(ResultsFlux$StatutNavettes)
# ResultsFlux$StatutNavettes <- as.factor(ResultsFlux$StatutNavettes)
# levels(ResultsFlux$StatutNavettes) <- c("Dominated", "Dominant", "Intermediary")
# colnames(ResultsFlux) <- c("CODGEO", "CODGEO_new", "StatutNavettes")



# Définition d'un vecteur comportant les identifiants des communes nouvelles à tester
# mesCommunes <- c("61324","73150","73006")
# mesCommunes <- c("73150","73006")
# mesCommunes <- c("73150")
# mesCommunes <- unique(testNormandie$CODGEO_new)
mesCommunes <- unique(ResultsFlux$CODGEO_new)

results <- data.frame() # On prépare le tableau avec les résultats

for (i in mesCommunes) { # Pour chaque identifiant de commune nouvelle, 
  Pranalyse <- subset(ResultsFlux, CODGEO_new == i) # On extrait les données concernant les communes fusionnantes d'une CN
  a <- sum(Pranalyse$StatutNavettes %in% "Dominant") # Nombre de pôles dominants dans la commune nouvelle
  b <- sum(Pranalyse$StatutNavettes %in% "Intermediary") # Nombre de pôles intermédiaires dans la commune nouvelle
  c <- sum(Pranalyse$StatutNavettes %in% "Dominated") # Nombre de pôles dominés dans la commune nouvelle
  d <- Pranalyse$StatutNavettes[Pranalyse$CODGEO == Pranalyse$CODGEO_new]

  resultspartiels <- c(i, a, b, c, d)
  results <- rbind(results,resultspartiels, stringsAsFactors= FALSE) # Combiner les résultats, par lignes
  
  rm(i, Pranalyse, a, b, c, d, resultspartiels)
  
}
colnames(results) <- c("CODGEO_new", "NbrDominant", "NbrIntermediary", "NbrDominated", "StatutChefLieu")


write.table(results, paste0("sorties/Synthese_Dominants_", test, ".txt"), sep="\t", row.names=FALSE)