{"word_count": 2637, "character_count": 23840}
Allez, un petit billet que je rêve de faire depuis des années, mais je viens enfin de trouver un prétexte, grâce à François de la formation Data Science pour l’Actuariat (à qui je vais emprunter l'application qui va suivre). Il y a plusieurs années, j'avais évoqué dans un billet les difficultés de travailler sur les données spatiales, surtout quand on est au pole. Le pole est un point de singularité dans la représentation classique (latitude, longitude) de la sphère. A l’époque, j'avais propose un bricolage simpliste (mettre les points de singularité dans un coin), mais il y a plus propre, à l'aide des quaternions d'Hamilton. Pour les présenter rapidement, et le lien avec les données spatiales, sur le globe terrestre, je vais reprendre la page wikipedia qui explique de manière incroyablement claire l’idée générale, et plus particulièrement les liens des quaternions avec les rotations (qui servent de base dans la représentation de la sphère qu'est le globe terrestre).Chaque rotation en dimension trois consiste à tourner d'un certain angle [latex]\alpha[/latex] autour d'un certain axe [latex]\vec{v}[/latex]. Pour un angle petit mais non nul (s'il est nul, on parle de rotation identité), l'ensemble des rotations possibles est une petite sphère entourant la rotation identité, où chaque point de la sphère représente un axe pointant dans une direction particulière. Des rotations d'angles de plus en plus grands s'éloignent progressivement de la rotation identité. Aussi, au voisinage de la rotation identité, l'espace abstrait des rotations ressemble à l'espace ordinaire en trois dimensions (qui peut également être vu comme un point central entouré de sphères de différents rayons).On peut assimiler les différentes directions à partir du pôle (c'est-à-dire les différents méridiens) aux différents axes de rotations et les différentes distances au pôle Nord aux différents angles&#160;: on a ainsi une analogie de l'espace des rotations. Mais la surface de la sphère est en deux dimensions alors que les axes de rotation utilisent déjà trois dimensions. L'espace des rotations est donc modélisé par une sphère de dimension 3 dans un espace à 4 dimensions (une hypersphère). Il faut alors penser la sphère ordinaire comme à une section de l'hypersphère (de la même façon qu'un cercle est une section de sphère). On peut prendre la section pour représenter, par exemple, uniquement les rotations d'axes dans le plan [latex]xy[/latex]. Et on peut légitimement penser maintenant aux rotations comme à des points de la sphère en dimension 4.Bon, rentrons dans les détails. On paramètre la surface d'une sphère à l'aide de deux coordonnées, comme la latitude et la longitude. Mais cela pose des soucis importants aux pôles (comme je l'avais note dans un ancien billet). Le théorème de la boule chevelue montre en fait qu'il n'existe aucun système de coordonnées à deux paramètres&#160; évitant cette dégénérescence. On va donc plonger la sphère dans l'espace à trois dimensions, en la paramétrant a l'aide de trois coordonnées cartésiennes (ici&#160;[latex]w[/latex],&#160;[latex]x[/latex] et&#160;[latex]y[/latex]). Par convention, on place le pôle Nord à [latex](w,y,z) = (1, 0, 0)[/latex], le pôle Sud à [latex](w,y,z) = (-1, 0, 0)[/latex] et l'équateur sera le cercle d'équations [latex]w = 0[/latex] et [latex]x^2+y^2=1[/latex]. Un point [latex](w,x,y)[/latex] de la sphère représente une rotation de l'espace ordinaire autour de l'axe horizontal dirigé par le vecteur [latex display="true"]{\displaystyle {\vec {v}}={\begin{pmatrix}x\\y\\0\end{pmatrix}}} [/latex] et d'angle [latex display="true"]{\displaystyle \alpha =2\cos ^{-1}w=2\sin ^{-1}{\sqrt {x^{2}+y^{2}}}}[/latex]C'est l’idée générale.Pour parler un peu des quaternions, pour rappel, un plan en dimension 2 peut être paramétré en utilisant les nombres complexes, en introduisant un symbole abstrait [latex] \mathbf {i}[/latex] qui vérifie la règle [latex] \mathbf {i}^2=-1[/latex]. On peut faire la meme chose en dimension 4, en introduisant des symboles abstraits [latex] \mathbf {i}[/latex], [latex] \mathbf {j}[/latex] et [latex] \mathbf {k}[/latex]. La partie imaginaire [latex] {\displaystyle b\mathbf {i} +c\mathbf {j} +d\mathbf {k} }[/latex] d'un quaternion se comporte comme un vecteur [latex display="true"]{\displaystyle {\vec {v}}={\begin{pmatrix}b\\c\\d\end{pmatrix}}} [/latex] d'un espace vectoriel à trois dimensions.Définissons le quaternion[latex display="true"] {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+\frac{\vec {v}}{\|v\|}\sin(\alpha /2)} [/latex]ou bien[latex display="true"] {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+{\vec {u}}\sin(\alpha /2)} [/latex] où [latex]{\displaystyle {\vec {u}}} [/latex] est un vecteur unitaire. Soit également [latex] {\displaystyle {\vec {v}}} [/latex] un vecteur ordinaire de l'espace en 3 dimensions, considéré comme un quaternion avec une coordonnée réelle nulle. On pourrait que le produit de quaternions[latex display="true"]{\displaystyle \mathbf{q}{\vec {\nu}}\mathbf{q}^{-1}} [/latex]renverrait le vecteur [latex] {\displaystyle {\vec {\nu}}} [/latex] tourné d'un angle [latex]{\displaystyle \alpha } [/latex] autour de l'axe dirigé par [latex] {\displaystyle {\vec {u}}} [/latex]. Et c'est effectivement ce qui se passe. Cette opération est connue comme la conjugaison par [latex] \mathbf {q}[/latex].Aussi, la multiplication de quaternions correspond à la composition de rotations, car si [latex] \mathbf {p}[/latex] et [latex] \mathbf {q}[/latex] sont des quaternions représentant des rotations, alors la rotation (conjugaison) par [latex] \mathbf {pq}[/latex] est[latex display="true"]{\displaystyle \mathbf{pq}q{\vec {\nu}}(\mathbf{pq})^{-1}=\mathbf{pq}{\vec {\nu}}\mathbf{q}^{-1}\mathbf{p}^{-1}=\mathbf{p}(\mathbf{q}{\vec {\nu}}\mathbf{q}^{-1})\mathbf{p}^{-1}} [/latex]ce qui revient à tourner (conjuguer) par [latex] \mathbf {q}[/latex], puis par [latex] \mathbf {p}[/latex].Le quaternion inverse d'une rotation correspond à la rotation inverse, car [latex]{\displaystyle q^{-1}(q{\vec {v}}q^{-1})q={\vec {v}}} [/latex]. Et assez naturellement, le carré d'un quaternion - noté [latex] \mathbf {q}^2[/latex] - correspond à la rotation de deux fois le même angle autour du même axe. Plus généralement, [latex] \mathbf {q}^{ {n}}[/latex] correspond à une rotation de [latex] {n}[/latex] fois l'angle autour du même axe que [latex] \mathbf {q}[/latex]. Par convention, on peut considérer un réel arbitraire [latex] {r}[/latex], ce qui permet de calculer des rotations intermédiaires de façon fluide entre des rotations de l'espace.Un petit exemple. Considérons la rotation [latex]f[/latex] autour de l'axe dirigé par [latex] {\displaystyle {\vec {v}}=\mathbf {i} +\mathbf {j} +\mathbf {k} } [/latex] et d'angle 120°, soit [latex] 2\pi/3[/latex].La norme de [latex] {\displaystyle {\vec {v}}} [/latex] est [latex] \sqrt{3} [/latex], le demi-angle est [latex]\pi/3[/latex] (ou 60°), le cosinus de ce demi-angle est [latex]1/2[/latex], et le sinus est [latex]\sqrt{3}/2[/latex]. Nous devons donc conjuguer avec le quaternion unitaire&#160;[latex display="true"] {\displaystyle \mathbf{q}=\cos {\frac {\pi }{3}}+\sin {\frac {\pi }{3}}\cdot {\frac {1}{\sqrt {3}}}{\vec {v}}}[/latex] soit[latex display="true"] \mathbf{q}={\frac {1}{2}}+{\frac {\sqrt {3}}{2}}\cdot {\frac {1}{\sqrt {3}}}{\vec {v}}={\frac {1}{2}}+{\frac {\sqrt {3}}{2}}\cdot {\frac {\mathbf {i} +\mathbf {j} +\mathbf {k} }{\sqrt {3}}} [/latex]qui peut finalement s’écrire simplement[latex display="true"]{\frac {1+\mathbf {i} +\mathbf {j} +\mathbf {k} }{2}}[/latex]Si [latex]f[/latex] est la fonction de rotation,[latex display="true"]{\displaystyle f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )=\mathbf{q}(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )\mathbf{q}^{-1}} [/latex]On peut prouver que l'on obtient l'inverse d'un quaternion unitaire simplement en changeant le signe de ses coordonnées imaginaires. Autrement dit[latex display="true"]{\displaystyle \mathbf{q}^{-1}={\frac {1-\mathbf {i} -\mathbf {j} -\mathbf {k} }{2}}} [/latex]et donc [latex]f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )[/latex] s'ecrit[latex display="true"] {\displaystyle {\frac {1+\mathbf {i} +\mathbf {j} +\mathbf {k} }{2}}(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} ){\frac {1-\mathbf {i} -\mathbf {j} -\mathbf {k} }{2}}} [/latex]En appliquant les règles ordinaires de calcul avec les quaternions, on obtient[latex display="true"] {\displaystyle f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )=c\mathbf {i} +a\mathbf {j} +b\mathbf {k} } [/latex]Ah oui, et de la même manière qu'on peut associer une matrice [latex]2\times2[/latex] à un nombre complexe[latex display="true"]z=a+b\mathbf{i} ~\rightarrow~\begin{pmatrix}a&#38;-b\\b&#38;a\end{pmatrix} [/latex]on peut associer une matrice [latex]4\times4[/latex]&#160;à un quaternion[latex display="true"]\mathbf{q}=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k} ~\rightarrow~\begin{pmatrix}\quad a&#38;\quad -b&#38;\quad -c&#38;\quad -d\\\quad b&#38;\quad a&#38;\quad -d&#38;\quad c\\\quad c&#38;\quad d&#38;\quad a&#38;\quad -b\\\quad d&#38;\quad -c&#38;\quad b&#38;\quad a\end{pmatrix} [/latex](on peut aussi passer par unE matrice [latex]2\times2[/latex] à coefficients complexes, mais ça ne servirait qu'à compliquer, ici).On a vu qu'on pouvait associer une rotation à un quaternion, et un quaternion à une matrice. Si la rotation est d'axe [latex]\vec {OM}[/latex], où [latex]O[/latex] est le centre de la terre, et [latex]M[/latex] un point sur la surface terrestre, décrit par sa latitude et sa longitude, avec pour angle [latex]\alpha[/latex] (comme décrit dans un billet sur stackoverflow), on a la fonction R suivantequat = function(lat,long,ang=NA){  n = length(lat)  lat = lat/180*pi  long = long/180*pi  x = cos(lat) * cos(long)  y = cos(lat) * sin(long)  z = sin(lat)  if (is.na(ang)){    Q = matrix(c(x,y,z,rep(0,n)), ncol = 4)  } else {    Q = matrix(c(sin(ang/2*pi/180) * c(x,y,z), cos(ang/2*pi/180)), ncol =4)  }  return(Q)}avec la réciproque, permettant de passer d'un quaternion à une rotation, avec une description de l'axe comme auparavant (un point sur la sphère - sur le globe terrestre) et un angle (coordonnées polaires)polaire = function(Q, digits=2) {Q = Q/norme(Q)ang = round(acos(Q[4])*2*180/pi,digits)n = norme(Q[1:3])x = Q[1]/ny = Q[2]/nz = Q[3]/nlat = asin(z) * 180/piif (z**2 == 1){long = 0}  else {phi = (x+1i*y)/sqrt(1-z**2)long = Im(log(phi)) * 180/pi }c(lat,long,ang)}à condition de définir au préalable la norme du quaternionnorme = function(Q){  sqrt(sum(Q**2))}On peut aussi définir le produits de quaternion (toutes les opérations sont décrites dans la page wikipedia), qui sera noté [latex]\otimes[/latex] par la suitepdt_quat = function(Q1,Q2){  Q=rep(0,4)  Q[1:3]  = Q1[4]*Q2[1:3]+Q2[4]*Q1[1:3]+pdt_vect(Q1[1:3],Q2[1:3])  Q[4] = Q1[4]*Q2[4] - pdt_scal(Q1[1:3],Q2[1:3])  return(Q)}mais aussi un produit scalairepdt_scal = function(M,N){  return(sum(M*N))}un produit vectorielpdt_vect = function(M,N){  return(c(M[2]*N[3]-M[3]*N[2],           M[3]*N[1]-M[1]*N[3],           M[1]*N[2]-M[2]*N[1]))}et finalement l'inverse du quaternion [latex]\mathbf{q}^{-1}[/latex]inv_quat = function(Q){  (c(0,0,0,2*Q[4])-Q)/norme(Q)}On peut aussi demander les coordonnées d'un point [latex]M'[/latex] obtenu comme transformation d'un point [latex]M[/latex] par la rotation [latex]\mathbf{q}[/latex]rotation = function(M, Q){polaire(pdt_quat(pdt_quat(Q,M),inv_quat(Q)))[1:2]}Maintenant, on va pouvoir passer aux choses sérieuses....Je l'ai évoqué en introduction, les quaternions peuvent permettre de contourner certains problèmes, comme manipuler des objets (comme la calotte glaciaire) qui sont situés autour du pole (qui est un point de singularité dans la représentation par coordonnées polaires). Une autre application, présentée par François, est celle du déplacement des plaques tectoniques. En particulier, on peut utiliser un fichier de rotations des plaques tectoniques entre aujourd’hui et une certaine date dans le passé (cette idée se retrouve dans le projet gplates programmé avec des librairies python). Ou plus généralement entre deux dates, [latex]t_1[/latex] et [latex]t_2[/latex]. Le fichier a notre disposition contient ainsi des quaternions [latex]\mathbf{q}^{P_0}_{t,P}[/latex] pour une date [latex]t[/latex] et une plaque [latex]P[/latex] (défini comme un polygone décrit par une collection de latitudes et de longitudes), ou le deplacement de la plaque est décrit relativement a la plaque [latex]P_0[/latex]. Pour des soucis de calculs, on va supposer qu'on peut interpoler linéairement les quaternions,[latex display="true"]\mathbf{q}^{P_0}_{t,P}=(1-\lambda)\mathbf{q}^{P_0}_{t_1,P}+\lambda\mathbf{q}^{P_0}_{t_2,P}[/latex]avec[latex display="true"]\lambda=\frac{t-t_1}{t_2-t_1}[/latex]et que récursivement, on peut composer les rotations, au sens ou[latex display="true"]\mathbf{q}^{P_3}_{t,P_1}\mathbf{q}^{P_3}_{t,P_2}\otimes\mathbf{q}^{P_2}_{t,P_1}[/latex]Aussi, a partir de notre fichier de rotations, on peut creer une fonction qui calcule l’ensemble des quaternions de rotation, pour les plaques plates données. A priori un quaternion [latex]\mathbf{q}^{P_0}_{t,P}[/latex] ne sera calculé qu'une foisprojecteur = function(t,plates,rot){  ll = length(plates)  Q0 = matrix(rep(0,4),ll,4)  for (i in seq(ll)){     cur_plate = i    Q = c(0,0,0,1)    while(cur_plate &#62; 0){  df = rot[rot$Start = t &#38; rot$plate == plates[cur_plate],]  if (dim(df)[1] &#62; 0) {  Q1 = Q0[cur_plate,]  cur_plate = 0   if (norme(Q1)==0){  Q_St = quat(df$lat_St, df$lon_St, df$ang_St)  Q_End = quat(df$lat_End, df$lon_End, df$ang_End)  periode = df$End - df$Start  pct = 0  if (periode &#62; 0) pct = min(max(0,(t - df$Start)/periode),1)  Q1 = Q_St + pct * (Q_End - Q_St)  cur_plate = which(plates==df$anchor)  if (length(cur_plate) == 0) cur_plate = 0  }   Q = pdt_quat(Q1,Q)  } else {cur_plate = 0}  }      Q0[i,] = Q  }  return(Q0)}On peut maintenant appliquer ces outils a des données. Ici, trois bases issues du site http://paleobiodb.org seront exploitées, pour visualiser ou les dinosaures vivaient :- une base des collections recensant les sites de fouilles (et en particulier leur géolocalisation)- une base d’occurrences recensant les spécimens trouvés, par collection- une base des spécimens décrivant les spécimens de dinosaures.site="http://paleobiodb.org/data1.2/"req=".txt?datainfo&#38;rowcount&#38;max_ma=999&#38;min_ma=0"limit = "" #"&#38;limit=100"names = c("colls/list", "occs/taxa", "occs/list")destfile = rep('',length(names))for (i in 1:length(names)) {  destfile[i] = paste0("data/",sub("/","_",names[i]),".csv")  download.file(paste0(site,names[i],req,limit),destfile=destfile[i])}Pour les données de collectioncollection = read.csv(destfile[1],skip=17,sep = ",", header = TRUE)coll = select(collection, c(collection_no, lng, lat, max_ma, min_ma))remove(collection)names(coll) = c('no', 'lng', 'lat', 'max_ma', 'min_ma')et pour les données d’occurrenceoccurence = read.csv(destfile[3],skip=17, sep = ",", header = TRUE)occ = select(occurence, c(occurrence_no, collection_no, accepted_no))remove(occurence)names(occ) = c('no', 'coll_no', 'taxo_no')Les espèces que nous étudierons ici sont les dinosaures de l’ordre des Ornithischia et des Saurischia (qui incluent les grandes familles classiques de dinosaures - a ce que j'ai pu comprendre) taxonomie = read.csv(destfile[2],sep = ",", skip=20, header = TRUE)    taxo = select(taxonomie,c(orig_no, accepted_rank, accepted_name, parent_no, container_no))    remove(taxonomie)    names(taxo) = c('N0', 'rang0', 'nom0', 'parent', 'container')    taxo = taxo[!is.na(taxo$N0),]    taxo["N1"]=as.character(taxo$container)    taxo[taxo$container=='',]$N1 = as.character(taxo[taxo$container=='',]$parent)    taxo$container = as.factor(taxo$container)    taxo$parent = NULL    taxo$container = NULL    t0 = taxo    i = 0    while(dim(taxo[!is.na(taxo[paste0("N",i)]),])[1] &#62; 0){      i = i+1      colnames(t0) = c(paste0('N',i),paste0('rang',i), paste0('nom',i), paste0('N',i+1))      taxo = merge(taxo,t0,all.x=TRUE)      if (i==10) break    }    t1 = select(taxo,N0)    for (niveau in levels(taxo$rang0)){      if (niveau != ""){        t1[niveau]=""        for (j in seq(0,i)){          test = which(taxo[paste0("rang",j)]==niveau)           t1[test,niveau] = as.character(taxo[test,paste0("nom",j)])        }        }    }    taxo=t1[t1["unranked clade"%in%c("Ornithischia","Saurischia")],]    head(taxo,5)On peut alors fusionner nos basesM1 = merge(occ,taxo, by.x=c('taxo_no'),by.y=c('N0'), all=FALSE)M2 = merge(M1,coll,by.x=c('coll_no'),by.y=c('no'), all.x=TRUE)paste(dim(M2)[1], "specimens étudiés")ce qui donne 9771 spécimenshead(M2,5)coll_no taxo_no     no class           family             genus1    5195   55999 373398           Nodosauridae      Pawpawsaurus2   10755   55580 130209       Chaoyangsauridae    Chaoyangsaurus 3   10760   38561 144305        Dromaeosauridae                  4   10764   66066 130295        Caudipterygidae       Caudipteryx5   10764   66068 130294                        Protarchaeopteryx  infraclass kingdom        order   phylum                   species1                                 Chordata    Pawpawsaurus campbelli2                                 Chordata     Chaoyangsaurus youngi3                    Avetheropoda Chordata                          4                    Avetheropoda Chordata          Caudipteryx zoui5                    Avetheropoda Chordata Protarchaeopteryx robusta  subclass subfamily subgenus suborder subspecies superclass superfamily1                                                                       2                                                                       3                                                                       4                                                                       5                                                                         superphylum tribe unranked clade      lng      lat max_ma min_ma1                     Ornithischia -97.3000 32.86667  105.3  99.602                     Ornithischia 123.9667 42.93330  150.8 132.903                       Saurischia  21.0500 46.11667   70.6  66.004                       Saurischia 120.7333 41.80000  130.0 122.465                       Saurischia 120.7333 41.80000  130.0 122.46Voila pour les dinosaures. On peut maintenant chercher des informations sur les plaques tectoniques,chemin = "data/Shapefile"download.file('https://www.earthbyte.org/webdav/ftp/earthbyte/GPlates/SampleData_GPlates2.0/Individual/FeatureCollections/Coastlines.zip', 'data/coastlines.zip')coast_file = 'Matthews_etal_GPC_2016_Coastlines'unzip(zipfile='data/coastlines.zip', exdir= chemin, junkpaths = TRUE)continents = readOGR(dsn=chemin,layer=coast_file,verbose=TRUE)Matthews et al. (2016) a mis en ligne un fichier de rotations simulant la dérive des plaquesdownload.file('https://www.earthbyte.org/webdav/ftp/earthbyte/GPlates/SampleData_GPlates2.0/Individual/FeatureCollections/Rotations.zip', 'data/rot.zip')rot_file = 'Matthews_etal_GPC_2016_410-0Ma_GK07.rot'unzip(zipfile='data/rot.zip', files = c(paste0('Rotations/',rot_file)), exdir= 'data', junkpaths = TRUE)rot = paste0("data/",rot_file)On va corriger quelques anomaliesx = readLines(rot)y = gsub( "!101 !", "!", x )cat(y, file=rot, sep="\n")remove(x,y)et on charge les données en mémoire (pour faire ensuite notre visualisation)rot_file2 = read.csv(file=rot,header=FALSE,sep='', comment.char = '!')ll = dim(rot_file2)[1]rot_file3 = cbind(rot_file2[1:ll-1,],rot_file2[2:ll,])names(rot_file3) = c('plate','Start','lat_St', 'lon_St', 'ang_St', 'anchor','plate2','End','lat_End', 'lon_End', 'ang_End', 'anchor2') rot_file = rot_file3[rot_file3$plate==rot_file3$plate2,]rot_file$plate2 = NULLrot_file$anchor2 = NULLPLATES = sort(unique(rot_file$plate))On ne va garder que les dinosaures qui peuvent etre rattaches&#160;à une plaque tectoniqueX = M2%&#62;%select(lng, lat)Y = SpatialPoints(X,proj4string = continents@proj4string)plaques = over(Y,continents)$PLATEID1filtre = which(!is.na(plaques))X0       = X[filtre,]NBX = dim(X0)[1]print(paste0(NBX, " spécimens retenus"))ce qui laisse quand même 9617 spécimensPERIOD = M2%&#62;%select(max_ma, min_ma)PERIOD  = PERIOD[filtre,]plaques = plaques[filtre]plaque_id = rep(0,NBX)for (j in seq(1,NBX)){    plaque_id[j]=which(PLATES==plaques[j])}dataX = M2[filtre,]dataX["plaques"] = plaquesOn y est presque... on va maintenant remonter de -250 millions d'annees a aujourd'hui, en faisant des bonds de 10 millions d'anneesTMAX = 250TMIN = 0PAS = 10Pour toutes ces dates, on calcule les quaternionsROT = array(rep(0,4),c(TMAX/PAS,length(PLATES),4))for (t in seq(1,TMAX/PAS)){  ROT[t,,] = projecteur(t*PAS, PLATES, rot_file)}QM = list()plaque=list()xy = list()ll = length(continents@polygons)for (i in seq(1,ll)){  M = continents@polygons[[i]]@Polygons[[1]]@coords  xy[[i]] = M  QM[[i]] = quat(M[,2],M[,1])  plaque[[i]] = which(PLATES==continents$PLATEID1[i])}QX = quat(X0[,2],X0[,1])On va ensuite projeterprojete=list()Xt = list()X = X0cpt = 0for (TIME in seq(TMAX,TMIN,-PAS)){  # setTxtProgressBar(pb, -TIME)  cpt = cpt+1  projete[[cpt]] = continents  Xt[[cpt]] = X0  if (TIME &#62; 0 ) {    for (i in seq(1,ll)){      M = xy[[i]]      for (j in seq(1,dim(M)[1])){        M[j,] = rev(rotation(QM[[i]][j,],ROT[TIME/PAS,plaque[[i]],]))      }      inf_180 = which(M[,1] &#60; -90); inf_180_ = length(inf_180)      sup_180 = which(M[,1] &#62; 90); sup_180_ = length(sup_180)      if (inf_180_ &#62; 0 &#38; sup_180_ &#62; 0) {        if(sup_180_&#62;inf_180_) {          M[inf_180,] = t(t(M[inf_180,]) + c(360,0))        } else { M[sup_180,] = t(t(M[sup_180,]) - c(360,0))}      }       projete[[cpt]]@polygons[[i]]@Polygons[[1]]@coords = M    }      # setTxtProgressBar(pb, -TIME + PAS/2)    for (j in seq(1,NBX)){      if (PERIOD$max_ma[j] &#62; TIME &#38; PERIOD$min_ma[j] &#60;= TIME){        X[j,] = rev(rotation(QX[j,],ROT[TIME/PAS,plaque_id[j],]))      } else{        X[j,] = c(NA, NA)      }    }    filtre = which(!is.na(X$lng))    if (length(filtre) &#62; 0){      Xt[[cpt]]=SpatialPointsDataFrame(coords = X[filtre,], data = dataX[filtre,])      } else {      Xt[[cpt]]=X    }  }}On peut faire en première carte, 70 millions d’années avant notre eret = 1 + (TMAX-70)/PASleaflet(options = leafletOptions(minZoom = 1)) %&#62;%    addPolygons(data=projete[[t]], weight=2) %&#62;%   addMarkers(data = Xt[[t]],              popup = ~paste(sep = "", paste("espèce :", Xt[[t]]$species),                            paste("genre :", Xt[[t]]$genus),                            paste("famille :", Xt[[t]]$family),                            paste("ordre : ", Xt[[t]]["unranked clade",]) ),             icon = dinoIcon)Je mets ici une copie d’écran du leaflet ainsi créé(l’idée est qu'on peut zoomer, ce qui rend l'analyse plus interactive)Mais on peut aussi aller 200 millions d’années avant notre eret = 1 + (TMAX-200)/PASleaflet(options = leafletOptions(minZoom = 1)) %&#62;%    addPolygons(data=projete[[t]], weight=2) %&#62;%   addMarkers(data = Xt[[t]],              popup = ~paste(sep = "", paste("espèce :", Xt[[t]]$species),                            paste("genre :", Xt[[t]]$genus),                            paste("famille :", Xt[[t]]$family),                            paste("ordre : ", Xt[[t]]["unranked clade",]) ),             icon = dinoIcon)Amusant, non ? en tout cas, merci François pour cette jolie application des quaternions ! Et merci d'avoir suggéré d'utiliser autre chose que des points rouges sur une carte !dinoIcon = makeIcon(iconUrl = "https://www.ludeek.com/wp-content/uploads/2015/03/uploadfsdfsdf1426350179.1426350368774.png",                    iconWidth = 30, iconHeight = 50,                    iconAnchorX = 15, iconAnchorY = 25)