Frattalita' in misure di spessore (2)

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  

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

The Scaglia Toscana Formation of the Monti del Chianti:
new lithostratigraphic and biostratigraphic data
Enrico Pandeli (*, **), Milvio Fazzuoli (*), Fabio Sandrelli (***), Roberto Mazzei (***),
Simonetta Monechi (*), Marisa Nocchi (****), Ivan Martini (***) & Gigliola Valleri (*) Ital. J. Geosci., Vol. 137, No. 1 (2018), pp. 38-61, 12 figs., 1 tab., 1 plate, 2 extrapl. (doi: 10.3301/IJG.2017.16)
© Società Geologica Italiana, Roma 2018

 

 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

In generale tutti i grafici mostrano una curva spezzata con valori di D>1 sulle scale di maggiori dimensioni e valori nettamente inferiori a 1 per le scale di dettaglio. Alle scale maggiori emerge quindi una tendenza alla frattalita' ed autosomiglianza che si perde alle piccola scala

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

 

 

 

Esponente di Hurst su torbidita' di acquiferi apuani

In questo post vengono analizzati i valori di torbidita' della sorgente del Frigido (periodo 13/02/2020 01/12/2025 dati orari ) e Carbonera (26/07/2006 21/06/2013 dati ogni 3 minuti ) legate all'acquiferi carsico delle Alpi Apuane

 

Torbidita' del tempo Frigido

 

Il parametro Esponente di Hurst e' stato introdotto Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116(1), 770-808. 

E' una misura dell'inerzia del sistema idrogeologico (ovvero capacita' di stoccaggio) e della complessita' della rete carsica. Per i valore di 0.82 e 0.85 trovati si ha che il sistema non e' casuale (H=0.5), non e' deterministicamente semplice ed il comportamento e' correlato su differenti scale di tempo ed e' un sistema che ha "memoria" del proprio passato

Non e' possibile prevedere in maniera completamente deterministica il valore esatto del parametro indagato  

Frigido

Con modelli deterministici si avra' un discreto successo su scale di tempo brevi ma non sulla previsione a lunga scadenza. Il passaggio da regime deterministico a frattale e' indicato come tempo di crossover e si individua dai punti di flesso sull'andamento dei punti. Per Carbonera (grafico sottostante) il punto di cross over e'  tra 1000 e 10000 unita' che tradotto nel passo di campionamento si arriva ad intervallo da le 5 e le 50 ore

Carbonera

Valori di H vicini ad 1 indicano persistenza nei valori (sia che siano massimi che minimi) e questo puo' incidere anche sul valore massimo e non solo sulle code della distribuzione

La memoria non si esaurisce mai ma ha una legge di decadimento come una potenze (in questo caso 2*0.84 -2 = -0.32 

 

con questa legge e l'esponente calcolate la correlazione ha un peso superiore al  10% fino a 1000 intervalli temporali (rho = 0.10^(1/-0.32)

Il fatto che il calcolo di H risulti sostanzialmente identico su serie con diverso passo di campionamento e' un indice della frattalita' del fenomeno osservato (e' indipendente dalla scala del tempo di osservazione)

 

il valore di Hurst prossimo ad 1 indica che e' probabile che la torbidita' aumenti nel prossimo lag (passo di campionamento) e' fortemente dipendente da cio' che accade nel lag precedente . Cio' sarebbe in contrasto con il concetto di massimo relativo della curva di scarica ma la statistica deve essere considerata valida al mantenersi delle condizioni al contorno e delle relative forzanti. Con un valore di 0.82 la coda di esaurimento sara' lunga

 

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

# --- 1. FUNZIONE DEFINIZIONE per l'ANALISI R/S MANUALE ---

def compute_rs_data(series):

    """

    Calcola la statistica Rango Riscalato (R/S) per una serie temporale.

    Restituisce gli intervalli temporali (n) e i valori medi di R/S.

    """

    N = len(series)

    # Definisce gli intervalli temporali (n) in progressione geometrica fino a N/4

    n_values = np.unique(np.logspace(1, np.log10(N // 4), num=50, dtype=int))

    n_values = n_values[n_values >= 1]

   

    rs_values = []

   

    for n in n_values:

        if n == 0: continue

       

        num_blocks = N // n

        rs_per_block = []

       

        for i in range(num_blocks):

            block = series[i * n : (i + 1) * n]

           

            mean = np.mean(block)

            std_dev = np.std(block)

           

            if std_dev == 0:

                continue

            X = np.cumsum(block - mean)

            R = np.max(X) - np.min(X)

           

            rs_per_block.append(R / std_dev)

       

        if rs_per_block:

            rs_values.append(np.mean(rs_per_block))

        else:

            break

    return n_values[:len(rs_values)], np.array(rs_values)

# --- 2. CARICAMENTO E PREPARAZIONE DEI DATI ADATTATI ---

FILE_NAME = "frigido.csv"

COLUMN_NAME = "Torbidita NTU"

FILTER_N_THRESHOLD = 100

try:

    # ATTENZIONE: Caricamento con separatore PUNTO E VIRGOLA (;) e decimale VIRGOLA (,)

    df = pd.read_csv(FILE_NAME, sep=';', decimal=',')

   

    # Seleziona la colonna 'Torbidita NTU'

    discharge_series = df[COLUMN_NAME]

   

    # Conversione in numerico e gestione dei NaN (valori mancanti o non validi)

    discharge_series = pd.to_numeric(discharge_series, errors='coerce').dropna()

   

    # Converti in array NumPy

    discharge_data = discharge_series.values

   

except FileNotFoundError:

    print(f"Errore: Il file '{FILE_NAME}' non è stato trovato. Assicurati che sia nella directory corretta.")

    exit()

except KeyError:

    print(f"Errore: La colonna '{COLUMN_NAME}' non è stata trovata nel file. Controlla il nome esatto.")

    print(f"Colonne disponibili: {df.columns.tolist()}")

    exit()

if len(discharge_data) < 100:

     print("Errore: Non ci sono abbastanza punti dati validi per eseguire l'analisi R/S.")

     exit()

# --- 3. CALCOLO R/S ---

rs_n, rs_rs = compute_rs_data(discharge_data)

# Calcolo di H (Esponente di Hurst) per l'INTERVALLO COMPLETO

log_n_full = np.log10(rs_n)

log_rs_full = np.log10(rs_rs)

coefficients_full = np.polyfit(log_n_full, log_rs_full, 1)

H_full = coefficients_full[0]

c_full = 10**coefficients_full[1]

# --- 4. PLOT: ANALISI R/S INTERVALLO COMPLETO ---

plt.figure(figsize=(8, 5))

plt.plot(rs_n, c_full * rs_n**H_full, color="deepskyblue", label=f'Pendenza Frattale (H={H_full:.2f})')

plt.scatter(rs_n, rs_rs, color="navy", alpha=0.5, label='Dati R/S empirici')

plt.xscale('log')

plt.yscale('log')

plt.xlabel('Intervallo Temporale (n)')

plt.ylabel('Rango Riscalato R/S')

plt.title(f'Scaling Frattale della Torbidità (Full Range)')

plt.legend()

plt.savefig('Frigido_Full_R_S_Plot.png')

# --- 5. PLOT: INTERVALLO FILTRATO (n > 100) ---

filter_mask = rs_n > FILTER_N_THRESHOLD

X_filtered = rs_n[filter_mask]

Y_filtered = rs_rs[filter_mask]

# Calcolo di H per il solo intervallo filtrato

if len(X_filtered) > 1:

    log_n_filtered = np.log10(X_filtered)

    log_rs_filtered = np.log10(Y_filtered)

   

    coefficients_filtered = np.polyfit(log_n_filtered, log_rs_filtered, 1)

    H_filtered = coefficients_filtered[0]

    c_filtered = 10**coefficients_filtered[1]

   

    plt.figure(figsize=(8, 5))

    plt.scatter(X_filtered, Y_filtered, color="darkgreen", alpha=0.7, label=f'Dati R/S empirici (n > {FILTER_N_THRESHOLD})')

    plt.plot(X_filtered, c_filtered * X_filtered**H_filtered, color="red", linestyle='--',

             label=f'Pendenza Frattale Filtrata (H={H_filtered:.2f})')

    plt.xscale('log')

    plt.yscale('log')

    plt.xlabel('Intervallo Temporale (n)')

    plt.ylabel('Rango Riscalato R/S')

    plt.title(f'Scaling Frattale della Torbidità (Intervallo Temporale > {FILTER_N_THRESHOLD})')

    plt.legend()

    plt.savefig('Frigido_Filtered_R_S_Plot_gt_100.png')

   

    print("\n--- Risultati del Calcolo ---")

    print(f"Hurst Esponent (Torbidità - Intervallo Completo): {H_full:.4f}")

    print(f"Hurst Esponent (Torbidità - Filtrato n>{FILTER_N_THRESHOLD}): {H_filtered:.4f}")

else:

    print(f"Nota: Non ci sono abbastanza punti dati rimasti dopo il filtraggio (n > {FILTER_N_THRESHOLD}) per generare il grafico filtrato.")

 

 

Frattalita' in spessore di strato

Aggiormento: un tentativo piu' serio e' qui 

Sempre leggendo i libri dei post precedenti ho visto che e' stata trovata una correlazione frattale nello spessore degli strati turbiditici arenacei. Avendo fatto la tesi sulle Arenarie del Cervarola mi sono incuriosito ma ho misure di sezioni o affioramenti

Gli unici dati che ho trovato in merito a spessori di strato sono relativi agli stratotipo ed in particolare dell'affioramento di Punta Piccola per il Piacenziano (The Global Standard Stratotype-section and Point (GSSP) of the Piacenzian Stage (Middle Pliocene) by D. Castradori, D. Rio, F. J. Hilgen, and L. J. Lourens https://doi.org/10.18814/epiiugs/1998/v21i2/003.Non si tratta di una formazione turbiditica ma non avendo altri dati ci ho provato 

Questa e' la distribuzione di frequenza degli spessori..un limite e' che la distribuzione non e' oltre un ordine di grandezza 

 

 plottando log-log lo spessore degli strati con il numero di strati con spessore superiore alla soglia si ha che

 

usando lo spettro di potenza 

il valore β≈0.47 e' compreso tra il rumore bianco ed ed il fattore di scala 1/f dove β=1 (Rothman, D. H., Grotzinger, J. P., & Richter, F. M. (1994). Scaling in the Sedimentary Rock Record. Nature, 372(6505), 337–342. suggerisce che il valore atteso per le sequenze stratigrafiche e' β=1 altrimenti detto Pink Noise )

Cerchero' dati piu' accurati