Informatica per Scienze Geologiche LTa.a.2019-2020

Rappresentazione grafica modelloDTM del Friuli Venezia Giulia

Analisi numerica del DTM

Docente: Prof. Carla Braitenberg, Ricerca dei minimi in collaborazione con Dott. Alberto Pastorutti

Dipartimento Matematica e Geoscienze, Via Weiss 1, Università di TriesteE-mail: berg@units.it Tel. 040 5582258

Modelli der terreno utilizzati• Nella esercitazione utilizzeremo due

database:

• DTM in coordinate geografiche alta risoluzione e bassa risoluzione

• DTM_ortom_reg_regrid10m_WGS84_HD.mat

• DTM_ortom_reg_regrid10m_WGS84_LD.mat

• DTM alta risoluzione in coordinate cartesiane

• DTM_esercizio3.mat

Modello digitale del terreno FVG

Rappresentazione grafica del modello digitale del terreno

• Scopo dell’esercizio:

• Utilizzare il modello digitale del terreno in zona carsica

• - rappresentare il modello digitale del terreno in 3D

• - calcolare le isolinee di quota

%Plot_Topografia3D%caricamento dati digitali del terrenoload DTM_ortom_reg_regrid10m_WGS84_HD.mat;figure%surface plotsurf(X1,Y1,Z1)shading interp;

% visualizzazione 3D (azimuth -78°, elevation 62°)view([-78,62]);

% etichette sugli assixlabel('asse orizzantale X'); ylabel('asse verticale Y');zlabel('asse Z quota');

Grafico della topografia in 3D

Comando View

• view(az,el) and view([az,el]) set the viewing angle for a three-dimensional plot. The azimuth, az, is the horizontal rotation about the z-axis as measured in degrees from the negative y-axis. Positive values indicate counterclockwise rotation of the viewpoint. el is the vertical elevation of the viewpoint in degrees. Positive values of elevation correspond to moving above the object; negative values correspond to moving below the

object.

Grafico creato con lo script Plot_Topografia3D.m

Grafico delle isolinee della topografia • %Plot_DTM_isolinea

• Zmin=200; Zmax=500;

• step=2;

• v=Zmin:step:Zmax;

• sv=size(v);

• load DTM_ortom_reg_regrid10m_WGS84_HD.mat;

• %num2str converte i numeri in stringhe, altrimenti la concatenazione(di stringhe) non e' possibile

• figure('name',['esempio 2D isolinee Zmin: ' num2str(Zmin) ' Zmax: ' num2str(Zmax)])

• %contour rappresenta in grafico le isolinee in 2D

• [C, h]=contour(X1,Y1,Z1,v);

Grafico prodotto con script Plot_DTM_isolinea

Scegli una zona di particolare interesse: aggiungi comando axis. Qui il grafico e’

ristretto su una dolina.• %Plot_DTM_isolinea

• Zmin=180; Zmax=500;

• step=2;

• v=Zmin:step:Zmax;

• sv=size(v);

• load DTM_ortom_reg_regrid10m_WGS84_HD.mat;

•

• %num2str converte i numeri in stringhe, altrimenti la concatenazione(di

• %stringhe) non e' possibile

• figure('name',['esempio 2D isolinee Zmin: ' num2str(Zmin) ' Zmax: ' num2str(Zmax)])

• %contour rappresenta in grafico le isolinee in 2D

• [C, h]=contour(X1,Y1,Z1,v);

• %axis([xmin xmax ymin ymax])

• axis([13.74 13.76 45.715 45.725]);

13.74 13.742 13.744 13.746 13.748 13.75 13.752 13.754 13.756 13.758 13.7645.715

45.716

45.717

45.718

45.719

45.72

45.721

45.722

45.723

45.724

45.725

Segue l’analisi numerica del DTM. Di seguito utilizziamo il DTM proiettato in coordinate

metriche.

• Proiezione utilizzata: RDN2008 fuso 33

• Nome del file: DTM_esercizio3.MAT

Analisi di una superficie:profilo e minimi locali nella topografia

Possiamo trovare automaticamente tutte le doline in un modello del terreno?

A A′

Profilo topografico /1

load('DTM_esercizio3');contourf(X,Y,Z,25);

• Vogliamo estrarre il profilo lungo un segmento qualsiasi

• Il segmento è definito dalle coordinate di inizio e fine

• Otteniamo le coordinate cliccando sulla figura della mappa

• carichiamo il DTM, che contiene Z, X, Y

• contourf: linee di livello con riempimento, 25 livelli

• Questo modello del terreno è in coordinate proiettate, le unità sugli assi sono metri

• È presente un valore di altitudine ogni 10 metri, sia lungo x che lungo y

• Ovvero: la risoluzione è di 10 metri/pixel

Profilo topografico /2

• Vogliamo estrarre il profilo lungo un segmento qualsiasi

• Il segmento è definito dalle coordinate di inizio e fine

• Otteniamo le coordinate cliccando sulla figura della mappa

Scriviamo una funzione per ottenere questo, utilizzeremo:

• ginput per leggere le coordinate dei punti cliccati da una figura aperta

• interp2 per calcolare, interpolando, l’altitudine per ogni punto lungo il profilo– i punti del profilo non coincidono con i nodi della griglia del modello del terreno

Profilo topografico /3

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

• argomenti input:

– X,Y,Z del terreno

– passo con cui campioniamo il profilo

• ovvero: ogni quanti metri un punto

• argomenti output:

– profX, profY coord. dei punti del profilo

– profL vettore della distanza lungo il profilo

• è 0 : passo : lunghezza del profilo

– profZ vettore delle quote del profilo

Profilo topografico /4

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

• argomenti input:

– X,Y,Z del terreno

– passo con cui campioniamo il profilo

• ovvero: ogni quanti metri un punto

• argomenti output:

– profX, profY coord. dei punti del profilo

– profL vettore della distanza lungo il profilo

• è 0 : passo : lunghezza del profilo

– profZ vettore delle quote del profilo

• disegna una mappa con curve di livello riempite, 25 livelli, quindi usa ginput per salvare (2) punti in xclick (coordinate x) e yclick (coordinate y)

• il tracciato del profilo viene poi disegnato sulla mappa, in rosso ‘r’ e con spessore 'LineWidth',4

Profilo topografico /5

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

• adesso costruiamo il tracciato del profilo, cioè:

–il vettore «profL» che descrive la distanza lungo il profilo, da zero alla fine, con un elemento ogni step

–due vettori «profX» e «profY», che sono le coordinate di ogni punto che compone il profilo

x1,y1 = xclick(1),yclick(1)

x2,y2 = xclick(2),yclick(2)

Δx = x2-x1

Δy

= y

2-y 1

lunghezza

Δ𝑥 + Δ𝑦 

azimuth (rispetto all’asse x)

𝑡𝑎𝑛Δ𝑦

Δ𝑥

Profilo topografico /6

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

x1,y1

x2,y2

Δx = x2-x1

Δy

= y

2-y 1

lunghezza

Δ𝑥 + Δ𝑦 

azimuth (rispetto all’asse x)

𝑡𝑎𝑛Δ𝑦

Δ𝑥

profL(1) = 0

profL(end) = lunghezza

quali sono le coordinate di questo i-esimo elemento di profL?

𝑝𝑟𝑜𝑓𝑋 = 𝑝𝑟𝑜𝑓𝐿 cos (𝑎𝑧𝑖)𝑝𝑟𝑜𝑓𝑌 = 𝑝𝑟𝑜𝑓𝐿 sin(𝑎𝑧𝑖)

a cui poi sommiamo le coordinate di x1, y1 (l’inizio del profilo)profX(i)

pro

fY(i)

Profilo topografico /7

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

profL(1) = 0

profL(end) = lunghezza

Osservazione:

Il prodotto * tra profL (un vettore) e il coseno/seno dell’angolo azi (uno scalare) ha come risultato un vettore lungo quanto profL.

Così facendo, per ogni elemento di profL(lunghezza lungo il profilo) abbiamo la sua coordinata X e la sua coordinata Y.

profX(i)

pro

fY(i)

Profilo topografico /8

function [profX,profY,profL,profZ] = ... TracciaProfilo(X,Y,Z,passo)% input dei punticontourf(X,Y,Z,25);hold on[xclick,yclick] = ginput(2);plot(xclick,yclick,'r','LineWidth',4)% costruzione del segmento del profilodeltax = xclick(2)-xclick(1);deltay = yclick(2)-yclick(1);lunghezza = sqrt((deltax)^2 + (deltay)^2);azi = atan2((deltay),(deltax));profL = 0:passo:lunghezza;profX = xclick(1) + (profL * cos(azi));profY = yclick(1) + (profL * sin(azi));% interpolazione delle quote lungo il profiloprofZ = interp2(X,Y,Z,profX,profY);end

123456789

101112131415161718

• Usiamo la funzione interp2, che interpola valori su griglie meshgrid (ad ogni elemento della matrice Z corrisponde un elemento nelle matrici X e Y, con rispettivamente le sue coordinate x ed y)

interp2(X,Y,Z,profX,profY)

meshgridcoordinate

(m x n)

coordinate X

dei punti su cui interpolare

coordinate Y

dei punti su cui interpolare

matrice F(x,y)

(m x n)

nel nostro caso

la quota

Profilo topografico /9

% chiudi finestre aperte da vecchie figure% ed elimina eventuali variabili già esistenticlose allclear variables

% carica la topografia (X,Y,Z)load('DTM_esercizio3');

% chiamata alla funzione TracciaProfilo[profiloX,profiloY,profiloDist,profiloQuota] = ...

TracciaProfilo(X,Y,Z,5);

% figura: mappaAssiMappa = subplot(3,1,[1 2]);contourf(AssiMappa,X,Y,Z,25);hold on;plot(AssiMappa,profiloX,profiloY,'r','LineWidth',4);

% figura: profiloAssiProfilo = subplot(3,1,3);plot(AssiProfilo,profiloDist,profiloQuota,'r');

• Con uno script, come quello qui lato:

1. Carichiamo il modello del terreno

2. Chiamiamo la funzione TracciaProfilo

3. La funzione ci restituisce:

– il tracciato del profilo in coordinate

– i punti lungo il profilo

– la quota in ciascuno dei punti lungo il profilo

4. Creiamo una figura con due subplot su tre righe:

– sopra: la mappa con i contour (occupa due caselle subplot) e il tracciato del profilo

– sotto: il profilo (in x la distanza, in y la quota)

Profilo topografico /10

% chiudi finestre aperte da vecchie figure% ed elimina eventuali variabili già esistenticlose allclear variables

% carica la topografia (X,Y,Z)load('DTM_esercizio3');

% chiamata alla funzione TracciaProfilo[profiloX,profiloY,profiloDist,profiloQuota] = ...

TracciaProfilo(X,Y,Z,5);

% figura: mappaAssiMappa = subplot(3,1,[1 2]);contourf(AssiMappa,X,Y,Z,25);hold on;plot(AssiMappa,profiloX,profiloY,'r','LineWidth',4);

% figura: profiloAssiProfilo = subplot(3,1,3);plot(AssiProfilo,profiloDist,profiloQuota,'r');

Ricerca dei minimi locali /1

Vogliamo individuare i minimi locali su una superficie.

Nello specifico, le doline in un modello digitale del terreno (DTM) in zona carsica.

Le condizioni sono le seguenti:

1. La superficie è molto estesa e i minimi sono numerosi

2. La superficie è campionata regolarmente, in punti discreti(non è una funzione analitica)

3. È presente del rumore, ovvero dei minimi locali che non sono il nostro obiettivo. Esempio: piccole asperità nel terreno.

Ricerca dei minimi locali /2

Vogliamo individuare i minimi locali su una superficie.

Nello specifico, le doline in un modello digitale del terreno (DTM) in zona carsica.

Le condizioni sono le seguenti:

1. La superficie è molto estesa e i minimi sono numerosi

2. La superficie è campionata regolarmente, in punti discreti(non è una funzione analitica)

3. È presente del rumore, ovvero dei minimi locali che non sono il nostro obiettivo. Esempio: piccole asperità nel terreno.

Dobbiamo pensare a una strategia adatta:

• Non è possibile individuarle una ad una, a mano.

• È necessario scegliere un criterio per discriminare tra doline e rumore.

Ricerca dei minimi locali /3Vediamo velocemente una prova con la stessa strategia usata per il livello del pozzo Trebiciano, ora però siamo in 2 dimensioni. Quindi le due condizioni per un minimo locale sono le seguenti:

• Determinante della matrice hessiana > 0

• Derivata seconda > 0 (è sufficiente che 𝑓 > 0)

load('DTM_esercizio3.mat')[gx, gy] = gradient(Z);[gxx, ~] = gradient(gx);[gxy, gyy] = gradient(gy);detHess = gxx.*gyy - gxy.^2;minimi=find(and(detHess>0,gxx>0));minimiX = X(minimi);minimiY = Y(minimi);

con detHess>0 : molto rumore

con detHess>2 : meno rumore

Per la dimostrazione, vedere:http://calvino.polito.it/~taise/Esercizi%202006/2006-05-25%20-%20Appunti%20-%20Hessiana.pdf

matrice Hessiana

Ricerca dei minimi locali /4

contourf(X,Y,Z,50);hold on;scatter(minimiX,minimiY,6,'r');% «6» indica la dimensione dei punti, ‘r’ il colore rosso

Disegniamo i punti trovati su una mappa contourf:

Problemi:

• Non sempre i punti trovati corrispondono al nostro obiettivo (doline)

• Spesso sono presenti più punti per ciascuna dolina (fondo piatto)

P non è una dolina:

entro la finestra c’è un valore inferiore alla quota di P

Ricerca dei minimi locali /5Pensiamo a una strategia differente!

Un punto è una dolina se si verificano entrambe queste condizioni:

1. è il più basso tra tutti i punti entro una finestra rettangolare attorno a sé

2. la differenza tra la media della quota di tutti i punti nella finestra e quella del punto è

superiore a una soglia (da noi impostata)

Abbiamo quindi due parametri, la cui scelta è soggettiva:

• dimensione della finestra attorno al punto

• valore della soglia

P

Q

Q è una dolina:

la quota di Q coincide col minimo entro la finestra

topografia

R

R non è una dolina:

la quota di R coincide col minimo entro la finestra,

ma la differenza rispetto alla quota media entro la finestra è troppo piccola, inferiore a «soglia»

topografia

Ricerca dei minimi locali /6

function [U,V] = CercaDoline(A,finestra,soglia)D = zeros(size(A));for u=1+finestra:size(A,1)-finestra

for v=1+finestra:size(A,2)-finestrafmin = min(min(A(u-finestra:u+finestra,v-finestra:v+finestra)));fmed = mean(mean(A(u-finestra:u+finestra,v-finestra:v+finestra)));if A(u,v)==fmin && fmed-fmin>=soglia

D(u,v) = 1;end

endend[U,V] = find(D==1);end

Scriviamo quanto appena descritto come una funzione.

Argomenti in input: A matrice delle quote, finestra dimensione della finestra, soglia valore della soglia

Argomenti in output: U, V indici degli elementi che sono doline (grazie ad X e Y poi li tradurremo in coordinate)

finestra =

numero di elementi di cui ci allontaniamo dall’elemento (u,v), in ciascuna direzione

Osservazione: il ciclo for evita i bordi della nostra matrice, inizia a una finestra di distanza dal bordo.

Vedi: le due righe del ciclo for.

Ricerca dei minimi locali /7

load('DTM_esercizio3');[U,V] = CercaDoline(Z,5,0.5);dolineX = X(1,V);dolineY = Y(U,1);contour(X,Y,Z,50);hold on;scatter(dolineX,dolineY,6,'r');

Chiamiamo «CercaDoline» in uno script e disegniamo i risultati su una mappa.

finestra = 5, c’è un elemento ogni 10 metri, quindi 5x10x2 = 100 metri

soglia = 0.5 m

estraiamo dalla prima riga di X e dalla prima colonna di Y le coordinate che corrispondono a ciascun punto definito da U e V

• Sembra funzionare meglio.

• Cosa accade variando i due parametri?

Curva ipsometricae istogramma quota doline

for i=1:length(U)dolineZ(i)=Z(U(i),V(i));

endnb = 50; % numero intervalliintervalli = linspace(min(Z(:)),max(Z(:)),nb);subplot(2,1,1);histogram(Z,intervalli,'Normalization','probability');subplot(2,1,2);histogram(dolineZ,intervalli,'Normalization','probability');