Geologia e Python: come creare una mappa dei terremoti.

Abbiamo visto in alcuni post precedenti quanto Python sia adatto alla creazione di script per le Geoscienze.

La presenza di moltissime librerie fa di questo linguaggio di programmazione uno dei più consoni alla creazione di script geologici.

In questo post vedremo come creare una mappa dei terremoti usando un database dei terremoti, prelevato dal sito dell’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.

mappa-placche-tettoniche
Mappa delle placche tettoniche che compongono la crosta terreste. 

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à.

Approfondimenti

Se volete approfondire l’argomento trattato in questo post, vi consiglio i seguenti libri:

Questi testi possono essere acquistati con Amazon Prime, che prevede una prova gratuita di 30 giorni. Potete iscriversi a questo servizio cliccando sul banner sottostante.

Il codice Python per generare la mappa dei terremoti

Innanzitutto dobbiamo scaricare dal sito dell’INGV il file csv contenente i terremoti monitorati dalla reste dell’Istituto Nazionale di Geofisica e Vulcanologia.

In questo script i terremoti che verranno plottati sono quelli che sono avvenuti dall’ 1/1/1985 al 27/11/2018.

Tra le informazioni presenti nel file csv al momento ci servono sono le sole coordinate (latitudine e longitudine, espresse in 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')
La mappa dei terremoti generata dal nostro codice Python

Lascia un commento

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