Python e geologia: uno script per la manipolazione e rappresentazione delle giaciture geologiche

Le proiezioni stereografiche sono uno strumento molto importante per noi geologi; infatti attraverso gli “stereonet” possiamo rappresentare le osservazioni che effettiamo in campagna ed effettuare analisi strutturali.
In questo post descriverò uno script che ho sviluppato con Python 3.6 per manipolare e rappresentare delle giaciture geologiche.

Durante un’attività di campagna ho acquisito 45 giaciture geologiche di strato con la notazione dip/dip direction.
In seguito dovendole integrare con altre misure di strato prelevate da dei colleghi, rilevate però con la notazione strike/dip, ho importato in python il file .csv contente le mie misure, solo dopo aver importato nel preambolo dello script il modulo csv:


import csv

with open("./giaciture.csv",newline="",encoding="ISO-8859-1") as file:

lettore = csv.reader(file, delimiter=",")

dip_direction=[]

x=[(riga [1]) for riga in lettore]

for i in x:
i=int(i)

if i<=90:
print(dip_direction.append(i+90))
elif 90<i<=180:
print(dip_direction.append(i-90))
elif 180<i<=270:
print(dip_direction.append(i-90))
elif 270 <i<=360:
print(dip_direction.append(i-270))

print(dip_direction)

</p>

<p>Una volta richiamata la riga che contiene, nel file <em>giaciture.csv</em>, le misure della <em>dip direction</em> sono state convertite con un <strong>ciclo for</strong> nella notazione <em>strike</em> ed inserite nella lista strike=[].
Le giaciture geologiche presenti nelle liste strike=[] e dip=[] (quest’ultima creata riportando la riga dei dip del file giaciture.csv) possono essere rappresentate su degli stereonet grazie al modulo Python <strong>mplstereonet</strong>.

import csv
import matplotlib.pyplot as plt
import numpy as np
import mplstereonet

with open("./giaciture.csv",newline="",encoding="ISO-8859-1") as file:

lettore = csv.reader(file, delimiter=",")

strike=[]

x=[(riga [1]) for riga in lettore]

for i in x:
i=int(i)
if i<=90:
print(strike.append(i+90))
elif 90<i<=180:
print(strike.append(i-90))
elif 180<i<=270:
print(strike.append(i-90))
elif 270 <i<=360:
print(strike.append(i-270))
print(strike)

with open("./giaciture.csv",newline="",encoding="ISO-8859-1") as file:

lettore = csv.reader(file, delimiter=",")

dip=[]

y=[(riga [0]) for riga in lettore]

for i in y:
i=int(i)
if i!=90:
print(dip.append(i+0))
print(dip)

###Proiezione stereografica dei piani###
strike=strike
dip=dip

fig = plt.figure()
ax = fig.add_subplot(111, projection='stereonet')
ax.plane(strike, dip, 'g-', linewidth=0.7)

ax.grid()

plt.show()

### Proiezione Poli ###
fig = plt.figure()
ax = fig.add_subplot(111, projection='stereonet')
ax.pole(strike, dip, 'g^', markersize=5)
ax.grid()
plt.show()

### Diagramma di densità ###
fig = plt.figure()

ax = fig.add_subplot(111, projection='stereonet')

ax.pole(strike, dip, c='k', label='Pole of the Planes')
ax.density_contourf(strike, dip, measurement='poles', cmap='Blues')
#ax.set_title('Density coutour dei poli', y=1.10, fontsize=12)
ax.grid()

plt.show()

stereonet

stereonet

stereonet

Inoltre con il modulo mplstereonet possono essere anche eseguite analisi strutturali; ad esempio si può ricavare l’asse di una piega conoscendo le giaciture dei suoi fianchi.

Dalla geologia strutturale sappiamo che per determinare l’asse di una piega possono essere utilizzati i diagrammi β ed i diagrammi π.
Con i diagrammi β i piani vengono proiettati come grandi cerchi. In una piega cilindrica ogni superficie misurata sui fianchi contiene l’asse della piega: l’intersezione di due di queste superfici definisce una linea, detta asse β, che è parallela all’asse della piega. Se la piega è perfettamente cilindrica, i grandi cerchi rappresentanti le superfici passano tutti per β. Proiettando quindi le superfici misurate sui due fianchi della piega queste si intersecheranno in un unico punto corrispondente all’asse della piega.
Invece con i diagrammi π i piani vengono proiettati mediante i loro poli. Se la piega è cilindrica, i poli dei piani misurati cadono tutti su un grande cerchio (piano π) che rappresenta il piano ortogonale all’asse della piega (in poche parole l’asse della piega corrisponde al polo del piano che contiene tutti i poli delle superfici misurate). Il polo di questo piano, definito dall’asse π, è parallelo all’asse della piega.


import csv
import numpy as np
import mplstereonet
import matplotlib.pyplot as plt

with open("./giaciture.csv",newline="",encoding="ISO-8859-1") as file:

lettore = csv.reader(file, delimiter=",")

strikes=[]

x=[(riga [0]) for riga in lettore]

for i in x:
i=int(i)
if i!=90:
print(strikes.append(i+0))

print(strikes)

with open("./giaciture.csv",newline="",encoding="ISO-8859-1") as file:

lettore = csv.reader(file, delimiter=",")

dips=[]

y=[(riga [1]) for riga in lettore]

for i in y:
i=int(i)
if i!=90:
print(dips.append(i+0))

print(dips)

from collections import OrderedDict

fig = plt.figure()

# Diagramma beta

ax = fig.add_subplot(121, projection='stereonet')
ax.plane(strikes, dips, c='k', label='Piani (Fianchi della piega)')
strike, dip = mplstereonet.fit_girdle(strikes, dips)
ax.pole(strike, dip, c='r', label='Asse beta (Intersezione dei piani)')

# Digramma pi

ax = fig.add_subplot(122, projection='stereonet')
ax.pole(strikes, dips, c='k', label='Poli (Fianchi della piega)')
ax.plane(strike, dip, c='g', label='Fitted GC')
ax.pole(strike, dip, c='r', label='Asse pi (Polo del GC)')

for ax, title in zip(fig.axes[1::2], ['Diagramma beta', 'Diagramma pi']):
ax.set_title(title, y=1.10, fontsize=12)
ax.grid()

# This will avoid repetition in the legend:
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))

ax.legend(by_label.values(), by_label.keys(), loc='upper left')

plt.show()

stereonet

L’uso di script di questo tipo è molto utile, sia per la gestione di dati presenti in file in formato csv che per la costruzione di proiezioni stereografiche senza essere vincolati all’uso di software proprietari.

Potete trovare lo script completo nella sezione Download del sito.

Approfondimenti

Per approfondire l’argomento trattato in questo post si consigliano i seguenti testi:

Ciao, mi chiamo Antonio Nirta e sono un geologo. Mi sono laureato all’Università di Pisa e dal 2017 svolgo la libera professione.

Attraverso i post presenti in questo blog cercherò di trasmettere la mia passione per le Scienze Geologiche e di fornire informazioni che riguardano il bellissimo settore della Geologia.

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *