On cherche ici à confirmer le contexte tectonique divergent associé aux dorsales grâce aux données GPS, en exploitant des bases de données de différents types (SQL, CSV, XLSX) avec Python3 et les modules Pandas et Basemap.
Le module sqlite3 est dans la bibliothèque standard de python, vous n’avez donc pas besoin de l’installer.
Les modules suivants ne font pas partie de la bibliothèque standard de python (cliquez sur les liens pour accéder aux sites officiels de ces modules et aux procédures d’installations):
- Pandas permet de faire de l’analyse de données
- Matplotlib et Basemap nous permettront de créer le fond de carte et de le modifier
- Numpy est nécessaire au fonctionnement de matplotlib
- Koala permet de tester facilement des régressions linéaire , exponentielle, et puissance.
Des introductions relatives à ces modules sont disponibles: Pandas, Matplotlib, Numpy
Koala doit être téléchargé ici et placé dans le répertoire actif à partir duquel vous travaillez.
Une fois les modules installés, vous pouvez les importer de la manière suivante:
>>> from sqlite3 import *
>>> import pandas as pd
>>> import numpy as np
>>> from mpl_toolkits.basemap import Basemap
>>> import matplotlib.pyplot as plt
>>> from koala import Koala
Si vous n’avez pas d’UI (interface graphique ... c’est à dire un bureau avec des dossiers , des fichiers, voire même un joli chaton dans un panier fleuri en fond d’écran), ajoutez ces deux lignes au tout début de votre code:
>>> import matplotlib
>>> matplotlib.use('Agg') # no UI backend
Cette base de données au format SQL contient les positions en latitude et longitude des stations EISL, PAMA, et BRAZZ au cours du temps.
Vous pouvez télécharger cette base de données ici.
Connectez vous à la base:
>>> con = connect("gps.db")
Examinez le contenu de la base:
>>> liste_tables = pd.read_sql_query("select name from sqlite_master where type = 'table'", con)
>>> liste_tables
name
0 brazz
1 pama
2 eisl
Récupérez chacune des tables:
>>> brazz = pd.read_sql_query("select * from brazz", con)
>>> eisl = pd.read_sql_query("select * from eisl", con)
>>> pama = pd.read_sql_query("select * from pama", con)
Nous allons utiliser le module Koala afin de reéaliser des régressions linéaires et de récupérer les équations de ces droites de régressions.
Convertissons notre objet DataFrame en objet Koala:
>>> kbrazz = Koala(brazz)
Examinons le contenu de la table kbrazz:
>>> kbrazz.head()
temps lat long
0 1995.1759 0 0
1 1995.1786 0.0992637499621001 1.17157402028476
2 1995.1814 -0.7112473790162 2.24455692153047
3 1995.1841 -0.859168912375599 0.17355308905353
4 1995.1869 -0.6836366309012 0.11077272029
Utilisez ensuite la méthode lin afin de tester une régression linéaire. Cette méthode requiert 2 arguments: les deux séries de données (X et Y) pour lesquelles on souhaite tester la régression:
>>> lat_brazz = kbrazz.lin("temps", "lat")
La méthode lin de l’objet Koala retourne un dictionnaire contenant toutes les informations relatives à la régression:
>>> lat_brazz
AttrDict({'equation': 'y = 1.2228585497565216 x + -2440.3739520157933', 'R': 0.995359092755059, 'a': 1.2228585497565216, 'b': -2440.3739520157933, 'f': 'y = 1.22x +-2440.37', 'r': 'R = 1.0', 'graph': <module 'matplotlib.pyplot' from '/usr/local/lib/python3.6/site-packages/matplotlib/pyplot.py'>})
Si vous possédez une UI, la méthode lin affiche automatiquement le graphique (plot + droite de régression). Dans tous les cas, vous pouvez sauvegarder ce graphique grâce à la commande savefig:
>>> plt.savefig("latbrazz.png")
On peut éventuellement annoter le graphique en y ajoutant par exemple l’équation et le coefficient de corrélation:
>>> plt.text(2000, 0, lat_brazz.f)
Text(2000,0,'y = 1.22x +-2440.37')
>>> plt.text(2000, -1, lat_brazz.r)
Text(2000,-1,'R = 1.0')
>>> plt.savefig("latbrazz2.png")
Faisons de même pour le déplacement de brazz en longitude
>>> lon_brazz = kbrazz.lin("temps", "long")
>>> plt.savefig("lonbrazz.png")
Puis conservons les coefficients directeurs des droites de régressions en les affectant à des variables. Ces coefficients directeurs correspondent à une variation de latitude (ou de longitude) sur une variation de temps, c’est à dire une vitesse de déplacement (en cm/an).
v_lat_brazz = lat_brazz.a
v_lon_brazz = lon_brazz.a
Le vecteur déplacement global est la somme des vecteurs déplacement latitude et longitude. Sa norme est égale à la racine carrée de la somme des carrés des vecteurs déplacement latitude et longitude:
n_brazz = round((v_lon_brazz**2 + v_lat_brazz**2)**0.5,2)
Faisons de même pour EISL et PAMA:
>>> keisl = Koala(eisl)
>>> lat_eisl = keisl.lin("temps", "lat")
>>> plt.savefig("lateisl.png")
>>> lon_eisl = keisl.lin("temps", "long")
>>> plt.savefig("loneisl.png")
>>> v_lat_eisl = lat_eisl.a
>>> v_lon_eisl = lon_eisl.a
>>> n_eisl = round((v_lon_eisl**2 + v_lat_eisl**2)**0.5,2)
>>> kpama = Koala(pama)
>>> lat_pama = kpama.lin("temps", "lat")
>>> plt.savefig("latpama.png")
>>> lon_pama = kpama.lin("temps", "long")
>>> plt.savefig("lonpama.png")
>>> v_lat_pama = lat_pama.a
>>> v_lon_pama = lon_pama.a
>>> n_pama = round((v_lon_pama**2 + v_lat_pama**2)**0.5,2)
Puis conservons toutes ces données dans un dictionnaire; nous en aurons besoin plus tard pour tracer les vecteurs sur le fond de carte.
vecteur = {"brazz":(v_lon_brazz, v_lat_brazz, n_brazz), "eisl":(v_lon_eisl, v_lat_eisl, n_eisl), "pama":(v_lon_pama, v_lat_pama, n_pama)}
On ferme l’objet plt, pour ne pas surcharger le fond de carte que l’on va créer:
plt.close()
On utilise ici le modèle etopo qui permet d’avoir les reliefs.
>>> map = Basemap(projection='aeqd', lon_0 = -85, lat_0 = 10, width = 18000000, height = 15000000)
>>> map.etopo()
>>> map.drawcoastlines()
>>> plt.show()
>>> plt.savefig("fond.png")
Rajoutons les méridiens et les parallèles:
>>> map.drawmeridians(range(0,360,60), dashes=[3,1], linewidth=0.5, labels=[0,0,0,1], zorder=2)
>>> map.drawparallels(range(-40,50,10), dashes=[4,2], linewidth=0.5, labels=[1,0,0,0], zorder=1)
>>> plt.show()
>>> plt.savefig("fondmp.png")
L’argument dashes pemret de gérer la longueur des tirets et l’espacement.
L’argument labels permet de gérer la position des graduations (latitudes et/ou longitudes)
La dorsale medio-atlantique et les fosses sont bien visibles; en revanche la dorsale pacifique n’est pas très nette. Nous allons ajouter les foyers sismiques et les volcans, de manière à mieux distinguer les frontières de plaques.
Le site https://earthquake.usgs.gov donne accès aux séismes récents (1 heure, 1 semaine, ou 1 mois) répertoriés et localisés dans un fichier csv que l’on peut récuperer avec pandas:
>>> seismes = pd.read_csv("https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.csv")
Cependant , selon les mois, l’activité sismique de la dorsale Pacifique n’est pas toujours suffisante pour bien reperer cette dorsale, donc j'utilise cette base dont je ne connais pas la source, mais qui semble tout a fait cohérente (au moins en ce qui concerne la localisation... et c’est tout ce dont nous avons besoin). Il faut la télécharger depuis Gist et l’importer localement.
>>> seismes = pd.read_csv("earthquakes.csv")
Listons les colonnes de cette base de données:
>>> seismes.columns
Index(['time', 'latitude', 'longitude', 'depth', 'mag', 'magType', 'nst',
'gap', 'dmin', 'rms', 'net', 'id', 'updated', 'place', 'type'],
dtype='object')
Et conservons toutes les lignes mais seulement les colonnes latitude et longitude:
>>> df = seismes.loc[:,["latitude", "longitude"]]
Plaçons alors les foyers sismiques sur notre fond de carte:
>>> for i in range(len(df)):
>>> x, y = df.longitude[i], df.latitude[i]
>>> try:
>>> x, y = map(x,y)
>>> map.plot(x, y, marker='*', color='y', markersize=2)
>>> except:
>>> pass
>>>
>>> plt.show()
>>> plt.savefig("seismes.png")
(On aurait éventuellement pu garder la profondeur depth et adapter la boucle for ci dessus afin d’attribuer des couleurs distinctes aux foyers sismiques en fonction de leurs profondeurs afin de "visualiser" le plan de Benioff)
La base de données utilisée ici provient du site data.world. Ce site recense des milliers de bases de données (dans différents formats: sql, csv, xlsx...) qui sont dans l’ensemble très bien documentées. Il est nécessaire de s’inscrire sur le site pour pouvoir y faire des recherches. La base est au format xlsx , et ce lien permet d’y accéder sans être inscrit sur data.world.
>>> volcans = pd.read_excel("https://query.data.world/s/jejujyib5byfslvmbixubz6leb3f4d")
Cette base contient beaucoup de colonnes qui ne nous seront pas utiles, et en plus elle recense les éruptions "significatives" depuis -4360 BC ... donc il va falloir filtrer tout ca afin de ne conserver que les volcans les plus "récents" (actifs depuis 1900) pour ne pas trop surcharger notre carte.
>>> vf = volcans.loc[:,["Year", "Latitude", "Longitude"]]
>>> for v in range(len(vf)):
>>> if vf.Year[v] > 1900:
>>> x, y = vf.Longitude[v], vf.Latitude[v]
>>> try:
>>> x, y = map(x,y)
>>> map.plot(x, y, marker='*', color="r", markersize=2)
>>> except:
>>> pass
>>> plt.show()
>>> plt.savefig("vulcan.png")
>>> stations = [("EISL", -109.57, -27.14), ("PAMA", -149.57, -17.56), ("BRAZZ", -47.80, -15.90)] # liste de tuple (nom, lon, lat)
>>> for s in stations:
>>> x, y = map(s[1], s[2])
>>> map.plot(x, y, marker='+', markersize=10, color="k")
>>> plt.text(x, y, s[0], ha = "left", va = "top", fontsize = 8)
>>> plt.show()
>>> plt.savefig("gpscarte.png")
Il y’a un bug: les arguments ha (horizontal alignement) et va (vertical alignement) sont "inversés", donc pour placer le nom de la station en bas à droite du repère de la station il faut faire comme si on voulait le placer en haut à gauche !!!
On utilise la méthode quiver qui permet de placer des vecteurs. Cette méthode requiert 4 arguments obligatoires:
-
x et y : coordonnées de la base du vecteur
-
u, v : les deux composantes du vecteurs
Les autres arguments sont optionnels, ils permettent notamment de configurer l’échelle et l’apparence du vecteur.
Pour des raisons de place, je n’ai pas utilisé la même échelle pour les différents vecteurs: EISL et BRAZZ sont représentés au 1/8 tandis que PAMA est au 1/15. Attention, l’unité est le pouce (inches) et non pas le centimètre.
>>> x, y = map(stations[0][1], stations[0][2])
>>> map.quiver(x, y, v_lon_eisl, v_lat_eisl, scale_units='inches', scale=8, headwidth=12, headlength=15, units='inches', width=0.01)
>>> x, y = map(stations[1][1], stations[1][2])
>>> map.quiver(x, y, v_lon_pama, v_lat_pama, scale_units='inches', scale=14, headwidth=12, headlength=15, units='inches', width=0.01)
>>> x, y = map(stations[2][1], stations[2][2])
>>> map.quiver(x, y, v_lon_brazz, v_lat_brazz, scale_units='inches', scale=8, headwidth=12, headlength=10, units="inches", width=0.01)
>>> for s in stations:
>>> x, y = map(s[1], s[2])
>>> plt.text(x, y," v = " + str( vecteur[str.lower(s[0])][2]) + " cm/an", ha="left", va="bottom", fontsize=6)
>>> plt.show()
>>> plt.savefig("vecteur.png")
On observe que les stations EISL et PAMA qui sont situées de part et d’autre de la dorsale Pacifique se déplacent dans des directions opposées, ce qui confirme le contexte tectonique divergent associé aux dorsales.