Informatica per Scienze Geologiche LTa.a.2017-2018
Esercitazione Analisi dati PozzoTrebiciano
Parte B2Parte B2
Slides preparate da Dott. Pastorutti
Docente: Prof. Carla Braitenberg, Dipartimento Matematica e Geoscienze, Via Weiss 1, Università di Trieste
E-mail: [email protected] Tel. 040 5582258
Esercizio livello acqua nell’Abisso di Trebiciano
• Anni di osservazione: 2006-2008, un campione ogni ora. Unita’ di misura: m sopra il livello del mare.
• Caricare il file. Nome dell’array: liv_Trebiciano
load('livello_AbissoTrebiciano_Tdec.mat');
• Determinare valore medio, valore massimo e data del valore massimo valore massimo
• Comandi necessari:
• max(livello_tdec): cerca il massimo sulla quinta colonna, comprendente il livello.
• mean(livello_tdec)
• std(livello_tdec)
%Max_TrebicianoWater
%read the data file and find extreme value.
load('livello_AbissoTrebiciano_Tdec.mat');
liv=livello_tdec;
% contains the workspace liv_Trebiciano with year, month, day hour, value
n=length(liv(:,1));
% extreme value
[xmax, kmax] = max(liv(:,6));[xmax, kmax] = max(liv(:,6));
tmax=liv(kmax,2:4);
disp([' max value at ' ,num2str(tmax),' value= ', num2str(liv(kmax,6)),'m']);
media=mean(liv(:,6)); scarto= std(liv(:,6));
disp(['average ' ,num2str(media),'m',' std= ', num2str(scarto),'m' ] );
Continua il programma: qui la parte che crea il grafico
figure;
plot(liv(:,1),liv(:,6))
title('Livello acqua nel pozzo Abisso Trebiciano');
xlabel('tempo in anni');
ylabel('livello (m)');ylabel('livello (m)');
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
%% valori estremiload('livello_AbissoTrebiciano_Tdec_CORRETTO.mat');liv=livello_tdec;n=length(liv(:,1));[xmax, kmax] = max(liv(:,6));tmax=liv(kmax,2:4);disp(['max value at ' ,num2str(tmax),...
• Lo script fatto finora: • Risultato:
Parte B/2
' value= ', num2str(liv(kmax,6)),'m']);media = mean(liv(:,6));scarto = std(liv(:,6));disp(['average ' ,num2str(media),'m',...
' std= ', num2str(scarto),'m' ] );
%% graficofigura1 = figure;assi1 = axes('Parent',figura1);plot(assi1,liv(:,1),liv(:,6))hold ontitle('Livello acqua nel pozzo Abisso Trebiciano');xlabel('tempo in anni');ylabel('livello (m)');
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
...[xmax, kmax] = max(liv(:,6));...
• Massimo assoluto: possiamo usare xmax, kmax già calcolati prima
scatter(assi1,liv(kmax,1),xmax)
x y
NB: la figura ottenuta prima non va chiusa,
altrimenti:
x y
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
function [M,I] = maxsort(in,n)M = zeros(1,n); I = M;for i=1:n
[M(i),I(i)] = max(in);in(I(i)) = 0;
endend
• cercare i primi n valori più grandi: definiamo una funzione
1234567 end
• Salvate la function in un file con lo stesso nome della chiamata, quindi maxsort.m
«Ottieni il valore M e l’indice I degli n valori più grandi di in.»
Nella riga 2 vengono creati M ed I, vettori lunghi n, composti di soli zeri.
Le righe da 3 a 6 appartengono a un ciclo for: ad ogni iterazione viene calcolato l’n-esimo massimo e il suo indice, quindi quell’elemento di in viene sostituito con zero.
7
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
[massimi, indici] = maxsort(liv(:,6),100);scatter(assi1,liv(indici,1),massimi);
• Usiamo la nuova funzione «maxsort» nel nostro script
i primi 100 valori
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
[massimi, indici] = maxsort(liv(:,6),1000);scatter(assi1,liv(indici,1),massimi);
• Usiamo la nuova funzione «maxsort» nel nostro script
• Non è una strategia efficace.• Non è una strategia efficace.
• Potremmo calcolare i massimi in finestre di tempo, ma e’ una strategia soggettiva: devo definire inizio e lunghezza della finestra.
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
• Un altro metodo: calcoliamo le derivate prima e seconda
• Cerchiamo i punti in cui
– la derivata prima è zero (entro una tolleranza)
– la derivata seconda è negativa (concavità verso l’alto)
• Produciamo tre grafici sovrapposti, usando subplot
Calcoliamo la derivata (numerica) con diff:Calcoliamo la derivata (numerica) con diff:
in questa maniera usiamo il passo corretto tra due campioni, che non è costante
step_livelli = liv(2:end,1) - liv(1:end-1,1);der1 = diff(liv(:,6)) ./ step_livelli;der2 = diff(der1) ./ step_livelli(1:end-1);der1 = der1(1:end-1);
dopo aver calcolato diff, il risultato è più corto di un elemento. Rendiamo la derivata 1° lunga quanto la derivata 2°, togliendo l’ultimo elemento
(fermandoci a end-1)
Dal help di matlab: diff(X), for a vector X, is [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)].
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
• Approssimiamo la derivata prima «uguale a zero» entro una tolleranza (più o meno tol_der1)
• Approssimiamo la derivata seconda «negativa» entro una tolleranza (minore di tol_der2)
• Teniamo solo i massimi locali trovati che sono superiori a tol_liv
tol_der1 = 1.25e3;tol_der2 = -0.001;tol_liv = 20;massimi_locali = find(and(and(and(der1<=tol_der1,der1>=-tol_der1),...
der2<=tol_der2),...der2<=tol_der2),...liv(1:end-2,6)>=tol_liv));
Ricerca dei massimi localinella serie livello_AbissoTrebiciano
• Tre grafici, con subplot (livello, prima derivata, seconda derivata)
figura2 = figure;AssiLivelli = subplot(4,1,[1 2]);plot(AssiLivelli,liv(:,1),liv(:,6))hold onscatter(AssiLivelli,liv(massimi_locali,1),liv(massimi_locali,6));
AssiDer1 = subplot(4,1,3);AssiDer1 = subplot(4,1,3);plot(AssiDer1,liv(1:end-2,1),der1)
AssiDer2 = subplot(4,1,4);plot(AssiDer2,liv(1:end-2,1),der2)
linkaxes([AssiLivelli,AssiDer1,AssiDer2],'x');