Geologia e Python: come creare una mappa dei terremoti.

Python è un linguaggio all’avanguardia per scrivere script per le Geoscienze. La presenza di moltissime librerie e la sua facilità di apprendimento fanno di Python un linguaggio di programmazione più adatti per la creazione di script nel settore geologico.

In questo post vedremo come creare una mappa dei terremoti usando un database dei terremoti, scaricabile dal sito di INGV (Istituto Nazionale di Geofisica e Vulcanologia), e Basemap.

Dove avvengono i terremoti?

Generalmente i terremoti si concentrano lungo i margini delle placche tettoniche. Le placche tettoniche, durante il loro moto relativo, possono scontrarsi o allontanarsi generando così zone di compressione e zone di distensione. In queste zone si accumulano enormi quantità di energia, che vengono rilasciate quando avviene un terremoto. L’energia rilasciata da un terremoto viene misurata con la magnitudo.

Basemap

Basemap è un toolkit della libreria matplotlib e serve per tracciare dati 2D su mappe in Python. È simile per funzionalità al toolbox per mappe di matlab e alle funzionalità per mappe di IDL, GrADS.

PyNGL e CDAT sono altre librerie che forniscono funzionalità simili in Python.

Basemap di per sé non traccia alcunché, ma fornisce le funzionalità per trasformare le coordinate in uno dei 23 diversi tipi di proiezioni per mappe (usando la libreria C PROJ.4). Viene in seguito usata matplotlib per tracciare contorni, immagini,vettori, linee o punti nelle coordinate trasformate. Sono forniti insiemi di dati per profili costieri, fiumi e confini politici (da Generic Mapping Tools), insieme con metodi per tracciarli. La libreria GEOS viene usata internamente per unire le caratteristiche delle coste e dei confini politici nella regione desiderata per la proiezione su mappa.

Basemap fornisce funzionalità per leggere dati nei formati netCDF e Shapefile, così come direttamente via HTTP usando OpeNDAP. Questa funzionalità viene fornita attraverso il client PyDAP e un’interfaccia Python alla libreria C Shapefile.

Basemap è orientato alle esigenze degli studio sia delle scienze della terra, particolarmente oceanografi e meteorologi. L’autore ha scritto originariamente Basemap per avere un aiuto nelle sue ricerche (previsioni del clima e meteorologiche), dato che al tempo CDAT era l’unico altro strumento in Python per tracciare dati su proiezioni di mappe. Col passare degli anni, le funzionalità di Basemap si sono evolute mano a mano che scienziati di altre discipline (come biologia, geologia e geofisica) hanno richiesto e contribuito con nuove funzionalità.

I terremoti: approfondimenti

Se vuoi approfondire l’argomento trattato in questo post, ti consiglio i seguenti libri:

Pensare in Python

Acquista su Amazon

Capire la Terra. Con aggiornamento online

Acquista su Amazon

Sei digiuno di Python e vuoi acquisire nozioni su questo potente linguaggio di programmazione? Ti consiglio il corso “Python 3.8 Guida Completa: da Principiante a Esperto” (LINK) presente su Udemy. Sicuramente saper usare questo linguaggio di programmazione e le sue librerie è molto utile, anche perchè possiamo creare degli scripts per le nostre analisi invece che comprare software propietari.

Python per visualizzare la distribuzione dei terremoti

Innanzitutto dobbiamo scaricare dal sito dell’INGV il file CSV contenente i dati sui terremoti monitorati dalla reste dell’Istituto Nazionale di Geofisica e Vulcanologia, considerando solo quelli che sono avvenuti dal giorno 1 gennaio 1985 al 27 novembre 2018.

Tra le informazioni presenti nel file CSV al momento ci servono sono le sole coordinate (latitudine e longitudine, espresse nel sistema geografico di riferimento WGS84) dell’epicentro del terremoto.

import pandas as pd

Dobbiamo creare un dataframe contenente le informazioni dei terremoti presenti nel file csv:

df = pd.read_csv("terremoti_ingv.txt", sep="|")

Una volta creato il dataframe selezioniamo le colonne contenenti le coordinate dei sismi :

lat=df.iloc[:,2]  long=df.iloc[:,3]

Selezionate le coordinate creiamo due liste, longitudine e latitudine, che verranno popolate con due cicli for.

longitudine=[] 
latitudine=[]

for i in long:     
    longitudine.append(i) 

for i in lat:     
    latitudine.append(i)

A questo punto possiamo costruire la mappa.

##----costruzione mappa base----##

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
# sono i valori di lat/long degli angoli inferiori sinistro e superiore destro della mappa.
# lat_0 lon_0 sono le coordinate del centro mappa

m = Basemap(llcrnrlon=4.30,
            llcrnrlat=33.24,
            urcrnrlon=28.56,
            urcrnrlat=48.04,
            lat_0=12.49,
            lon_0=41.92,
            resolution='h',
            epsg=4326)

m.drawcoastlines() #disegna le linee di costa
m.drawcountries() #disegna i continenti
m.bluemarble() #disegna i mari/oceani

Quindi possiamo aggiungere alla mappa i meridiani ed i paralleli:

parallels = np.arange(0.,81,10.)
m.drawparallels(parallels)
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians)

Possiamo colorare i continenti, con i loro corsi d’acqua, ed i mari/oceani con questi due righe:

m.fillcontinents(color='gray',lake_color='aqua')
m.drawmapboundary(fill_color='aqua')

Infine possiamo plottare i terremoti e salvare la figura.

###Plott dei terremoti###
x,y = m(longitudine, latitudine)
m.plot(x, y, 'ro', markersize=2)

###Salvataggio della mappa###
plt.savefig('mappa_terremoti.png')

Se sei digiuno di Python e vuoi acquisire nozioni su questo potente linguaggio di programmazione? Ti consiglio questo corso presente su Udemy:

Sicuramente saper usare questo linguaggio di programmazione e le sue librerie è molto utile, anche perché possiamo creare degli scripts per le nostre analisi invece che comprare software proprietari.

Infine ti consiglio questo libro, che tratta l’analisi di dati con Python.

Python for Data Analysis: Data Wrangling With Pandas, Numpy, and Jupyter

Acquista su Amazon