Tracer un cercle ou un carré sur la Terre
Kezhan SHI
02 Novembre 2015
1 Introduction
En lisant des articles sur internet l’autre jour, j’ai vu qu’une start-up dénommée what3words se développait rapidement en déployant une technique de géolocalisation. De prime abord, c’était assez farfelu : elle consiste tout simplement à donner trois mots, pour désigner un carré de 3m x 3m sur la surface de de la Terre. L’argument principal est que c’est un moyen beaucoup plus facile à retenir et à communiquer pour les hommes, que les coordonnées GPS avec plein de chiffres.
L’article parlait aussi de l’industrie géospatiale qui pourrait générer des millards d’euros d’ici quelques années. Du coup, je me suis repenché sur quelques techniques sur R :
- Géocoder une adresse quelconque sur la Terre
- Afficher les points sur un plan
- Compter les points à un endroit donné pour étudier leur répartition
Sur le dernier point, j’avais utilisé un simple maillage en fonction des latitude/longitude, pour étudier l’implantation géographique de certaines enseignes en France Et si on devait tracer des cercles ?
2 Eléments de base
Quelques éléments pour commencer :
Comme la Terre est approximativement sphérique, on utilise des coordonnées angulaires (latitude/longitude) pour situer un point sur la surface.
Comme on aime bien travailler sur une surface plate, on fait alors des projections. Il y a plusieurs types de projection.
La projection utilisée couramment est la projection Mercator. Google Maps ou OpenStreetMap utilisent cette projection (Plus précisément le Web Mercator, qui consiste à diminuer la déformation lors qu’on zoome). Néamoins, lors qu’on est sur une vue globale de la Terre, on voit exactement les déformations d’une projection Mercator classique. Ainsi, on peut lire de nombreux articles sur internet qui essaie de donner une meilleure vision des superficies des pays. Les exemples les plus frappants sont l’Afrique et le Groenland. Vous pouvez visiter le site thetruesize pour déplacer les pays.
Avec R, on peut utiliser le package leaflet pour visualiser facilement OpenStreetMap.
library(leaflet)
m <- leaflet(width="100%") %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addMarkers(lng=174.768, lat=-36.852, popup="The birthplace of R")
m # Print the map3 Distance du grand cercle
La distance à vol d’oiseau entre deux points se calcule avec la fonction suivante :
dis <- function(long1, lat1, long2, lat2) {
# long1=long1*pi/180
# lat1=lat1*pi/180
# long2=long2*pi/180
# lat2=lat2*pi/180
R <- 6371000
delta.long <- (long2 - long1)
delta.lat <- (lat2 - lat1)
a <- sin(delta.lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta.long/2)^2
c <- 2 * asin(sqrt(a))
d = R * c
return(d)
}On peut également exprimer les lat/lon d’un point final, à partir d’un point de départ, pour une distance et une direction donnée. Vous pourrez voir les routes orthodromique et loxodromique.
# Direction, l'angle à partir du nord, en tournant à droite
t=pi/2
# Point de départ
phi1=48*pi/180;lambda1=0*pi/180
# Distance
distance=seq(1,9000000,1200)
dest=function(t,phi1,lambda1,distance){
# Rayon de la Terre
R=6371000
delta=distance/R
phi=asin( sin(phi1)*cos(delta) + cos(phi1)*sin(delta)*cos(t))
lambda=lambda1 + atan2( sin(t)*sin(delta)*cos(phi1),
cos(delta)-sin(phi1)*sin(phi))
return(list(phi=phi,lambda=lambda))
}
m=leaflet(width="100%") %>%
addTiles() %>%
addPolygons(
lng=dest(t,phi1,lambda1,distance)$lambda*180/pi,
lat=dest(t,phi1,lambda1,distance)$phi*180/pi,fillColor = "transparent"
)%>%
addCircleMarkers(lng=lambda1*180/pi, lat=phi1*180/pi,radius=1,color="red")
mOn peut ainsi faire la tour de la Terre.
# Direction, l'angle à partir du nord, en tournant à droite
t=pi/4
# Point de départ
phi1=48*pi/180;lambda1=0*pi/180
# Rayon de la Terre
R=6371000
# Distance
distance=seq(1,2*pi*R,1200)
m=leaflet(width="100%") %>%
addTiles() %>%
addPolygons(
lng=dest(t,phi1,lambda1,distance)$lambda*180/pi,
lat=dest(t,phi1,lambda1,distance)$phi*180/pi,fillColor = "transparent"
)%>%
addCircleMarkers(lng=lambda1*180/pi, lat=phi1*180/pi,radius=1,color="red")
m4 Cercle sur la Terre
La définition du cercle est assez intuitive sur une surface sphérique : étant donné un point, appelé le centre du cercle, l’ensemble des points à distance égale du centre forment un cercle.
Pour tracer le cercle, on suit ainsi cette définition.
t=seq(0,2*pi,0.01)
phi1=0*pi/180;lambda1=0*pi/180
distance=9000000
R=6371000
m=leaflet(width="100%") %>%
addTiles() %>%
addPolygons(
lng=dest(t,phi1,lambda1,distance)$lambda*180/pi,
lat=dest(t,phi1,lambda1,distance)$phi*180/pi,
fillColor = "transparent"
)%>%
addCircleMarkers(lng=lambda1*180/pi, lat=phi1*180/pi,radius=1,color="red")
m5 Carré sur la Terre
On peut définir un carré sur une surface sphérique de plusieurs façons :
- On peut le définir par rapport à la distance au centre, comme pour un cercle.
- On peut dire que ce sont les 4 géodésiques qui relient les 4 points
- On peut également dire ce sont les 4 courbes qui relient les 4 points de façon à ce que les angles soient droits.
(D’ailleurs, ça me fait penser à une question, sur un triangle rectangle : étant donné un triangle rectangle, le côté opposé à l’angle rectangle a une longueur de 6, et la hauteur issue du sommet rectangle est de 5, quelle est l’aire du triangle ?)
5.1 Carré 1
Pour le tracer, il y a plusieurs façons, comme j’avais commencé le paramétrage angulaire, j’ai gardé ce paramétrage, en faisant varier la distance selon l’angle.
t=seq(0,2*pi,0.001)
phi1=30*pi/180;lambda1=0*pi/180
pc=rep(c(1/cos(seq(0,pi/4,0.001)),rev(1/cos(seq(0,pi/4,0.001)))[-1]),4)
distance=5000000*pc
R=6371000
m=leaflet(width="100%") %>%
addTiles() %>%
addPolygons(
lng=dest(t,phi1,lambda1,distance)$lambda*180/pi,
lat=dest(t,phi1,lambda1,distance)$phi*180/pi,
fillColor = "transparent"
)%>%
addCircleMarkers(lng=lambda1*180/pi, lat=phi1*180/pi,radius=1,color="red")
m5.2 Carré 2
t=seq(pi/4,7*pi/4,pi/2)
phi1=30*pi/180;lambda1=0*pi/180
distance=5000000*sqrt(2)
b=function(phi1,lambda1,phi2,lambda2){
dl=lambda2-lambda1
b= atan2(sin(dl)*cos(phi2),
cos(phi1)*sin(phi2)-sin(phi1)*cos(phi2)*cos(dl))
return (b)
}
coins=dest(t,phi1,lambda1,distance)
tlr=b(coins$phi[4],
coins$lambda[4],
coins$phi[1],
coins$lambda[1])
trb=b(coins$phi[1],
coins$lambda[1],
coins$phi[2],
coins$lambda[2])
blt=b(coins$phi[3],
coins$lambda[3],
coins$phi[4],
coins$lambda[4])
brl=b(coins$phi[2],
coins$lambda[2],
coins$phi[3],
coins$lambda[3])
distance_tlr=c(seq(1,dis(coins$lambda[1],coins$phi[1],
coins$lambda[4],coins$phi[4]),3000))
distance_tlb=c(seq(1,dis(coins$lambda[4],coins$phi[4],
coins$lambda[3],coins$phi[3]),3000))
polylng=c(dest(tlr,coins$phi[4],coins$lambda[4],distance_tlr)$lambda*180/pi,
dest(trb,coins$phi[1],coins$lambda[1],distance_tlb)$lambda*180/pi,
dest(brl,coins$phi[2],coins$lambda[2],distance_tlb)$lambda*180/pi,
dest(blt,coins$phi[3],coins$lambda[3],distance_tlr)$lambda*180/pi)
polylat=c(dest(tlr,coins$phi[4],coins$lambda[4],distance_tlr)$phi*180/pi,
dest(trb,coins$phi[1],coins$lambda[1],distance_tlb)$phi*180/pi,
dest(brl,coins$phi[2],coins$lambda[2],distance_tlb)$phi*180/pi,
dest(blt,coins$phi[3],coins$lambda[3],distance_tlr)$phi*180/pi)
m=leaflet(width="100%") %>%
addTiles() %>%
addCircleMarkers(
lng=coins$lambda*180/pi,
lat=coins$phi*180/pi,
radius=1,color="red"
)%>%
addPolygons(
lng=polylng,lat=polylat,
fillColor = "transparent"
)%>%
addCircleMarkers(lng=lambda1*180/pi, lat=phi1*180/pi,radius=1,color="red")
m6 Maillage
Comment on pourrait réaliser un maillage qui permet de découper la surface sphérique avec des aires égales ? Finalement, ce n’est pas une question si triviale. Si on est proche de l’équateur, on pourrait faire un maillage simple selon les latitude/longitude pour une surface pas trop grande. Si la surface est proche des pôles, et pas trop grande, on pourrait peut-être changer la base des coordonnées sphériques. Si on garde les coordonnées classiques, il faudrait décaler le quadriallage au fur et à mesure qu’on monte en latitude.
7 Références
Copyright © 2016 Blog de Kezhan Shi