# Moving the North Pole to the Equator

I am still working with @3wen on visualizations of the North Pole. So far, it was not that difficult to generate maps, but we started to have problems with the ice region in the Arctic. More precisely, it was complicated to compute the area of this region (even if we can easily get a shapefile). Consider the globe,

worldmap <- ggplot() +
geom_polygon(data = world.df, aes(x = long, y = lat, group = group)) +
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45)

and then, add three points in the northern hemisphere, and plot the associated triangle

P1 <- worldmap + geom_polygon(data = triangle, aes(x = long, y = lat, group = group),
fill ="blue", alpha = 0.6, col = "light blue", size = .8)+
geom_point(data = triangle, aes(x = long, y = lat, group = group),colour = "red")+

for some given projection, e.g.

coord_map("ortho", orientation=c(61, -74, 0))

This can be done with the following function

proj1=function(x=75){
triangle <- data.frame(long=c(-70,-110,-90*(x<90)+90*(x>90)),
lat=c(60,60,x*(x<90)+(90-(x-90))*(x>90)),group=1, region=1)
worldmap <- ggplot() +
geom_polygon(data = world.df, aes(x = long, y = lat, group = group)) +
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45)
P1 <- worldmap + geom_polygon(data = triangle, aes(x = long, y = lat, group = group),
fill ="blue", alpha = 0.6, col = "light blue", size = .8)+
geom_point(data = triangle, aes(x = long, y = lat, group = group),colour = "red")+
coord_map("ortho", orientation=c(61, -74, 0))
print(P1)
}

or

I am not sure if I understand why the projection of the triangle is not convex on the graph above, but say it’s not a big deal, here. Actually, our problem is that our interest is on regions (polygons, from a geometrical point of view) that do contain the North Pole. And here, it starts to be messy. I can easily move the upper point on the other side of the globe, but the polygon is not correct,

I do understand that it should be a problem, non-trivial, but it means that it should not be that simple to compute the area of a polygon (a region) that contains the North Pole. Which is exactly what we did observe in our computation. And I believe that one heuristic interpretation is related to the following graph

My skills in geometry are extremely poor. So do not expect that I will go through the code of the function that compute the area of a polygon ! Actually, my idea is the following : if the problem is that the North Pole is in the region, let’s consider some rotation, to shift the North Pole on the Equation. The code here is, from latitudes and longitude, to get new latitudes and longitudes, after a rotation around the y-axis (the North Pole will go down, along Greenwhich meridian) is

rotation=function(Z,theta){
lon=Z[,1]/180*pi; lat=Z[,2]/180*pi
x=cos(lon)*cos(lat)
y=sin(lon)*cos(lat)
z=sin(lat)
pt1=cbind(x,y,z)
M=matrix(c(cos(theta),0,-sin(theta),0,1,0,sin(theta),0,cos(theta)),3,3)
pt2=t(M%*%t(pt1))
lat=asin(pt2[,3])*180/pi
lon=atan2(pt2[,2],pt2[,1])*180/pi
return(cbind(lon,lat))}

With a rotation from $0$ (no change) to $\pi/2$ (the North Pole on the equator), we get

From now on, it is possible to compute the area of any region containing the North Pole ! One should simply apply the rotation function on all datebases generated from shapefiles (and then the opposite rotation to get a proper location) ! We can then compute the centroid of the ice region, for example,

r.glace=glace
r.glace[,1:2]=rotation(glace[,1:2],pi/2)
M=matrix(NA,length(unique(glace$id)),3) j=0 for(i in unique(glace$id)){j=j+1
Polyglace <- as(r.glace[glace$id==i,c("long","lat")],"gpc.poly") M[j,1]=area.poly(Polyglace) M[j,2:3]=centroid(r.glace[r.glace$id==i,c("long","lat")])
}
Z=c(weighted.mean(M[,2],M[,1]),weighted.mean(M[,3],M[,1]))
rotation(rbind(Z),-pi/2)[1,])

And we get

and below, we can visualize all the locations of the centroid of the ice region in the past 25 years

# Pavage et magie, rotations et symétries

(oui, mon fils voulait une rime dans le titre… promis on va tenter de faire une comédie musicale à partir de mes billets). La semaine passée, un article intéressant était publié (http://www.bbc.co.uk/news/…) sur sciences et magie. Cet article m’a fait penser à une réflexion que je m’étais faite lors de notre visite au MoMaths, cet hiver. Pour revenir sur l’histoire que j’avais raconté il y a quelques mois, sur la fin de la visite, les enfants avaient joué avec les pavages, à la Escher (les pièces étaient magnétiques, et s’assemblaient sur un grand tableau blanc)

Le soucis c’est qu’à aucun moment je n’ai pu voir d’explication du  pourquoi ça marche ? ou pour être plus explicite comment a été construit le singe pour que tout s’emboîte aussi bien. Car il y a un truc… La magie ça n’existe pas ! ce sont juste des maths (c’est un peu ce que racontent Persi Diaconis et Ron Graham dans Magical Mathematics, même si j’avoue ne pas l’avoir fini, faute de temps, et d’avoir laissé des chapitres en attente). Mon fils travaillant les rotations et les symétries, à l’école, j’ai voulu qu’on prenne le temps de revenir sur les pavages. Mais avant de me lancer avec lui, j’ai essayé de m’entraîner un peu tout seul (dimanche après midi, pendant le temps calme).

La difficulté (ou la magie) est qu’il y a des symétries dans tous les sens. Regardons un exemple simple (ou disons pas trop compliqué) avec des chinois (identifiés par les chapeaux dits chinois). En regardant rapidement, on note que le triangle ci-dessous est reproduit à l’infini, par un jeu de rotations. Il suffit donc d’avoir un tel triangle. Le soucis est que le découpage (parfaitement symétrique) n’est pas trivial (il est délicat de faire un personnage en faisant un découpage dans un triangle, même si c’est la méthode retenue dans l’activité mentionnée sur https://irem.univ-lille1.fr/…)

Le secret part de l’observation suivante: sur le graphique de gauche, en dessous, les trois points (qui sont des intersections des trois personnages, ou disons des trois couleurs) sont les sommets d’un triangle équilatéral. C’est là que la magie va opérer, car c’est cette figure (avec une autre, juxtaposée) que l’on va tenter de reproduire. En fait, un personnage (disons un petit chinois bleu) c’est deux triangles collés (on doit pouvoir parler de parallélogramme isocèle).

On ne le voit pas très bien, mais pourtant, au départ, on avait un parallélogramme isocèle (dont les contours sont tracés en noir) qui était  tout bleu. Bien entendu, on a découpé des parties et en les a envoyé ailleurs (par rotation autour d’un sommet). Le plus simple est de construire la figure, à la main. On part donc du parallélogramme suivant, tout bleu :

Oui, j’ai repris les fonctions de base utilisées l’autre jour pour expliquer R aux enfants,

source("http://freakonometrics.free.fr/RforKIDZ.R")

En fait, j’ai souhaité définir une fonction créant un polynôme à partir d’une liste de points (ou de sommets).

polyg=function(L,couleur="blanc",tour=couleur,ep=1,st=1,...){
b=unlist(L)
n=length(b)
polygon(b[seq(1,n,by=2)],b[seq(2,n,by=2)],col=couleurc(couleur),border=couleurc(tour),lwd=ep,lty=st,...)}

La fonction pour faire une rotation d’une liste de points est ici

rotation=function(L,angle=90,centre=L[[1]]){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]-centre[1],b[seq(2,n,by=2)]-centre[2])
Y=matrix(c(cos(angle/180*pi),sin(angle/180*pi),-sin(angle/180*pi),cos(angle/180*pi)),2,2)%*%t(X)
Z=Y+centre
L2=L
for(i in 1:(n/2)){L2[[i]]=Z[,i]}
return(L2)}

On a alors ici les quatre sommets

A=c(4,0)
B=c(4,10)
C=c(4+sqrt(3)*5,5)
D=c(4+sqrt(3)*5,15)

et on peut tracer le parallélogramme avec

L1=list(A,B,D,C)
polyg(L1,"bleu",tour="noir",ep=3)

On peut alors jouer à découper des parties, puis à les pivoter (on va faire des rotations, quoi). Par exemple, le prend le triangle en haut à droite, et je pivote de 240 degrés autour d’un des sommets du parallélogramme. C’est la figure ci-dessus à droite.  Le triangle, il est défini par le code suivant

E=c(10,10)
F=c(7,10+sqrt(3))
L1=list(B,E,F)
polyg(L1,tour="rouge",ep=3)

ce qui enlève la portion – en traçant un triangle blanc – et pour faire la rotation, j’utilise le code suivant

L2=rotation(L1,240)
polyg(L2,"bleu","rouge",ep=3)

(la seconde ligne permet de tracer le triangle bleu). Il va de soi qu’on peut découper des objets plus complexes que des triangles, voire non connexes. Par exemple, je peux découper la forme verte ci-dessous, puis la pivoter par rapport au sommet en haut à gauche, de 240 degrés.

G=c(10,2*sqrt(3)+10)
H=c(4+sqrt(3)*5,2*sqrt(3)+10)
F1=c(8,10+2/sqrt(3))
F2=c(8,10+4/sqrt(3))
L1=list(B,E,F1,F2,B,G,H,D,G)
polyg(L1,"vert",tour=NA)
L2=rotation(L1,240)
polygL(L2,"bleu",tour=NA)

(il n’est pas nécessaire de revenir au point B, mais peu importe). Ou la figure jaune, qui va pivoter autour du sommet en bas à droite, de 120 degrés.

F3=c(8,10-2/sqrt(3))
F4=c(8,10-4/sqrt(3))
G2=c(10,-2*sqrt(3)+10)
H2=c(4+sqrt(3)*5,-2*sqrt(3)+10)
L1=list(H,G,F2,F1,E,F3,F4,G2,H2,H)
polyg(L1,"jaune",tour=NA)
L2=rotation(L1,120,centre=C)
polygL(L2,"bleu",tour=NA)

Avec des ciseaux, ça sera plus simple, mais le but est aussi de voir que programmer, c’est facile. En gros, on a deux fonctions de base: une qui trace un polynôme et l’autre qui fait des rotations.

Si on fait les deux en même temps (à peu près), on enlève toute la partie en haut à droite, et on la répartie: une partie à gauche, et une partie en bas à droite.  On a alors une forme qui ressemble à un chinois. Non ? On peut aller creuser un peu sous les bras, pour améliorer les mains, mais on y est presque, me semble-t-il…

On peut continuer à améliorer, par exemple en enlevant des parties, puis en faisant des rotations par rapport à un des sommets initiaux. Comme le petit triangle jaune, qui va pivoter de 120 degrés par rapport au sommet en bas à droite. Puis le triangle vert par rapport au sommet en bas à gauche (on note qu’on pivote par rapport aux sommets du triangle équilatéral qu’on avait repéré, initialement)

S=c(8-.4,-2*sqrt(3)+10-.4)
L1=list(F4,G2,S,F4)
polyg(L1,"jaune")
L2=rotation(L1,120,centre=C)
polyg(L2,"bleu")

Bon, là on fait une pause, parce qu’on a, par morceaux, notre chinois ! On un paquet de lettre (voire un paquet de sommets qu’il va falloir nommer, suite à nos rotations). On les nomme avec un code qui ressemble à (oui, on aurait pu pu les nommer dès le départ, lors des rotations)

rE=rotation(list(E),240,centre=B)[[1]]
rF1=rotation(list(F1),240,centre=B)[[1]]
rH=rotation(list(H),240,centre=B)[[1]]

Une fois qu’on a identifié nos sommets, on a notre chinois,

chinois=list(B,rE,rF1,LG,S2,LF2,rH,A,r2H,r2G,r2S2,
r2F2,r2F1,r2E,rF3,rF4,r2S,rG2,rH2,C,H2,G2,S,F4,F3,E,B)

Bon, ensuite, reste à paver, ce qui est est le plus simple. Par exemple, on peut faire des rotations,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")

voire définir une translation,

translation=function(L,h){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]+h[1],b[seq(2,n,by=2)]++h[2])
L2=L
for(i in 1:(n/2)){L2[[i]]=X[i,]}
return(L2)}

puis translater notre chinois,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")
polyg(translation(r1chinois,15/sqrt(3),15),couleur="jaune")

Il manque l’autre, le vert. On peut vérifier qu’il se superpose parfaitement :

Bon, ben ça y est, on a constitué notre pavage. Reste à travailler un peu le code, essentiellement sur la forme, et on pourra jouer avec les enfants, à faire les lézards, des singes, et toutes sortes de pavages possibles, à la Escher. En notant qu’on peut utiliser d’autres formes de bases plus complexes, comme expliqué sur http://xavier.hubaut.info/…