Introduction
La question est très simple, si vous montez au sommet de la Tour Eiffel à Paris (pas la Tour Eiffel à Las Vegas) et que vous regardez vers l'ouest. Quel est le premier pays que votre regard "croisera"?
Pour vous aider un peu, en regardant vers l'ouest votre regard va donc voir l'ouest parisien, mais si l'on imagine plus loin derrière l'horizon on traversera la Bretagne puis l'Atlantique jusqu'à arriver aux côtes du continent américain. Le premier pays (autre que le France) que votre regard croisera est donc un pays d'Amérique. Mais lequel ? Réfléchissez 1 minute avant de continuer la lecture afin d'être bien certain de votre réponse.
L'objectif de cet article est de calculer, étant donné un point de référence (ici la Tour Eiffel) quel est le premier pays (autre que la France évidemment) que votre regard croisera selon la direction dans laquelle vous regardez. Vers l'est on doit certainement tomber sur l'Allemagne, vers le sud-ouest surement sur l'Espagne, etc... C'est plutôt simple à visualiser pour les pays limitrophes mais nous verrons que les choses sont moins intuitives lorsqu'il y a de grandes distances à traverser.
Méthode
Pour effectuer les calculs il nous faut les positions et la forme de chaque pays. Il faut alors utiliser des fichiers de type shapefile qui représentent chaque pays sous forme de polygones plus ou moins détaillés (selon la précision de la carte) et en définissent les contours sous forme de coordonnées. L'analyse consiste alors à vérifier si la droite (la direction) suivie par notre regard traverse le polygone en question et si tel est le cas, à quelle distance (puisque nous cherchons le "polygone" le plus proche).
La méthode retenue consiste à choisir évidemment un point de départ et une direction (un angle). Ensuite, deux possibilités s'offrent a vous :
- Transformer la représentation du fichier de polygones en coordonnées azimutales centrées sur votre point de départ (ce qui conserve les angles) et ensuite étudier les intersections entre votre droite et les polygones comme vous le feriez dans un plan orthogonal en utilisant les équations de droite et les intersections. C'est la méthode retenue. La représentation azimutale est par exemple utilisée pour le logo des nations unies :
La représentation est ici centrée sur le pôle nord mais il est possible facilement de centrer cette représentation sur n'importe quel point. A partir de ce point, les angles sont parfaitement conservés.
- Vous pouvez aussi garder les coordonnées initiales, en revanche dans ce cas-là, l'intersection entre la droite formée par votre regard et le polygone ne pourra se résoudre qu'en utilisant des équations beaucoup plus complexes à cause de la rotondité de la Terre. Cela reste faisable néanmoins.
Étape 1 : On charge le shapefile (sauvé préalablement dans mon working directory)
# On charge les shapefiles des pays library('rgdal')
library('ggplot2')
shp10 <- readOGR(dsn = "ne_10m_admin_0_countries", layer = "ne_10m_admin_0_countries") shp50 <- readOGR(dsn="ne_50m_admin_0_countries", layer="ne_50m_admin_0_countries") shp110 <- readOGR(dsn = "ne_110m_admin_0_countries", layer = "ne_110m_admin_0_countries") shp2 <- fortify(shp50) colnames(shp2)[1:2] = c("Longitude", "Latitude") # On a 10654 points de shapes appartenant a 176 Pays pour shp110 # On a 99599 points pour shp50 # on a 548121 points pour shp10
Finalement c'est le fichier de plus petite taille qui a été retenu. Il est composé de 10 000 points ce qui constitue une résolution malgré tout suffisante pour notre analyse. Le temps de calcul étant de toute façon proportionnel au nombre de points, on ne veut pas y passer trop de temps. Malgré tout il n'est pas nécessaire d'avoir une très grande précision dans le détail des polygones puisque le point qui nous intéresse est uniquement de savoir si on croise ou pas le pays. Il sera toujours possible d'ergoter sur quelques cas particuliers mais dans la grande majorité le résultat sera suffisamment satisfaisant.
Étape 2 : Définition des villes d’intérêt, on crée juste un data.frame avec les informations initiales. On parcourra une à une chaque ligne.
VILLE=data.frame(NOM=c("PARIS","BREST","CALAIS","MARSEILLE","STRASBOURG","AJACCIO"), LON=c(2.2069771,-4.5095871,1.6557295,5.2400702,7.6918828,8.6356297), LAT=c(48.8587741,48.3637114,50.9718019,43.2805546,48.5692058,41.9227051))
Et enfin l’étape 3, la boucle sur chacune des villes. J'annonce à l'avance que ce code peut être amélioré. Il est très basique dans le sens que l'on teste toutes les paires de coordonnées d'un polygone ainsi que la majorité des polygones. Il s'agit évidemment de la partie la plus consommatrice de temps. Il existe énormément de polygones qui n'auraient pas besoin d'être testés (par exemple, depuis la France, tous les polygones représentant les pays au sud du Maghreb ne seront jamais "visibles" car toujours "bloqués" par les pays d'Afrique du Nord. Idem pour la plupart des pays asiatiques. En fait c'est le cas pour tous les pays "enclavés" (c'est à dire entourés d'autres pays). Seuls les pays ayant une frontières terrestre avec la France ou bordant la mer ou un océan peuvent nous convenir. Tous les autres pourraient être éliminés.
Il existe d'autre versions de la résolutions de ce problème que vous pourrez trouver un peu partout sur le net). Par exemple là c'est du Python. Ici, ça sera en R.
for(N in 1:nrow(VILLE)){ print(VILLE[N,1]) Longitude0 <- VILLE[N,"LON"] Latitude0 <- VILLE[N,"LAT"] # Transfo azimutale (a faire a chaque fois que l'on change de centre) shp110a <- spTransform(shp110, CRS(sprintf("+proj=aeqd +lat_0=%f +lon_0=%f",Latitude0,Longitude0))) # Conversion de l'objet shpa <- fortify(shp110a) colnames(shpa)[1:2] = c("Longitude", "Latitude") # On splitte par groupe # attention un pays peut appartenir a +ieurs gpe : expl france + corse + guyane) sshpa <- split(shpa,as.numeric(shpa$group)) # On va stocker le pays le plus proche pour chaque valeur d'angle BESTCOUNTRY = c() # Liste des angles a parcourir, ici tous les 1° SEQ = seq(0,360,by=1) t0 <- Sys.time() for (angle in SEQ){ cat(sprintf("\n-----------\nAngle %i\n-----------\n",angle)) # On convertit l'angle qui est en degree theta = pi * angle / 180 RESangle = array(Inf,dim=length(sshpa)) # On parcourt les pays 1 a 1 pylogones 1 a 1 for(i in 1:length(sshpa)){ # On skippe l'id de la france if(i%in% c(120,121,122)) next dat <- unique(sshpa[[i]][,1:2]) colnames(dat) = c("long","lat") # En fonction de l'angle, certain pays n'ont pas a etre testés # Par exemple si on regarde vers le NE # Tous les points des coordonnées negatives n'ont pas a etre analysé # Réduit le nb de calculs a faire if(angle%%360 >=0 & angle%%360 <=180) if(all(dat$lat<0)) next if(angle%%360 >=180) if(all(dat$lat>=0)) next if((angle-90)%%360 >=0 & (angle-90)%%360 <=180) if(all(dat$long>0)) next if((angle-90)%%360 >=180) if(all(dat$long<=0)) next # Si le pays est visible on l'affiche # On parcourt toutes les paires de points Dist = Inf tmp1=tmp2 =c () for(j in 1:nrow(dat)){ if(j==nrow(dat)){ p1 = dat[j,c("long","lat")] p2 = dat[1,c("long","lat")] }else{ p1 = dat[j,c("long","lat")] p2 = dat[j+1,c("long","lat")] } # On calcule l'intersection entre notre point (0,0) et le segment p1 p2 # Equation de la droite passant par 0,0 et de pente "angle" a1 = tan(theta) b1 = 0 # Equation du segment formé par les deux points consecutifs a2 = (p2[2]-p1[2]) / (p2[1]-p1[1]) b2 = p1[2] - a2 * p1[1] # Intersection des deux droites xi = (b2-b1)/(a1-a2) yi = a1 * xi + b1 # On teste si cette intersection se trouve a l'interieur du segment isIN = xi < max(p1[1],p2[1]) & xi > min(p1[1],p2[1]) & yi < max(p1[2],p2[2]) & yi > min(p1[2],p2[2]) # L'intersection est peut-etre sur le segment mais il faut verifier qu'on soit dans le bon sens # expl : avec un angle a 180 degree on s'attend a ce que xi < 0 if(isIN){ if(angle%%360 >=0 & angle%%360 <=180) if(round(yi,3)<0) isIN = FALSE if(angle%%360 >=180) if(round(yi,3)>0) isIN = FALSE if((angle-90)%%360 >=0 & (angle-90)%%360 <=180) if(round(xi,3)>0) isIN = FALSE if((angle-90)%%360 >=180) if(round(xi,3)<0) isIN = FALSE } # Et on calcule la distance de cette intersection a 0,0 if(isIN) if(sqrt(xi^2 + yi^2) < Dist) Dist = sqrt(xi^2 + yi^2) tmp1[j]=isIN tmp2[j]=unlist(sqrt(xi^2 + yi^2)) } RESangle[i] = as.numeric(Dist) } # On regarde quel est l'id du polygone le plus proche # et que l'on a traversé gpebest = order(RESangle)[1] # On retrouve le pays auquel il appartient idbest = as.numeric(sshpa[[gpebest]]$id[1]) print(shp110$ADMIN[idbest+1]) BESTCOUNTRY[match(angle,SEQ)] = shp110$ADMIN[idbest+1] print( Sys.time() - t0) } save(BESTCOUNTRY,SEQ,file=sprintf("dat%s.RData",VILLE$NOM[N])) }
Il s'agit donc d'un code assez "bourrin" qui cherche les coordonnées de l'intersection entre la droite (en fait, une demi-droite) définie par notre regard et une arête d'un polygone. Si cette intersection se trouve bien à l’intérieur du segment définit par l'arête, alors notre regard a bien croisé cette arête et donc ce polygone, il ne reste plus qu'à en calculer la distance et stocker tout cela. Sans parallélisation du code et en utilisant la cartographie à 10 000 points, il faut environ 1h30 de calcul pour obtenir les 360 résultats (un résultat par degré) donc il ne faut pas être trop pressé. La partie importante de ce code réside dans la conversion en coordonnées azimutale, ce qui implique que le point où vous vous trouvez se situe forcément aux coordonnées (0,0) et que les coordonnées des polygones des pays ont été redéfinies en fonction de cette nouvelle origine. Cela simplifie énormément les calculs et comme la fonction de transformation existe déjà cela enlève la partie la plus complexe de l'analyse.
Résultats
Une fois obtenu le résultat, la représentation graphique est laissé à l'inspiration de chacun. J'en ai fait une version assez basique mais très colorée. (De toute façon peu importe les couleurs que vous choisirez, il y a 28 pays différents, donc 28 couleurs a représenter, inutile d’espérer de trouver une palette de couleurs avec 28 teintes réellement différentiables a tous les coups)
 |
| Regarde maman, on voit vachement bien le Venezuela d'ici |
Dans cette représentation, le nord est en haut et le reste des points cardinaux est à l'avenant.
Analyse
En regardant ce premier résultat pour Paris, on arrive bien a imaginer qu'en se postant au sommet de la tour Eiffel et en regardant vers l'est le prochain pays que notre regard "rencontrera" sera l'Allemagne, ou la Suisse, le Luxembourg et la Belgique selon l'angle choisit. Je met des guillemets à "rencontrera" parce qu'on est bien d'accord qu'arrivé à cette étape vous vous doutez bien qu'on ne "voit" pas vraiment l'autre pays depuis la tour Eiffel, déjà parce que la Terre étant ronde la visibilité maximale autour de la tour Eiffel et d'environ 65 km et aussi parce que même si la Terre était plate, l'air de Paris serait trop pollué pour espérer voir quoique ce soit. (En revanche si la Terre est plate les calculs doivent surement être plus simples)
Revenons à nos moutons. Si vous êtes habitués aux cartes et aux projections traditionnelles, vous devez surement vous dire, comme moi, qu'en regardant vers l'ouest depuis Paris on devrait tomber sur le Canada ou les USA, en tout cas pas la République Dominicaine !!!
Preuve (Google Maps) :
Il existe donc deux possibilités : soit mon code est faux, soit la Terre n'est pas plate. Comme mon code n'est (évidemment) pas faux, la Terre n'est donc pas plate (CQFD) et il faut donc utiliser Google Earth pour tracer une vraie ligne droite sur le globe.
Preuve (Google Earth) :
 |
| On part bien vers l'ouest depuis Paris |
 |
| Et paf, on est en République Dominicaine ! |
Étonnant non? Vous ne l'aviez pas vue venir celle-là? On peut ensuite regarder les différents résultats obtenus selon la ville de départ sélectionnée :
Pour Brest :
C'est l'endroit d'où on voit le mieux le Brésil.
Pour Calais :
Apparemment, de Calais on semble très bien voir le Royaume-Uni et la Belgique (et aussi accessoirement l'Espagne mais il me semble que c'est un peu plus loin). Il est amusant de constater que toutes ces villes du Nord peuvent apercevoir (à des degrés divers) des pays qui semblent très au sud comme le Brésil, ce qui n'est possible que parce que la vue n'est pas obstruée par l'Espagne lorsque l'on se situe suffisamment au nord.
Pour Marseille :
Côté est c'est évidemment l'Italie qui prédomine, au sud l'Algérie, pour la partie est / sud-ouest on a aussi l'Espagne. Ces trois pays représentent respectivement 127°, 60° et 57° soit donc un total de 244°, les 2/3 du paysage donc.
Pour Strasbourg :
Le Pac-Man allemand.
Pour Ajaccio :
Et pour notre dernier essai, le Pac-Man italien, a noter une petite incursion de la Libye.
Aucun commentaire:
Enregistrer un commentaire