Que faire des données spatiales ? C'est à cette question que nous tenterons de répondre dans cette étude. En effet, les données spatiales ou géographiques gagnent à être représentées sur une carte pour permettre de mieux en saisir le sens.
Les données travaillées dans cette études sont issues du site data.tours-metropole. Elles sont contenues dans la base de données Sirene où sont recensées les informations sur des millions d'entreprises et établissements actifs ou non, de tous les secteurs et dans toute la France. Dans cette étude, on se restreindra aux informations sur les entreprises et établissements de l'agglomération de Tours.
Nous consacrerons ce notebook à la représentation spatiale à l'aide des packages ggmap et leaflet à partir des données géographiques contenues dans cette table de données. En particulier, nous essayerons de représenter la densité d'entreprise par zone géographique. Que pourra-t-on conclure de ces représentations cartographiques ? Où les établissements et entreprises sont-ils le plus nombreux ? Existe-t-il des disparités en termes de dynamisme économique au sein de l'agglomération tourangelle ?
La table de données contient 127089 entrées et 108 variables. Nous avons donc beaucoup de variables et accessoirement beaucoup de données manquantes. La plupart de nos variables sont qualitatives. On va donc commencer par sélectionner les variables sur lesquelles on va travailler et préparer notre base de données pour la suite de l'étude.
Après étude approfondie des variables disponibles, on va sélectionner 18 variables parmis les 108 disponibles.
Nous les recodons pour plus de maniabilité.
Les variables conservées sont les suivantes :
SIRET : Le numéro Siret afin de pouvoir se débarasser des doublonsDATE_CREATION_ETAB : La date de création de l'établissementEFFECTIF_ETAB : La tranche d'effectif de l'établissementCP_ETAB : Le code postal de l'adresse de l'établissementVILLE_ETAB : La ville où est localisé l'établissementETAT_ADMIN_ETAB : L'état administratif de l'établissement (actif ou cessé)EMPLOYEUR_ETAB : Le caractère employeur ou non de l'établissementALTITUDE : L'altitude moyenne de la ville de l'établissementSUPERFICIE : La superficie de la ville de l'établissementPOP : La population de la ville de l'établissementEFFECTIF_UL : L'effectif de l'unité légaleCAT_JUR_UL : La catégorie juridique de l'unité légaleEMPLOYEUR_UL : Le caractère employeur de l'unité légaleSECTION_UL : La section de l'activité de l'unité légaleSSECTION_UL : La sous-section de l'activité de l'unité légaleNAT_JUR_UL : La nature juridique de l'unité légaleGEOLOCALISATION_ETAB : Les coordonnées GPS de l'établissementENSEIGNE_ETAB : L'enseigne de l'établissementOn change le type de certaines variables qui ne sont pas au bon format et on observe la structure de nos données.
Pour terminer, nous allons recoder la variable VILLE_ETAB car nous avons des tirets dans certains cas pour le nom des villes. Nous allons donc remplacer tous les tirets par des espaces dans les modalités de cette variable.
list<-str_replace_all(data$VILLE_ETAB, "-", " ")
data$VILLE_ETAB<-c(list)Notre table de données est maintenant prête à être exploitée.
L'établissement désigne l'entreprise ou l'organisme considéré. Nous avons 29367 établissements employant des salariés et 97716 n'en employant pas. Parmis les établissements, on peut observer les effectifs de salariés dans le graphique suivant.
Pour une grande partie des établissements, l'effectif n'est pas renseigné dans la table de données. Globalement, il y a de très peu d'établissements qui emploient beaucoup de salariés.
Il y a 47722 établissements encore actifs et 79367 établissements dont l'activité a pris fin.
Les informations sur l'unité légale nous permettent de connaître le champ d'activité des établissements.
Dans chaque section, il existe des sous-sections. Par exemple, pour la section "Commerce", nous avons quatre sous-sections.
Nous possédons également des informations sur la nature juridique des unités légales.
La table de données contient également la population, la superficie et l'altitude de la ville où est situé l'établissement.
Notre analyse descriptive des données étant réalisée, nous allons à présent nous concentrer sur les données géographiques contenues dans notre table de données.
Notre base de données contient des informations géographiques sur l'emplacement des établissements. Grâce à ces informations, nous allons pouvoir placer les entreprises sur une carte.
Pour plus de commodité, on commence par séparer la colonne GEOLOCALISATION_ETAB en deux variables :
gps1 : la latitudegps2 : la longitudedata<-data %>%
separate(col = "GEOLOCALISATION_ETAB",
into = paste0("gps", 1:2), sep = ",",
extra = "merge")La ville de Tours est caractérisée par les coordonnées GPS suivantes :
On s'attend donc à trouver des coordonnées de cet ordre là dans nos données.
Le graphique précédent nous confirme que les données géographiques dont nous disposons sont correctes.
Les observations des coordonnées GPS des différentes communes de l'aglomération Tourangelle sont cohérentes.
Avant toute autre chose, il faut une clé API Google pour récupérer les cartes depuis Google Maps. Pour cela, il faut visiter le site suivant : Obtenir une clé API Google. On s'enregistre sur le site et on obtient donc une chaîne de caractères qui constitue notre clé API. Il en faut une par ordinateur.
Pour plus d'informations, on peut consulter l'aide dans l'onglet "help" de RStudio avec les mots clé suivants : "register_google". On utilise ensuite la commande suivante pour s'enregistrer (notons que la spécification write=TRUE permet de rester enregistré).
register_google(key = "[XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX]", write = TRUE)Nous avons donc des données géographiques exploitables. On va les représenter sur une carte de l'agglomération de Tours. Nous allons dans cette partie utiliser le package ggmap. Le principe est de récupérer une carte à partir de coordonnées que l'on fait passer dans la fonction get_map et en l'occurrence get_stamenmap pour récupérer des cartes de la source "stamen".
Nous cherchons à afficher une carte. On va l'importer depuis google grâce à la commande suivante. On va spécifier les limites de notre carte grâce aux arguments left, right,... et à l'aide des graphiques précédents. On doit aussi spécifier le zoom (10 pour l'échelle ville) et le type de carte.
La commande ggmap permet d'afficher la carte que l'on vient d'importer.
map<-get_stamenmap(c(left = 0.45, bottom = 47.25, right = 0.9, top = 47.5), zoom=10, source = "stamen", maptype ="terrain")
ggmap(map)Nous obtenons donc une carte donc nous avons pu fixer les limites. Comme ggmap fonctionne avec ggplot2, on va pouvoir utiliser ses commandes pour ajouter des objets sur nos fonds de carte.
Les informations géographiques dont nous disposons nous permettent de placer des localisations précises sur une carte statique. Cependant, à la vue du nombre de données dont nous disposons (126903), placer les entreprises sur une carte n'est pas la meilleure idée que nous puissions avoir. Pour mieux comprendre nos données, étudier la densité d'entreprises par secteur semble plus pertinent.
Dans ce graphique, nous allons observer l'intensité des activités commerciales. On divise cette section en 4 sous-sections avec l'option facet_wrap. La transparence du nuage de points facilite la lisibilité. On utilise finalement geom_density_2d pour mieux se rendre compte de la concentration des entreprises.
data_commerce <- data %>%
filter(SECTION_UL == "Commerce")
ggmap(map3)+ geom_point(data=data_commerce, aes(x=gps2, y=gps1), alpha=0.3, colour="#74C476")+
geom_density2d(data = data_commerce, aes(x = gps2, y = gps1), size = 0.4, colour="#006D2C")+
facet_wrap(~SSECTION_UL)+
theme(legend.position="right", axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0, 0, -1, -1), 'lines'),
plot.title = element_text(color="#74C476", size=8))+
labs(y = "",
x = "",
title = "Intensité des activités commerciales")+ scale_fill_gradient()+ stat_ellipse(data=data_commerce, aes(x=gps2, y=gps1), colour="#006D2C")L'ajout d'une représentation de la dispersion ainsi que l'utilisation de facet_wrap nous permet donc d'y voir plus clair qu'avec un simple nuage de points.
Dans cette carte, nous allons montrer l'intensité du nombre d'industries manufacturières au sein de l'agglomération Tourangelle par petite zone géographique. On utilise cette fois qmplot pour gérer le zoom.
On utilisera geom_bin2d pour représenter le nombre d'industries par zone avec une coloration par densité.
data_IM<-data %>%
filter(
SECTION_UL %in% c("Industrie manufacturiere"))
ggmap(map8)+
geom_bin2d(data=data_IM, aes(x=gps2, y=gps1), alpha=0.9) +
labs(y = "Latitude",
x = "Longitude", title="Industrie manufacturière")+
theme(legend.position="right", axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0, 0, -1, -1),'lines'),
plot.title = element_text(color="#74C476", size=8))+
scale_fill_gradient(low="#74C476", high="#006D2C")Pour cette carte, nous allons tenter de représenter graphiquement la concentration des activités immobilières. On utilise la commande stat_density_2d pour avoir une meilleure représentation.
data_AI<-data %>%
filter(
SECTION_UL %in% c("Activites immobilieres"))
ggmap(map12)+geom_point(data=data_AI, aes(x=gps2, y=gps1), colour="#74C476", alpha=0.1)+
labs(y = "Latitude",
x = "Longitude", title="Activités immobilières")+
stat_density_2d(data=data_AI, aes(fill = ..level.., x=gps2, y=gps1), geom = "polygon", alpha = .6, color = NA) +
scale_fill_gradient2("count", low ="#74C476", mid = "#006D2C", high = "black", midpoint = 600)+
theme(legend.position="right", axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0, 0, -1, -1),'lines'),
plot.title = element_text(color="#74C476", size=8))Dans cet exemple, nous allons voir comment représenter la densité lorsque l'on a plusieurs catégories à représenter. On va créer un graphique de densité pour la latitude et un pour la longitude. On va les combiner avec un nuage de points sur une carte.
data_ED<-data %>%
filter(
SECTION_UL %in% c("Administration publique"))
a<-ggmap(map8)+geom_point(aes(x = gps2, y = gps1, colour = ETAT_ADMIN_ETAB), data = data_ED, size = 0.6) + stat_ellipse()+
theme(legend.position=c(0.25, 0.9), axis.line = element_blank(),
plot.margin = unit(c(0, 0, -1, -1), 'lines'),
plot.title = element_text(color="#74C476", size=8), legend.background = element_rect(fill="white", size=0.1, linetype="solid", colour="#74c493"))+
labs(y = "latitude",
x = "longitude",
title = element_blank())+scale_colour_manual(values=c("#74C476","#006D2C"))
xdensity <- ggplot(data_ED, aes(gps2, fill=ETAT_ADMIN_ETAB)) +
geom_density(alpha=.5) +
scale_fill_manual(values = c("#74C476","#006D2C")) +
labs(y = "",
x = "longitude", title="Administrations publiques")+
theme(legend.position="none", axis.line = element_blank(), axis.text.x=element_text(size=5), axis.text.y=element_text(size=5),
plot.margin = unit(c(0, 0, 0, 0), 'lines'),
plot.title = element_text(color="#74C476", size=8))
ydensity <- ggplot(data_ED, aes(gps1, fill=ETAT_ADMIN_ETAB)) +
geom_density(alpha=.5) +
scale_fill_manual(values = c("#74C476","#006D2C")) + coord_flip()+
labs(y = "",
x = "latitude")+
theme(legend.position="none", axis.line = element_blank(), axis.text.x=element_text(size=5), axis.text.y=element_text(size=5),
plot.margin = unit(c(0, 0, 0, 0), 'lines'))
blankPlot <- ggplot()+geom_blank(aes(1,1))+
theme(
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.line = element_blank()
)
grid.arrange(xdensity, blankPlot, a, ydensity,
ncol=2, nrow=2, widths=c(4, 1.4), heights=c(1.4, 4))Dans cet exemple, nous utilisons geom_point avec une coloration des données par état administratif. On ajoute également des représentations de la densité pour comprendre mieux où se situent les administrations publiques fermées.
Le package ggmap a des avantages et inconvénients. D'une part, comme il fonctionne avec ggplot2, on peut utiliser toutes les commandes de ce package pour customisez nos cartes. D'autre part, il n'est pas aisé de gérer le zoom et l'importation des fonds de cartes nécéssite une clé API. Seule la source "Stamen" fonctionne actuellement, ce qui restreint le nombre de fonds de carte disponibles. Nous verrons dans la suite une alternative à ggmap.
Le package leaflet a l'avantage de ne pas souffrir des mêmes inconvénients que ggmap. Sa dimmension interactive permet de choisir le zoom à sa convenance. De fait, comme un zoom adapté est possible, on peut représenter des nuages de points assez denses sur une carte et permettre à l'utilisateur de la carte la meilleure lisibilité possible.
Le package leaflet fonctionne sur le même principe que ggplot2 avec la superposition de diverses couches d'informations. Par contre, il a ses commandes propres.
addProviderTiles, addTiles : Choisir la couche de fondaddMarkers, addAwesomeMarkers, addCircleMarkers : Ajouter des marqueurs personnalisésaddPopups : Ajouter une bulle avec de l’informationaddPolygons, addCircles, addRectangles : Ajouter des surfaces en forme de polygones, de cercles ou rectanglesaddLegend : Ajouter une légendeaddMeasure : Ajouter un bouton qui permet de mesurer des distances et des surfacesaddEasyButton : Ajouter un bouton avec un niveau de zoom prédéfiniLes fonds de cartes disponibles sont très nombreux et consultables ici.
Dans cette carte, on va représenter l'emplacement des établissements de l'immobilier avec une coloration selon l'état administratif. On choisit le fond de carte Esri.WorldImagery.
library(leaflet)
cof <- colorFactor(palette = c("#74C476","#006D2C"), domain = NULL)
leaflet(data_AI) %>%
addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircleMarkers(~gps2, ~gps1,
popup=data_AI$SECTION_UL,
#label = ~htmlEscape(df_categ()$EFFECTIF_ETAB),
weight = 4,radius = 5,
color=~cof(data_AI$ETAT_ADMIN_ETAB),
stroke = F, fillOpacity = 0.7) %>%
addLegend(position = "bottomright", values = ~ETAT_ADMIN_ETAB, pal = cof, title = "Etat")Dans cette carte, on représente les effectifs d'entreprises et établissements par ville. Le rayon du cercle est calculé à partir de cet effectif. Plus précisément, le rayon du cercle est l'effectif divisé par 500 pour garder de la visibilité. On ajoute aussi un label avec le nom de la ville et le nombre d'établissements.
data_ville<-cbind(data[,c("VILLE_ETAB", "gps1", "gps2")])
data_ville=data_ville %>% distinct(VILLE_ETAB, .keep_all = TRUE)
df9<-data.frame(table(data$VILLE_ETAB))
names(df9)<-c("VILLE_ETAB", "FREQ")
dt<-inner_join(data_ville, df9, by=c("VILLE_ETAB"))leaflet(dt) %>% addTiles() %>% addProviderTiles("CartoDB.Voyager") %>%
addCircleMarkers(~gps2, ~gps1, label = ~ paste(VILLE_ETAB, "," , FREQ, " établissements"),
radius = ~FREQ/500,
color = "#74C476",
stroke = FALSE, fillOpacity = 0.5
)