Aujourd’hui, j’ai dû déterminer quels étaient les points géographiques (des stations météorologiques) se situant dans un cercle de centre et de rayon donnés. L’année dernière, j’avais tenté quelque chose, mais le résultat ne me plaisait pas vraiment. Je suis tombé sur une fonction qui m’avait échappé : gcDestination, du pakage maptools. Elle permet de retourner les coordonnées d’un point selon une distance et un angle donné.

Trois packages sont nécessaire pour pouvoir exécuter le code qui suit : 

library(maps)
library(maptools)
library(rgeos)

Pour illustrer le propos, penchons-nous sur la Bretagne.

Donnons-nous 10 points en Bretagne :

longitudes <- c(-2.373247, -2.523580, -3.112924, -1.474192, -4.151028, -4.367183, -2.143235, -1.528023, -3.733716, -3.311343)
latitudes <- c(47.69491, 48.30560, 48.39880, 48.24052, 48.03038, 48.53544, 48.49801, 48.01217, 47.87922, 48.75708)
lesPoints <- data.frame(lon=longitudes,lat=latitudes)

À présent, récupérons la carte de la France. J’utilise les méthodes que j’ai apprises sur le tutoriel de Baptiste Coulmont (section 2) et sur le billet d’Arthur Charpentier.

# La carte de la France
france <- map('france', plot=FALSE,fill=TRUE)
detpartement_bzh &lt;- france$names[c(25,30,32,41)] #bretagne
# La carte de la Bretagne
bretagne <- map('france',regions=detpartement_bzh, fill=TRUE, col="transparent", plot=FALSE)
# Les noms des departements
# c.f. http://stackoverflow.com/questions/13316185/r-convert-zipcode-or-lat-long-to-county
IDs <-sapply(strsplit(bretagne$names,":"),function(x) x[1])
# Convertir en objet SP
bretagne.sp <- map2SpatialPolygons(bretagne,IDs=IDs,proj4string=CRS("+proj=longlat +datum=wgs84"))

On va repérer des points particuliers pour chaque département. C’est l’occasion d’introduire la fonction gCentroid, du package rgeos, qui permet de trouver le barycentre d’un polygone.

barycentre <- gCentroid(bretagne.sp,byid=TRUE)

Le rendu visuel est le suivant :

Barycentre Bretagne

La fonction qui vient sert à créer 360 points autour d’une référence.

# A partir d'un point (unPoint) repere par ses coordonnees (longitude, latitude), cree un polygone circulaire de rayon donne (distance)
cercleDistant <- function(unPoint,distance){
# Creation de 360 points distincts sur le cercle de centre
# unPoint et de rayon distance
do.call("rbind",
lapply(0:360,function(bearing){
res <- gcDestination(
lon=unPoint[1],
lat=unPoint[2],
bearing=bearing,
dist=distance,
dist.units="km",
model="WGS84")
rownames(res) <- NULL
return(data.frame(res))
})
)
}

Comme l’idée n’est de donner un exemple, on isole le barycentre des Côtes d’Armor :

# L'indice du barycentre des Cotes d'Armor
indice.armor <- which(rownames(coordinates(barycentre))%in%"Cotes-Darmor")
# Les coordonnees du barycentre des Cotes d'Armor
barycentre.armor <- coordinates(barycentre[indice.armor])

Créons un cercle de rayon 40km autour de ce point :

# Le cercle de centre barycentre.armor et de rayon 40km
cercle.armor <- cercleDistant(barycentre.armor,distance=40)

Il ne reste plus qu’à repérer les points de la data.frame initiale (celle qui contient les 10 couples de coordonnées) qui sont dans ce cercle :

# Indice des points qui sont dans le cercle cercle.armor
indices.points.armor <- which(point.in.polygon(lesPoints$lon,lesPoints$lat,cercle.armor$long,cercle.armor$lat)==1)
# Les points dans le cercle
points.cercle.armor <- lesPoints[indices.points.armor,]

Et pour finir, un petit affichage graphique :

map(bretagne)
points(barycentre.armor,col="red",pch=19)
points(cercle.armor,col="dodger blue",cex=.5,pch=1)
points(lesPoints,col="gold",pch=19)
points(points.cercle.armor,col="purple",pch=19)

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Time limit is exhausted. Please reload CAPTCHA.