Aggiornamento bibliografico (che mi ero completamente saltato all'Universita')
Sediment Accumulation Rates and the Completeness of Stratigraphic Sections
Peter M. Sadler
The Journal of Geology, Vol. 89, No. 5 (Sep., 1981), pp. 569-584
Estimation of completeness of stratigraphical sections using
empirical data and theoretical models
empirical data and theoretical models
Peter M.Sadler & David J. Strauss
Journal of the Geological Society, London, Vol. 147, 1990,pp. 471-485
The nature of stratigraphical record (3 ed)
D. Ager (1993)
Wiley and Sons
capitolo 3 More gaps than records pag.43
Fractals models in Earth Sciences
Korvin Gabor (1992)
Elsevier Science Ltd
cap.2 Fractals in Flatland : a romance in <2 dimensions
-----------------------------------------------------------------------------
Seconda puntata dell'esperimento iniziato qui
Su suggerimento di uno degli autori ho provato ad utilizzare i dati di questo articolo in cui sono presenti dettagliate stratigrafie del Scaglia Toscana
Sono state considerati considerati gli affioramenti Sugame (48 dati Middle Eocene-Lower Oligocene) , Castello di Brolio (251 dati Upper Cretaceous-Lower Oligocene) , Cintoia West (150 dati Upper Cretaceous-Middle Upper Eocene ) , Montegrossi (208 dati) e Lucolena North (Upper Cretaceous--Lower Oligocene 63 dati). Come si vede non tutte le serie ricoprono esattamente lo stesso periodo temporale e le stess formazioni. Si evidenzia che alcune serie non presentano una serie deposizionale continua
Il motivo di questo andamento e' perche' i cicli tettonici e climatici influenzano solo le scali maggiori (ma potrebbe essere anche un problema di sotto campionamento degli strati meno potenti)
Questo e' il codice gentilmente fornito da Gemini
import numpy as np
import matplotlib.pyplot as plt
import os
def box_counting_dimension_1d(data, max_box_size=None, plot=True):
"""
Calcola la dimensione frattale di una serie di dati 1D utilizzando il metodo del box-counting.
Args:
data (np.array): La serie di dati 1D (spessore di strato) come array NumPy.
max_box_size (int, optional): La dimensione massima delle scatole da testare.
plot (bool): Se True, visualizza il log-log plot.
Returns:
float: La dimensione frattale calcolata (pendenza del log-log plot).
"""
# 1. Preparazione dei dati (la logica rimane la stessa)
data = np.asarray(data)
# Rimuovi eventuali valori non finiti (NaN, Inf) che potrebbero derivare dal caricamento
data = data[np.isfinite(data)]
if len(data) < 2:
print("Errore: Il set di dati è troppo piccolo o vuoto dopo la pulizia.")
return np.nan
min_val = np.min(data)
max_val = np.max(data)
if max_val == min_val:
return 1.0
L = len(data)
if max_box_size is None:
max_box_size = L // 4
# Raccogli le dimensioni delle scatole (epsilon) da testare
r_values = 2**np.arange(1, int(np.log2(max_box_size)) + 1)
# Scalatura dell'asse Y per mappare i punti in uno spazio 2D approssimativamente quadrato
y_scaled = (data - min_val) / (max_val - min_val) * (L - 1)
points = np.column_stack((np.arange(L), y_scaled))
log_r = []
log_N = []
# 3. Box-Counting
for r in r_values:
# Calcola le coordinate delle scatole
box_coords = np.floor(points / r).astype(int)
# Conta le scatole uniche
N_r = len(np.unique(box_coords, axis=0))
if N_r > 0:
log_r.append(np.log(1/r))
log_N.append(np.log(N_r))
log_r = np.array(log_r)
log_N = np.array(log_N)
if len(log_r) < 2:
print("Errore: Non abbastanza punti per calcolare la pendenza. Prova con una serie di dati più lunga.")
return np.nan
# Nuova sezione: Selezione dei primi 4 punti
num_punti_fit = 3
if len(log_r) < num_punti_fit:
print(f"Attenzione: Trovati solo {len(log_r)} punti. Eseguo la regressione su tutti i punti disponibili.")
log_r_fit = log_r
log_N_fit = log_N
else:
# Seleziona i primi 'num_punti_fit' punti (corrispondenti alle scale più piccole/fini)
#log_r_fit = log_r[-num_punti_fit:]
#log_N_fit = log_N[-num_punti_fit:]
log_r_fit = log_r[:num_punti_fit]
log_N_fit = log_N[:num_punti_fit]
# Esegue la regressione lineare solo sui punti selezionati
coeffs, residuals, rank, singular_values, rcond = np.polyfit(log_r_fit, log_N_fit, 1, full=True)
fractal_dimension = coeffs[0]
# ----------------------------------------------------
# 5. Plot (Opzionale)
if plot:
plt.figure(figsize=(10, 6))
# Plot di TUTTI i punti
plt.plot(log_r, log_N, 'bo', label='Punti Dati Totali')
# Plot della linea di regressione utilizzando solo l'intervallo di punti selezionato
p = np.poly1d(coeffs)
# Plot della linea SOLO sul range dei punti usati per il fit (opzionale, per chiarezza)
plt.plot(log_r_fit, p(log_r_fit), 'r-',
label=f'Fit (Primi {len(log_r_fit)} punti): D = {fractal_dimension:.4f}')
# Evidenzia i punti usati per il fit
plt.plot(log_r_fit, log_N_fit, 'ro', markerfacecolor='none', markersize=10,
label=f'Punti usati per la regressione')
plt.title('Sugame')
plt.xlabel(r'$\log(1/\epsilon)$')
plt.ylabel(r'$\log(N(\epsilon))$')
plt.legend()
plt.grid(True, which="both", linestyle="--")
plt.show()
return fractal_dimension
# --- Esempio di Utilizzo AGGIORNATO ---
NOME_FILE_CSV = "sugame.csv"
print("--- Calcolo della Dimensione Frattale ---")
if not os.path.exists(NOME_FILE_CSV):
print(f"ERRORE: File '{NOME_FILE_CSV}' non trovato.")
print("Assicurati che il file esista nella stessa directory dello script.")
else:
try:
# Carica i dati dal file CSV
# usecols=0 specifica di usare solo la prima colonna
# dtype=int assicura che i dati vengano letti come interi
spessore_strato_dati = np.genfromtxt(
NOME_FILE_CSV,
delimiter=',',
dtype=int,
skip_header=0,
filling_values=np.nan # Sostituisce i valori mancanti con NaN (che puliremo dopo)
)
if spessore_strato_dati.ndim > 1:
# Se per qualche motivo numpy legge più colonne, forziamo solo la prima
spessore_strato_dati = spessore_strato_dati[:, 0]
print(f"Dati caricati con successo da '{NOME_FILE_CSV}'.")
print(f"Numero di punti dati: {len(spessore_strato_dati)}")
# Calcola e stampa il risultato
D = box_counting_dimension_1d(spessore_strato_dati, plot=True)
print(f"\nRisultato:")
print(f"Dimensione Frattale (D) calcolata: {D:.4f}")
except Exception as e:
print(f"Si è verificato un errore durante il caricamento o l'analisi dei dati: {e}")
Ben sapendo che non e' propriamente corretto ho provato a vedere cosa succedeva sommando tutti i dati di tutte le serie
