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 \alpha autour d’un certain axe \vec{v}. 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 : 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 xy. 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 é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 w, x et y). Par convention, on place le pôle Nord à (w,y,z) = (1, 0, 0), le pôle Sud à (w,y,z) = (-1, 0, 0) et l’équateur sera le cercle d’équations w = 0 et x^2+y^2=1. Un point (w,x,y) de la sphère représente une rotation de l’espace ordinaire autour de l’axe horizontal dirigé par le vecteur {\displaystyle {\vec {v}}={\begin{pmatrix}x\\y\\0\end{pmatrix}}} et d’angle {\displaystyle \alpha =2\cos ^{-1}w=2\sin ^{-1}{\sqrt {x^{2}+y^{2}}}}
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 \mathbf {i} qui vérifie la règle \mathbf {i}^2=-1. On peut faire la meme chose en dimension 4, en introduisant des symboles abstraits \mathbf {i}, \mathbf {j} et \mathbf {k}. La partie imaginaire {\displaystyle b\mathbf {i} +c\mathbf {j} +d\mathbf {k} } d’un quaternion se comporte comme un vecteur {\displaystyle {\vec {v}}={\begin{pmatrix}b\\c\\d\end{pmatrix}}} d’un espace vectoriel à trois dimensions.
Définissons le quaternion {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+\frac{\vec {v}}{\|v\|}\sin(\alpha /2)} ou bien {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+{\vec {u}}\sin(\alpha /2)} où {\displaystyle {\vec {u}}} est un vecteur unitaire. Soit également {\displaystyle {\vec {v}}} 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{\displaystyle \mathbf{q}{\vec {\nu}}\mathbf{q}^{-1}} renverrait le vecteur {\displaystyle {\vec {\nu}}} tourné d’un angle {\displaystyle \alpha } autour de l’axe dirigé par {\displaystyle {\vec {u}}} . Et c’est effectivement ce qui se passe. Cette opération est connue comme la conjugaison par \mathbf {q}.
Aussi, la multiplication de quaternions correspond à la composition de rotations, car si \mathbf {p} et \mathbf {q} sont des quaternions représentant des rotations, alors la rotation (conjugaison) par \mathbf {pq} est{\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}} ce qui revient à tourner (conjuguer) par \mathbf {q}, puis par \mathbf {p}.
Le quaternion inverse d’une rotation correspond à la rotation inverse, car {\displaystyle q^{-1}(q{\vec {v}}q^{-1})q={\vec {v}}} . Et assez naturellement, le carré d’un quaternion – noté \mathbf {q}^2 – correspond à la rotation de deux fois le même angle autour du même axe. Plus généralement, \mathbf {q}^{ {n}} correspond à une rotation de {n} fois l’angle autour du même axe que \mathbf {q}. Par convention, on peut considérer un réel arbitraire {r}, 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 f autour de l’axe dirigé par {\displaystyle {\vec {v}}=\mathbf {i} +\mathbf {j} +\mathbf {k} } et d’angle 120°, soit 2\pi/3.
La norme de {\displaystyle {\vec {v}}} est \sqrt{3} , le demi-angle est \pi/3 (ou 60°), le cosinus de ce demi-angle est 1/2, et le sinus est \sqrt{3}/2. Nous devons donc conjuguer avec le quaternion unitaire {\displaystyle \mathbf{q}=\cos {\frac {\pi }{3}}+\sin {\frac {\pi }{3}}\cdot {\frac {1}{\sqrt {3}}}{\vec {v}}} soit \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}}} qui peut finalement s’écrire simplement{\frac {1+\mathbf {i} +\mathbf {j} +\mathbf {k} }{2}}
Si f est la fonction de rotation,{\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}}
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{\displaystyle \mathbf{q}^{-1}={\frac {1-\mathbf {i} -\mathbf {j} -\mathbf {k} }{2}}} et donc f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} ) s’ecrit {\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}}} En appliquant les règles ordinaires de calcul avec les quaternions, on obtient {\displaystyle f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )=c\mathbf {i} +a\mathbf {j} +b\mathbf {k} } Ah oui, et de la même manière qu’on peut associer une matrice 2\times2 à un nombre complexez=a+b\mathbf{i} ~\rightarrow~\begin{pmatrix}a&-b\\b&a\end{pmatrix} on peut associer une matrice 4\times4 à un quaternion\mathbf{q}=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k} ~\rightarrow~\begin{pmatrix}\quad a&\quad -b&\quad -c&\quad -d\\\quad b&\quad a&\quad -d&\quad c\\\quad c&\quad d&\quad a&\quad -b\\\quad d&\quad -c&\quad b&\quad a\end{pmatrix} (on peut aussi passer par unE matrice 2\times2 à 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 \vec {OM}, où O est le centre de la terre, et M un point sur la surface terrestre, décrit par sa latitude et sa longitude, avec pour angle \alpha (comme décrit dans un billet sur stackoverflow), on a la fonction R suivante
quat = 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)
} |
quat = 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]/n
y = Q[2]/n
z = Q[3]/n
lat = asin(z) * 180/pi
if (z**2 == 1){
long = 0
}
else {
phi = (x+1i*y)/sqrt(1-z**2)
long = Im(log(phi)) * 180/pi
}
c(lat,long,ang)
} |
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]/n
y = Q[2]/n
z = Q[3]/n
lat = asin(z) * 180/pi
if (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 quaternion
norme = function(Q){
sqrt(sum(Q**2))
} |
norme = 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é \otimes par la suite
pdt_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)
} |
pdt_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 scalaire
pdt_scal = function(M,N){
return(sum(M*N))
} |
pdt_scal = function(M,N){
return(sum(M*N))
}
un produit vectoriel
pdt_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]))
} |
pdt_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 \mathbf{q}^{-1}
inv_quat = function(Q){
(c(0,0,0,2*Q[4])-Q)/norme(Q)
} |
inv_quat = function(Q){
(c(0,0,0,2*Q[4])-Q)/norme(Q)
}
On peut aussi demander les coordonnées d’un point M' obtenu comme transformation d’un point M par la rotation \mathbf{q}
rotation = function(M, Q){
polaire(pdt_quat(pdt_quat(Q,M),inv_quat(Q)))[1:2]
} |
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, t_1 et t_2. Le fichier a notre disposition contient ainsi des quaternions \mathbf{q}^{P_0}_{t,P} pour une date t et une plaque P (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 P_0. Pour des soucis de calculs, on va supposer qu’on peut interpoler linéairement les quaternions,\mathbf{q}^{P_0}_{t,P}=(1-\lambda)\mathbf{q}^{P_0}_{t_1,P}+\lambda\mathbf{q}^{P_0}_{t_2,P}avec\lambda=\frac{t-t_1}{t_2-t_1}et que récursivement, on peut composer les rotations, au sens ou\mathbf{q}^{P_3}_{t,P_1}\mathbf{q}^{P_3}_{t,P_2}\otimes\mathbf{q}^{P_2}_{t,P_1}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 \mathbf{q}^{P_0}_{t,P} ne sera calculé qu’une fois
projecteur = 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 > 0){
df = rot[rot$Start = t & rot$plate == plates[cur_plate],]
if (dim(df)[1] > 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 > 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)
} |
projecteur = 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 > 0){
df = rot[rot$Start = t & rot$plate == plates[cur_plate],]
if (dim(df)[1] > 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 > 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&rowcount&max_ma=999&min_ma=0"
limit = "" #"&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])
} |
site="http://paleobiodb.org/data1.2/"
req=".txt?datainfo&rowcount&max_ma=999&min_ma=0"
limit = "" #"&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 collection
collection = 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') |
collection = 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’occurrence
occurence = 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') |
occurence = 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] > 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) |
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] > 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 bases
M1 = 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") |
M1 = 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écimens
head(M2,5)
coll_no taxo_no no class family genus
1 5195 55999 373398 Nodosauridae Pawpawsaurus
2 10755 55580 130209 Chaoyangsauridae Chaoyangsaurus
3 10760 38561 144305 Dromaeosauridae
4 10764 66066 130295 Caudipterygidae Caudipteryx
5 10764 66068 130294 Protarchaeopteryx
infraclass kingdom order phylum species
1 Chordata Pawpawsaurus campbelli
2 Chordata Chaoyangsaurus youngi
3 Avetheropoda Chordata
4 Avetheropoda Chordata Caudipteryx zoui
5 Avetheropoda Chordata Protarchaeopteryx robusta
subclass subfamily subgenus suborder subspecies superclass superfamily
1
2
3
4
5
superphylum tribe unranked clade lng lat max_ma min_ma
1 Ornithischia -97.3000 32.86667 105.3 99.60
2 Ornithischia 123.9667 42.93330 150.8 132.90
3 Saurischia 21.0500 46.11667 70.6 66.00
4 Saurischia 120.7333 41.80000 130.0 122.46
5 Saurischia 120.7333 41.80000 130.0 122.46 |
head(M2,5)
coll_no taxo_no no class family genus
1 5195 55999 373398 Nodosauridae Pawpawsaurus
2 10755 55580 130209 Chaoyangsauridae Chaoyangsaurus
3 10760 38561 144305 Dromaeosauridae
4 10764 66066 130295 Caudipterygidae Caudipteryx
5 10764 66068 130294 Protarchaeopteryx
infraclass kingdom order phylum species
1 Chordata Pawpawsaurus campbelli
2 Chordata Chaoyangsaurus youngi
3 Avetheropoda Chordata
4 Avetheropoda Chordata Caudipteryx zoui
5 Avetheropoda Chordata Protarchaeopteryx robusta
subclass subfamily subgenus suborder subspecies superclass superfamily
1
2
3
4
5
superphylum tribe unranked clade lng lat max_ma min_ma
1 Ornithischia -97.3000 32.86667 105.3 99.60
2 Ornithischia 123.9667 42.93330 150.8 132.90
3 Saurischia 21.0500 46.11667 70.6 66.00
4 Saurischia 120.7333 41.80000 130.0 122.46
5 Saurischia 120.7333 41.80000 130.0 122.46
Voila 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) |
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 plaques
download.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) |
download.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 anomalies
x = readLines(rot)
y = gsub( "!101 !", "!", x )
cat(y, file=rot, sep="\n")
remove(x,y) |
x = 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 = NULL
rot_file$anchor2 = NULL
PLATES = sort(unique(rot_file$plate)) |
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 = NULL
rot_file$anchor2 = NULL
PLATES = sort(unique(rot_file$plate))
On ne va garder que les dinosaures qui peuvent etre rattaches à une plaque tectonique
X = M2%>%select(lng, lat)
Y = SpatialPoints(X,proj4string = continents@proj4string)
plaques = over(Y,continents)$PLATEID1
filtre = which(!is.na(plaques))
X0 = X[filtre,]
NBX = dim(X0)[1]
print(paste0(NBX, " spécimens retenus")) |
X = M2%>%select(lng, lat)
Y = SpatialPoints(X,proj4string = continents@proj4string)
plaques = over(Y,continents)$PLATEID1
filtre = 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écimens
PERIOD = M2%>%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"] = plaques |
PERIOD = M2%>%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"] = plaques
On y est presque… on va maintenant remonter de -250 millions d’annees a aujourd’hui, en faisant des bonds de 10 millions d’annees
TMAX = 250
TMIN = 0
PAS = 10 |
TMAX = 250
TMIN = 0
PAS = 10
Pour toutes ces dates, on calcule les quaternions
ROT = 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]) |
ROT = 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 projeter
projete=list()
Xt = list()
X = X0
cpt = 0
for (TIME in seq(TMAX,TMIN,-PAS)){
# setTxtProgressBar(pb, -TIME)
cpt = cpt+1
projete[[cpt]] = continents
Xt[[cpt]] = X0
if (TIME > 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] < -90); inf_180_ = length(inf_180)
sup_180 = which(M[,1] > 90); sup_180_ = length(sup_180)
if (inf_180_ > 0 & sup_180_ > 0) {
if(sup_180_>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] > TIME & PERIOD$min_ma[j] <= 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) > 0){
Xt[[cpt]]=SpatialPointsDataFrame(coords = X[filtre,], data = dataX[filtre,])
} else {
Xt[[cpt]]=X
}
}
} |
projete=list()
Xt = list()
X = X0
cpt = 0
for (TIME in seq(TMAX,TMIN,-PAS)){
# setTxtProgressBar(pb, -TIME)
cpt = cpt+1
projete[[cpt]] = continents
Xt[[cpt]] = X0
if (TIME > 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] < -90); inf_180_ = length(inf_180)
sup_180 = which(M[,1] > 90); sup_180_ = length(sup_180)
if (inf_180_ > 0 & sup_180_ > 0) {
if(sup_180_>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] > TIME & PERIOD$min_ma[j] <= 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) > 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 ere
t = 1 + (TMAX-70)/PAS
leaflet(options = leafletOptions(minZoom = 1)) %>%
addPolygons(data=projete[[t]], weight=2) %>%
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) |
t = 1 + (TMAX-70)/PAS
leaflet(options = leafletOptions(minZoom = 1)) %>%
addPolygons(data=projete[[t]], weight=2) %>%
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 ere
t = 1 + (TMAX-200)/PAS
leaflet(options = leafletOptions(minZoom = 1)) %>%
addPolygons(data=projete[[t]], weight=2) %>%
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) |
t = 1 + (TMAX-200)/PAS
leaflet(options = leafletOptions(minZoom = 1)) %>%
addPolygons(data=projete[[t]], weight=2) %>%
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) |
dinoIcon = makeIcon(iconUrl = "https://www.ludeek.com/wp-content/uploads/2015/03/uploadfsdfsdf1426350179.1426350368774.png",
iconWidth = 30, iconHeight = 50,
iconAnchorX = 15, iconAnchorY = 25)
Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite
