Utilizza il clustering basato sulla densità e l’analisi di sopravvivenza per stimare quando si verificano i terremoti
introduzione
Sebbene abbiamo una conoscenza abbastanza buona di dove e perché si verificano i terremoti, capire quando accadrà un terremoto è molto impegnativo. In questo articolo utilizzeremo i dati storici sui terremoti e una combinazione di clustering e analisi di sopravvivenza per rispondere alle seguenti domande:
“Come sono distribuiti spazialmente gli eventi sismici?”
“Se si verifica un terremoto, qual è la probabilità che si verifichi un altro terremoto entro la prossima ora?”
“Ci sono differenze regionali nel tempo che intercorre tra gli eventi sismici?”
Ottenere i dati
Utilizzeremo i dati dell’USGS (1) che registra gli eventi sismici (i dati sui terremoti sono di pubblico dominio e vengono forniti per gentile concessione dell’US Geological Survey). Hanno un bel API ciò significa che possiamo ottenere dati sui terremoti in modo semplice in Python:
import http.client
import pandas as pd
import jsonurl = '/fdsnws/event/1/query'
query_params = {
'format': 'geojson',
'starttime': "2020-01-01",
'limit': '10000',
'minmagnitude': 3,
'maxlatitude': '47.009499',
'minlatitude': '32.5295236',
'maxlongitude': '-114.1307816',
'minlongitude': '-124.482003',
}
full_url = f'https://earthquake.usgs.gov{url}?{"&".join(f"{key}={value}" for key, value in query_params.items())}'
print('defined params...')
conn = http.client.HTTPSConnection('earthquake.usgs.gov')
conn.request('GET', full_url)
response = conn.getresponse()
Nella nostra richiesta API abbiamo alcuni parametri:
- limite — Il numero massimo di eventi sismici che desideriamo
- Ora di inizio — Quando dovrebbe verificarsi il primo evento sismico
- magnitudine min — La magnitudo minima di un evento sismico
- latitudine massima — Latitudine massima
- minlatitudine — Latitudine minima
- longitudine massima — Lunghezza massima
- minlongitudine — Lunghezza minima
Imposteremo la magnitudo minima su 3, poiché questa è tipicamente quella che può essere percepita in superficie, e forniremo le coordinate di latitudine e longitudine per creare un riquadro di delimitazione attorno allo stato americano della California. La California è la patria dei famosi Faglia di Sant’Andreaquindi avremo molti eventi sismici da analizzare.
Per gestire la risposta GeoJSON utilizziamo il seguente codice:
if response.status == 200:
print('Got a response.')
data = response.read()
print('made the GET request...')
data = data.decode('utf-8')
json_data = json.loads(data)
features = json_data('features')
df = pd.json_normalize(features)if df.empty:
print('No earthquakes recorded.')
else:
df(('Longitude', 'Latitude', 'Depth')) = df('geometry.coordinates').apply(lambda x: pd.Series(x))
df('datetime') = df('properties.time').apply(lambda x : datetime.datetime.fromtimestamp(x / 1000))
df('datetime') = df('datetime').astype(str)
df.sort_values(by=('datetime'), inplace=True)
else:
print(f"Error: {response.status}")
Ciò restituirà un Pandas DataFrame in cui ogni riga è un evento di terremoto, con colonne che descrivono le proprietà dell’evento (ad esempio, latitudine, longitudine, magnitudo, ID ecc.).
Chi ha una mentalità geospaziale riconoscerà che le coordinate delle chiamate API rappresentano a rettangolo di selezione dei terremoti (cioè trova i terremoti in un’area delimitata), tuttavia ci interessa solo terremoti in California. Pertanto dovremo filtrare tutti i terremoti che non si sono verificati in California (ad esempio, eventi sismici in Nevada e Oregon). Per fare questo utilizzeremo il comodo pacchetto OSMNX per ottenere un file Multipoligono del confine della California:
import osmnx
import geopandas as gpdplace = "California, USA"
gdf = osmnx.geocode_to_gdf(place)
# Get the target geometry
gdf = gdf(("geometry", "bbox_north", "bbox_south", "bbox_east", "bbox_west"))
Successivamente convertiremo le nostre coordinate del terremoto in Geometrie dei punti utilizzando Shapelye poi eseguire a unione spaziale per filtrare i nostri terremoti non californiani:
from shapely.geometry import Point
# Convert to a GeoDataFrame with Point geometry
geometry = (Point(xy) for xy in zip(df('Longitude'), df('Latitude')))
earthquake_gdf = gpd.GeoDataFrame(df, geometry=geometry, crs='EPSG:4326')# Filter to keep only points within the California bounding box
points_within_california = gpd.sjoin(earthquake_gdf, gdf, how='inner', predicate='within')
# Extract latitude, longitude etc.
df = points_within_california(('id', 'Latitude', 'Longitude', 'datetime', 'properties.mag'))
Ora vediamo la nostra prima mappa dei terremoti californiani codificati a colori in base alla magnitudo:
Eccellente, ora abbiamo una mappa dei terremoti della California!
Clustering spaziale dei terremoti
Osservando i terremoti tracciati sulla mappa in precedenza puoi vedere che gli eventi sono disposti spazialmente in un orientamento lineare SE-NW, dove puoi vedere circa 2-3 chiari raggruppamenti di terremoti a forma lineare. Ciò ha senso perché:
- I terremoti si verificano a causa dei movimenti lungo le faglie, che sono caratteristiche lineari (cioè crepe nella crosta terrestre).
- L’orientamento degli eventi sismici corrisponde all’orientamento della zona della faglia di Sant’Andrea.
Il nostro obiettivo è ora raggruppare questi eventi sismici in base alla loro posizione, in modo da poter produrre cluster spaziali associati alle faglie regionali.
Per il clustering utilizzeremo HDBSCAN (Hierarchical Density-Based Spatial Clustering of Applications with Noise), un algoritmo di clustering basato sulla densità utile per determinati tipi di set di dati, come quelli con cluster di forma irregolare e densità di cluster variabili. Ciò è rilevante per il nostro set di dati poiché possono verificarsi eventi sismici grappoli di forma irregolare (cioè lungo faglie di diverso orientamento) e in densità diverse (vale a dire, alcune zone di faglia sono più attive di altre). Anche HDBSCAN lo è robusti ai valori anomaliin questo caso possono trattarsi di terremoti che si verificano lontano dalle zone di faglia.
Attraverso alcuni esperimenti i seguenti parametri producono risultati accettabili (sono coinvolti un po’ di tentativi ed errori, ma guidati dalle nostre ipotesi precedenti su come i terremoti dovrebbero essere raggruppati):
# Fit HDBSCAN
clusterer = HDBSCAN(min_cluster_size=200, metric='haversine', min_samples=20, cluster_selection_epsilon=0.05)
result_df('cluster_label_hdbscan') = clusterer.fit_predict(data_scaled)# Find the number of unique cluster labels
num_clusters = result_df('cluster_label_hdbscan').nunique()
# Create a list of colors
colors = ('lightgray', 'red', 'blue', 'green')(:num_clusters)
# Create a Folium map centered at the mean coordinates
map_center = (result_df('Latitude').mean(), result_df('Longitude').mean())
mymap = folium.Map(location=map_center, zoom_start=6)
# Add markers for each data point with cluster color
for _, row in result_df.iterrows():
cluster_color = colors(row('cluster_label_hdbscan') + 1) # Map cluster label to color
folium.CircleMarker(
location=(row('Latitude'), row('Longitude')),
radius=2,
color=cluster_color,
fill=True,
fill_color=cluster_color,
fill_opacity=0.7,
popup=f'Cluster: {row("cluster_label_hdbscan")}'
).add_to(mymap)
mymap
Eccellente, ora abbiamo utilizzato HDBSCAN per raggruppare i nostri eventi sismici in tre aree, che corrispondono alle nostre precedenti conoscenze regionali. Si noti che alcuni terremoti sono considerati valori anomali e colorati in grigio.
Esaminando il tempo tra i terremoti
Poiché abbiamo la data e l’ora di ogni terremoto, possiamo calcolare il tempo tra ogni evento (ovvero, la durata temporale tra due terremoti), necessaria per la nostra analisi di sopravvivenza nella sezione successiva. Abbiamo già ordinato i nostri valori in base a ‘appuntamento’ colonna in modo che siano in ordine cronologico, ora dobbiamo calcolare la differenza oraria:
from datetime import datetime# Convert 'time' column to datetime objects
df('time') = pd.to_datetime(df('datetime'))
# Sort dataframe by time
df.sort_values(by=('time'), inplace=True)
# Calculate time elapsed between consecutive events
df('time_elapsed') = df('time').diff().shift(-1)
Dato che non ci sarà alcuna differenza di orario per l’ultimo terremoto (cioè il più recente), possiamo calcolare la differenza a partire da un orario che assumeremo come quello attuale:
# Assuming present time is '2024-01-08 16:06:00.000'
present_time = pd.to_datetime('2024-01-08 16:06:00.000')# For the last event of each group, replace NaN with time between event and present time
df('time_elapsed').fillna(present_time - df('time'), inplace=True)
Infine aggiungeremo una nuova colonna che indicherà se un terremoto è stato seguito da un altro terremoto. Anche se questo sembra un po’ ridondante, dato che solo l’ultimo terremoto ne è interessato, evidenzia un punto interessante sull’uso dell’analisi di sopravvivenza per questi tipi di compiti (di cui parleremo più avanti):
# Label events based on whether they happened or not
df('event_happened') = df('time_elapsed').apply(lambda x: 1 if pd.Timedelta(days=0) > 0 else 0)
Ecco un istogramma del tempo trascorso tra gli eventi sismici:
Possiamo vedere che la maggior parte dei terremoti si verificano in breve successione l’uno dall’altro (circa 2,5 ore l’uno dall’altro), con molti meno terremoti che si verificano a intervalli di tempo più lunghi (ovvero, il tempo massimo tra i terremoti in questo set di dati è di 350 ore).
Analisi di sopravvivenza degli intervalli di tempo dei terremoti
Abbiamo i nostri cluster di terremoti, così come il tempo che intercorre tra i terremoti, quindi ora utilizziamo l’analisi di sopravvivenza per informarci sulla probabilità che si verifichino terremoti. I nostri obiettivi saranno quindi:
- Utilizzare l’analisi di sopravvivenza per creare curve che mostrino la probabilità che si verifichi un terremoto nel tempo, dato che il terremoto si è appena verificato.
- Confronta le nostre curve di probabilità nelle nostre regioni sismiche già raggruppate, per vedere se ci sono somiglianze/differenze nell’attività sismica regionale.
Cos’è l’analisi di sopravvivenza?
Ci sono molti articoli ben scritti (2,3,4) che copre questo aspetto, ma in breve, l’analisi di sopravvivenza è una tecnica statistica utilizzata per analizzare i dati relativi al tempo trascorso all’evento, tipicamente nel contesto di studi biomedici o osservazionali. Si concentra sulla stima e sulla modellazione del tempo fino al verificarsi di un evento di interesse, come la morte, il fallimento o un altro risultato specifico, fornendo informazioni sulla probabilità che gli eventi si verifichino nel tempo e sull’impatto delle covariate sul verificarsi dell’evento.
Nel nostro caso si è verificato un terremoto e vogliamo conoscere la probabilità che si verifichi un altro terremoto nel tempo. Nota: presupponiamo che gli eventi sismici siano indipendenti, nonché che recuperiamo solo eventi con una magnitudo minima pari a 3, il che rappresenta un limite di questo lavoro poiché semplifica eccessivamente la dinamica dei terremoti.
IL Stimatore di Kaplan-Meier è un metodo non parametrico utilizzato nell’analisi di sopravvivenza per stimare la funzione di sopravvivenza, che ci fornirà le nostre probabilità. Utilizzeremo il pacchetto “linee di sicurezza” per adattare lo stimatore di Kaplan-Meier ai terremoti in ciascun cluster, per produrre una curva di probabilità:
import plotly.graph_objects as go
from lifelines import KaplanMeierFitter# Initialize Kaplan-Meier Fitter
kmf = KaplanMeierFitter()
# Create a Plotly Figure
fig = go.Figure()
color_map = {0: 'red', 1: 'blue', 2: 'green', -1: 'gray'}
# Fit and plot survival curves for each cluster excluding cluster -1
for label, group in result_df.groupby('cluster_label_hdbscan'):
if label != -1:
durations = group('time_elapsed').dt.total_seconds() / 3600 # Convert hours to days
event_observed = group('event_happened').values
kmf.fit(durations=durations, event_observed=a)
# Plot survival function
fig.add_trace(go.Scatter(
x=kmf.survival_function_.index,
y=1-kmf.survival_function_.KM_estimate,
mode='lines',
name=f'Cluster {label}',
line=dict(color=color_map(label))
))
# Customize the plot
fig.update_layout(
xaxis_title='Elapsed Time (hours)',
yaxis_title='Probability of experiencing an earthquake',
legend_title='Cluster Labels',
width=800, # Adjust width
height=500, # Adjust height
margin=dict(l=50, r=20, t=50, b=50), # Adjust margins
legend=dict(x=0.8, y=0.99), # Adjust legend position
)
# Show the plot
fig.show()
Prima di vedere i risultati, nota il parametro “event_observed” che utilizza la colonna “event_happened” creata in precedenza. Il nostro terremoto finale è un esempio di punto dati censurato a destra, poiché non si è verificato alcun terremoto successivo, ed è quindi trattato come un’osservazione parziale.
Ora per visualizzare le curve:
Ogni gruppo di terremoti ha la propria curva colorata, dove l’asse X rappresenta il tempo trascorso dopo un terremoto e l’asse Y rappresenta la probabilità che si verifichi un terremoto. Alcune osservazioni interessanti:
- La forma delle curve è simile.
- La pendenza iniziale delle curve varia in base alla posizione del cluster.
Concentriamoci sul momento iniziale successivo al verificarsi di un terremoto, poiché la maggior parte dei nostri dati mostra che i terremoti si verificano in una successione relativamente rapida:
Ciò suggerisce che il il tempo tra i terremoti non è lo stesso a livello regionale. Osservando il grafico possiamo dire che una volta che si è verificato un evento sismico, la probabilità che si verifichi un altro terremoto entro dieci ore è di circa il 35% per il cluster 0, il 60% per il cluster 1 e il 48% per il cluster 2. Ciò dimostra che i terremoti in la regione che appartiene al cluster 1 ha avuto un rapido successo rispetto alle altre due aree.
Possiamo ricontrollarlo osservando un istogramma del tempo trascorso tra i terremoti per ciascuna posizione dell’ammasso:
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import numpy as np# Filter data
filtered_df = result_df(result_df('cluster_label_hdbscan') != -1)
# Create subplots in a 2 by 2 grid
fig = make_subplots(rows=2, cols=2, subplot_titles=(f'Cluster {label}' for label in filtered_df('cluster_label_hdbscan').unique()),
shared_xaxes='all', shared_yaxes='all')
# Create histograms for elapsed time for each cluster
for i, (label, group) in enumerate(filtered_df.groupby('cluster_label_hdbscan')):
# Calculate histogram
hist_data, bin_edges = np.histogram(group('time_elapsed').dt.total_seconds() / 3600, bins=20)
# Add histogram trace to subplot
fig.add_trace(go.Bar(
x=bin_edges,
y=hist_data,
opacity=0.5,
name=f'Cluster {label}',
marker=dict(color=color_map(label))
), row=(i // 2) + 1, col=(i % 2) + 1)
# Customise subplot
fig.update_xaxes(title_text='Elapsed Time (hours)', row=(i // 2) + 1, col=(i % 2) + 1)
fig.update_yaxes(title_text='Frequency', row=(i // 2) + 1, col=(i % 2) + 1)
# Update layout
fig.update_layout(
showlegend=False,
margin=dict(l=50, r=50, t=50, b=50), # Adjust margins
height=600,
width=800,
)
# Show the plot
fig.show()
Come mostrato dalle curve di probabilità, il cluster 1 ha intervalli più lunghi tra gli eventi sismici rispetto ai cluster 2 e 3.
Conclusione
Questo lavoro è stato utilizzato raggruppamento spaziale E analisi di sopravvivenza per svelare modelli temporali e geografici all’interno dei dati sui terremoti. Utilizzando tecniche come HDBSCAN per il clustering spaziale e il Kaplan-Meier stimatore per l’analisi della sopravvivenza, abbiamo acquisito preziose informazioni sul variazioni regionali nel tempo tra i terremotiche è stato tradotto in curve di probabilità che possono essere utili per la valutazione del rischio e la preparazione nelle regioni a rischio sismico.
Grazie per aver letto!
Dichiarazione di non responsabilità: questo articolo deve essere utilizzato solo a scopo didattico.
Fonte: towardsdatascience.com
