2.3- Cartes et métadonnées
Comment trier l'information en fonction des codes INSEE, comment faire des mesures, convertir des longitudes/latitudes en distances
Le tout en langage R
1- Comptabiliser les points sur une aire, les points qui se superposent à un polygone.
library(raster)
# Contours administratifs
adm_fr <- getData('GADM', country='FRA', level=1)
plot(adm_fr)
# Simulons les points ! (rnorm sert à simuler des valeurs : pas d'inquiétude !)
x <- rnorm(400,2,2); y <- rnorm(400,45,2)
points(x,y,col="red",cex=0.5,pch=16)
# Comptabiliser les points répartis région par région
library("sp")
p<-SpatialPoints(data.frame(x,y))
# install.packages("GISTools")
library(GISTools)
poly.counts(p, adm_fr) -> res
setNames(res, adm_fr$NAME_1)
Résultats :
compil
Auvergne-Rhône-Alpes 18
Bourgogne-Franche-Comté 11
Bretagne 24
Centre-Val de Loire 12
Corse 10
Grand Est 14
Hauts-de-France 10
Île-de-France 12
Normandie 12
Nouvelle-Aquitaine 23
Occitanie 42
Pays de la Loire 7
Provence-Alpes-Côte d'Azur 21
Quelques lignes de code permettent de renvoyer le nombre de points par région et plus généralement le nombre d'objets superposés à un polygone.
2- Croiser des données, mesurer des distances... En construction
Établir une matrice des mesures des frontières communes entre chaque polygones.
La fonction boundary_matrix marche aussi bien avec les objets sp et sf.
library("raster")
adm <- getData('GADM', country='FRA', level=1)
# Mettre les longitudes-Latitudes en mètres (sinon ça marche pas)
crs.meters <- CRS("+init=epsg:28992") # testons aussi
# Remarque, CRS à vérifier, conversion incertaine...
adm2 <- spTransform(adm, crs.meters)
library("prioritizr")
boundary_matrix(adm2[adm2$NAME_1 %in% adm2$NAME_1[1:5],], str_tree = FALSE)-> matrice
colnames(matrice) <- adm$NAME_1[1:5]
rownames(matrice) <- adm$NAME_1[1:5]
print(round(matrice[,1:3]/1000,0))# Une partie de la matrice (en km)
Auvergne-Rhône-Alpes Bourgogne-Franche-Comté Bretagne
Auvergne-Rhône-Alpes 1624 492 .
Bourgogne-Franche-Comté 492 849 .
Bretagne . . 3671
Centre-Val de Loire 108 217 .
Corse . . .
On voit ainsi que la Corse ne partage pas de frontière avec la Bretagne (Oh, surprise !), en revanche, l'Auvergne&co partage 492 km de frontières avec la Bourgogne&co. Le périmètre de la Bretagne est de 3671 km.
Cette page est en construction, nous y adjoignions un code brute qui sera mis en forme prochainement.
Désolé pour le manque de clarté.
1) Comment lier deux tableaux (l'un avec des codes postaux de relevés), l'autre avec une correspondance INSEE/code postale
2) Comment mesurer une distance entre deux communes
3) Comment extraire les numéros des codes postaux pour coloriser les département en fonction d'une fréquence
4) En vrac :
RQGIS
https://jannes-m.github.io/RQGIS/index.html
Mettre des cartes en 3D - scipt 3D
###################################
# OUVRIR LES DONNEES
###################################
data <- read.csv2(file.choose(), encoding="UTF-8");head(data);dt <- data[-1,]
# Récupérer les codes postaux des parents et compter les individus par code
a <- table(dt$CP_parent,exclude=c("Code postal - adresse des parents","","999999","","Indiquez son code postal :_Notez \"99999\" si étranger."));code <- rownames(a);code
nb_code <- as.vector(a)
# Récupérer le code département correspondant à chaque code postal
code_dep <- floor(as.numeric(code)/1000)
###################################
# CHARGER UN FOND DE CARTE
###################################
#install.packages("raster");
library(raster)
adm_fr <- getData('GADM', country='FRA', level=2)
names(adm_fr); adm_fr$CCA_2 #summary(adm_fr)
# Récupérer les numéros de département de la France ou région Aquitaine
#raster_id_dep <- as.numeric(adm_fr$CCA_2[adm_fr$NAME_1=="Aquitaine"]);raster_id_dep
raster_id_dep <- as.numeric(adm_fr$CCA_2);raster_id_dep
raster_id_dep<- raster_id_dep[!is.na(raster_id_dep)]
###################################
# Pour chaque code département désiré dans rater_de_dep
# Faire le total d'individus
###################################
nb_dep <- c();id_unique =c()
for (i in unique(raster_id_dep)) {
nb_dep <- c(nb_dep , sum(nb_code[code_dep==i]))
id_unique<- c(id_unique,i)
}
#############################
# Céer un gradient de couleurs triées en fonction du vecteur nb_dep
#############################
colfunc<-colorRampPalette(c("white","red"))
colors <- (colfunc(100))
colors <- colors[rank(nb_dep)]
#############################
# TRACER LA CARTE
#############################
plot.new();par(fig=c(0.2,1,0,0.9), mar = c(6,4,0.5,2), new=TRUE)
plot(adm_fr,lwd=1)
#plot(adm_fr[adm_fr$NAME_1=="Aquitaine",],lwd=2);box()
for (i in c(1:length(raster_id_dep))) {
plot(adm_fr[adm_fr$CCA_2==id_unique[i],],add=T,col=colors[i])
}
par(fig=c(0,0.2,0,0.9), mar = c(2,0,0.5,2), new=TRUE)
# Tracer un graphique blanc à la place de la légende
plot(c(0,2),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = 'legend title',cex.main=0.5)
text(x=1.5, y = seq(0,1,l=5), labels = round(seq(0,max(nb_dep),l=5),0))
legend_image <- as.raster(matrix(rev(colfunc(20)), ncol=1))
rasterImage(legend_image, 0, 0, 1,1)
# POUR ALLER PLUS LOIN ET FAIRE DES LEGENDES AUTOMATIQUE AVEC DECOUPE PROPRE ET GRADIENT INDIQUE SUR LA CARTE VOIR LE LIEN SUIVANT
https://www.sylvaindurand.fr/cartographie-avec-R/
library(raster) ; library("sp") ; library("rgdal") ; library("maptools")
pointDistance(c(1, 44), c(1, 46), lonlat=TRUE)/1000
adm_fr <- getData('GADM', country='FRA', level=3)
plot(adm_fr)
getwd()
setwd("C:/Users/binet/Desktop/Train/Formation_R")
myfr <- readRDS("gadm36_FRA_5_sp.rds")
plot(myfr)
names(myfr)
myfr$CC_5
head(myfr)
# IMPORTANT
france@proj4string
install.packages("maptools") ; library("maptools")
france<- readOGR(dsn = ".", layer = "n_com_fla_000", p4s=CRS("+proj=longlat +ellps=WGS84")))
france<-readShapeSpatial("n_com_fla_000.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84"))
france2<-readOGR(dsn = ".", layer = "n_com_fla_000",p4s =CRS("+proj=longlat"))
x <- locator(n=1)
plot(france)2)
france2<-readShapeSpatial("n_com_fla_000.shp") # A remplacer par readOGR
plot(france2)
x <- locator(n=1)
install.packages("rgdal")
france1 <- spTransform(france, CRS("+proj=longlat +datum=WGS84"))
plot(france)
head(france)
names(france)
france$insee_comm
france$gid
plot(france[which(france$insee_comm==24035),])
summary(france[which(france$insee_comm==24035),])
coordinates(france)
coordinates(france[which(france$insee_comm==24035),])
france[which(france$insee_comm==24035),]@polygons
a <- france[which(france$insee_comm==24035),]@polygons[[1]]
summary(a)
x <- mean(coordinates(france[which(france$insee_comm==24035),]@polygons)[,1])
y <- mean(coordinates(dordogne@polygons[[1]]@Polygons[[1]])[,2])
a <- coordinates(france[which(france$insee_comm==24035),])
points(a[,1],a[,2],col="red")
# Exemple propre
# L'université Pau
a <- coordinates(france[which(france$insee_comm==64445),])
# L'étudiant de Bayonne
b <- coordinates(france[which(france$insee_comm==64102),])
pointDistance(as.vector(a), as.vector(b), lonlat=F)/1000
##############################
#
a <- coordinates(france1[which(france$insee_comm==64445),])
# L'étudiant de Bayonne
b <- coordinates(france1[which(france$insee_comm==64102),])
pointDistance(as.vector(a), as.vector(b), lonlat=F)/1000
plot(france)
x <- locator(n=1)
##################################
setwd("C:/Users/binet/Desktop/Train/Formation_R/NouvelleCarteGeo/")
france<-readShapeSpatial("n_com_fla_000.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84")) # A remplacer par readOGR()
plot(france)
x <- locator(n=1) ; x
************
setwd("C:/Users/binet/Desktop/Train/Formation_R/Carte_communes/")
france<-readShapeSpatial("communes-20181230.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84")) # A remplacer par readOGR
plot(france)
x <- locator(n=1) ; x
###########################
setwd("C:/Users/binet/Desktop/Train/Formation_R/Carte_communes/")
france<-readShapeSpatial("communes-20181230.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84")) # A remplacer par readOGR()
plot(france)
x <- locator(n=1) ; x
# L'université Pau
a <- coordinates(france[which(france$insee==64445),]);a
# L'étudiant de Bayonne
b <- coordinates(france[which(france$insee==64102),]);b
pointDistance(as.vector(a), as.vector(b), lonlat=T)/1000
#############################
# Ouvrir le tableau des étudiants
# Faire une boucle pour chaque valeurs de la colonne insee