library(knitr)
library(dplyr, warn.conflicts = FALSE)
options(dplyr.summarise.inform = FALSE)
library(sf, quietly=T)
library(mapsf)
library(RColorBrewer)4 Population, surface, densité
Le premier exemple de création d’indicateur concerne l’établissement de cartes de densité ce qui suppose de mobiliser à la fois les données individuelles (pour agréger la population) et les données géographiques (pour calculer la superficie).
4.1 Agrégation de la population
Une question intéressante concerne la marge d’erreur du calcul des populations selon que l’on part d’un échantillon plus ou moins dense (ici, à 1 pct ou 10 pct). Ce point méritera d’être discuté avec les spécialistes de statistique inférentielle du projet. Nous nous bornerons ici à calculer les indicateurs tirés de l’échantillon à 10%
On calcule pour chacune des unités territoriales de niveau I ou II la somme des individus pondérée par le poids des individus (PERWT) et on stocke le résultat dans un dossier commun à tous les recensments appelés all. Ceci évitera d’avoir à refaire le calcul dans les étapes suivantes.
## Chargement des données
rp <-readRDS(file = "ipums/rp/rp_fivecountries_samp10pct.RDS")
## Calcul des populations au niveau administratif I
pop1<-rp %>% group_by(GEOLEV1) %>%
summarise(POP=sum(PERWT)) %>%
mutate(GEOLEV1 = as.character(GEOLEV1))
saveRDS(pop1, "indic/all/POP_LEV1.RDS")
## Calcul des populations au niveau administratif II
pop2<-rp %>% group_by(GEOLEV2) %>%
summarise(POP=sum(PERWT)) %>%
mutate(GEOLEV2 = as.character(GEOLEV2))
saveRDS(pop2, "indic/all/POP_LEV2.RDS")4.2 Jointure avec les géométries
On va maintenant procéder à la jointure des données avec les fonds de carte puis stocker ces derniers enrichis par la population au format sf. On crée pour cela une procédure semi-automatique qu’on va applique au RP Bénin 2013 et qui sera ensuite répliquée sur les autres recensements étudiés.
4.2.1 Bénin, RP 2013
library(sf)
census="bj2013"
refmap="bj1979_2013"
# --------------- Niveau administratif I ----------------------
## Load population
pop1<-readRDS("indic/all/POP_LEV1.RDS")
## Read Geometry and create code
map1<-readRDS(paste0("ipums/geom/geo1_", refmap,".RDS")) %>% rename(GEOLEV1=GEOLEVEL1)
## Join
mapdon1<-left_join(map1,pop1) %>% select(code = GEOLEV1,
nom = ADMIN_NAME,
pays = CNTRY_NAME,
pop = POP,
geometry = geometry)%>%
filter(nom !="Unknown") %>%
st_make_valid()
## Save
saveRDS(mapdon1,paste0("indic/",census,"/","pop_lev1_",census,".RDS"))
# --------------- Niveau administratif II ---------------------- ## Load population
pop2<-readRDS("indic/all/POP_LEV2.RDS")
## Read Geometry and create code
map2<-readRDS(paste0("ipums/geom/geo2_", refmap,".RDS")) %>% rename(GEOLEV2=GEOLEVEL2)
## Join
mapdon2<-left_join(map2,pop2) %>% select(code = GEOLEV2,
nom = ADMIN_NAME,
pays = CNTRY_NAME,
pop = POP,
geometry = geometry)%>%
filter(nom !="Unknown") %>%
st_make_valid()
## Save
saveRDS(mapdon2,paste0("indic/",census,"/","pop_lev2_",census,".RDS"))
# Vérification
par(mfrow=c(1,2))
mf_map(mapdon1, type="base", col="lightyellow")
mf_map(mapdon1,type="prop", var="pop",inches = 0.1)
mf_layout(title = paste(census,"- Population au niveau I"),
credits = "Source : IPUMS, 2020 & INS Benin, 2013",
scale = FALSE, arrow=FALSE)
mf_map(mapdon2, type="base", col="lightyellow")
mf_map(mapdon2,type="prop", var="pop",inches = 0.1)
mf_layout(title = paste(census,"- Population au niveau II"),
credits = "Source : IPUMS, 2020 & INS Benin, 2013",
scale = FALSE, arrow=FALSE)
4.2.2 Burkina Faso, RP 2006
4.2.3 Mali, RP 2009
- N.B. : il a fallu supprimer un polygone vide concernant les valeurs “Unknown” et réparer la géométrie d’une unité en appliquant la fonction st_make_valid() du package sf.
4.2.4 Togo, RP 2010
- N.B. : il a fallu supprimer un polygone vide concernant les valeurs “Unknown” et réparer la géométrie d’une unité en appliquant la fonction st_make_valid() du package sf.
4.2.5 Sénégal, RP 2013
4.3 Calcul des superficies et densités
Maintenant que les variables de population sont incluses dans un fichier sf, on peut y ajouter la superficie à l’aide d’un simple calcul géométrique et en déduire la densité de population des unités territoriales. Le calcul des superficies est évidemment légèrement entaché d’erreur puisqu’il dépend de la précision du fonds de carte. Ainis, dans le cas Bénin, nous trouvons une superficie de 115768 km2 alors que la valeur donnée par la Banque Mondiale est de 114760 km2.
Les fichiers sont exportés au à la fois au format sf (en .RDS) et au format shapefile pour servir de point de départ à l’ensemble des analyses cartographiques ultérieures.
4.3.1 RP Bénin, 2013
library(sf)
census="bj2013"
# calcul des indicateurs au niveau I
map1<-readRDS(paste0("indic/",census,"/","pop_lev1_",census,".RDS"))
map1$sup<-as.numeric(st_area(map1)/(1000*1000))
map1$den<- map1$pop/map1$sup
saveRDS(map1, paste0("indic/",census,"/","map_lev1_",census,".RDS"))
st_write(map1, paste0("indic/",census,"/","map_lev1_",census,".shp"),delete_dsn = T, quiet=T)
# calcul des indicateurs au niveau II
map2<-readRDS(paste0("indic/",census,"/","pop_lev2_",census,".RDS"))
map2$sup<-as.numeric(st_area(map2)/(1000*1000))
map2$den<- map2$pop/map2$sup
saveRDS(map2, paste0("indic/",census,"/","map_lev2_",census,".RDS"))
st_write(map2, paste0("indic/",census,"/","map_lev2_",census,".shp"),delete_dsn = T, quiet=T)
# Vérification
mybreaks=c(0,5,10,20,40,80,160,320,10000)
mycols = brewer.pal(8,"YlOrBr")
par(mfrow=c(1,2))
mf_map(map1,type="choro",
var="den",
breaks=mybreaks,
pal=mycols,
leg_title = "hab./km2",
leg_val_rnd = 0)
mf_layout(title = paste(census,"- Densite de population au niveau I"),
credits = "Sources : IPUMS, 2020 & INS Bénin,2013",
scale = FALSE, arrow=FALSE)
mf_map(map2,type="choro",
var="den",
breaks=mybreaks,
pal=mycols,
leg_title = "hab./km2",
leg_val_rnd = 0,
lwd=0.3,
border = "gray80")
mf_map(map1,type="base",
add=T,
lwd=1,
col = NA,
border="black")
mf_layout(title = paste(census,"- Densite de population au niveau II"),
credits = "Sources : IPUMS, 2020 & INS Bénin,2013",
scale = FALSE, arrow=FALSE)
4.3.2 RP Burkina Faso, 2006
4.3.3 RP Mali, 2009
4.3.4 RP Sénégal, 2013
4.3.5 RP Togo, 2010
4.4 Cartogrammes démographiques
Nous procédons dans cette dernière étape à la création de cartogrammes où la surface des unité terrioriales est approximativement proportionelle à leur population à la date du recensement. L’opération de transformation est effectuée en dehors du logiciel R à l’aide du logiciel ScapeToad qui est l’un des plus performant. La transformation est réalisée sur le niveau administraif II puis celui-ci est agrégé dans R à l’aide du package sf pour reconstituer le niveau administratif I.
- N.B. : le calcul de l’anamorphose avec Scapetoad a été assez long et coûteux en calcul pour le Sénégal en raison de la complexité du tracé des limites. Un fonds généralisé est en général plus rapide à transformer.
4.4.1 RP Bénin, 2013
census = "bj2013"
# Load level2 and remove columns created by scapetoad
cartog2<-st_read(paste0("indic/",census,"/","cartog_lev2_",census,".shp"), quiet = T) %>%
select(code, nom, pays, pop, geometry)
st_write(cartog2,paste0("indic/",census,"/","cartog_lev2_",census,".shp"),delete_dsn = T , quiet = T)
# agregate level2 to level 1
cartog1 <- cartog2 %>% mutate(code2 = code, code=substr(code2,1,6))%>%
select(code, geometry) %>%
group_by(code) %>%
summarize() %>% st_as_sf() %>%
st_buffer( dist=0.001)
# load attribute from level1
don1<-st_read(paste0("indic/",census,"/","map_lev1_",census,".shp"), quiet = T) %>% st_drop_geometry()
cartog1 <- don1 %>% right_join(cartog1) %>% st_as_sf() %>% select(code, nom, pays, pop, geometry)
st_write(cartog1,paste0("indic/",census,"/","cartog_lev1_",census,".shp"),delete_dsn = T , quiet = T)
# Vérification
par(mfrow=c(1,2))
mf_map(cartog1, type="base", col="lightyellow")
mf_map(cartog1,type="prop", var="pop",inches = 0.1)
mf_layout(title = paste(census,"- Cartogramme au niveau I"),
credits = "Source : IPUMS, 2020 & INS Benin, 2013",
scale = FALSE, arrow=FALSE)
mf_map(cartog2, type="base", col="lightyellow")
mf_map(cartog2,type="prop", var="pop",inches = 0.1)
mf_layout(title = paste(census,"- Cartogramme au niveau II"),
credits = "Source : IPUMS, 2020 & INS Benin, 2013",
scale = FALSE, arrow=FALSE)
4.4.2 RP Burkina Faso, 2006
4.4.3 RP Mali, 2009
4.4.4 RP Sénégal, 2013
4.4.5 RP Togo, 2010
