Tutorial esercitazioni inquinamento elettromagnetico a bassa frequenza
	Aggiornato al 23 marzo 2015
	http://embiolab.ifac.cnr.it/~zoppetti/MasterGIS/TutorialMasterGIS/

Indice

• Introduzione
• Link utili
• Impostazioni preliminari
  • Generazione proiezione GAUSS-BOAGA
  • Proprietà del progetto
  • Generazione nuovo Database / Location / Mapset
• Visualizzazione Raster, vista panoramica e bookmarks
• Modello del terreno
  • Impostazione della regione di lavoro di GRASS
  • Importazione shapefile punti quotati e curve di livello.
  • Trasformazione degli shapefile in in layer vettoriali GRASS.
  • Trasformazione di un layer vettoriale di linee in punti e unione di due layer.
  • Generazione di un DTM Raster a partire da un layer vettoriale GRASS di punti.
  • Modifica della mappa di colori del raster GRASS che rappresenta il DTM.
  • Terrain analysis.
• Digitalizzazione dell'elettrodotto
  • Generazione di un nuovo shapefile puntuale dei sostegni
  • Importazione e visualizzazione sostegni in GRASS (etichette).
  • Determinazione quota dei sostegni in base a DTM e aggiunta coordinate X,Y,Z sostegni nella tabella degli attributi.
  • Generazione tracciato (vettore lineare) da sostegni (vettore di punti): approccio manuale
  • Generazione tracciato (vettore lineare) da sostegni (vettore di punti): approccio automatizzato (con script v.to.lines)
• Generazione di nuovi location / mapset e riproiezione di layer in diversi sistemi di coordinate
• Generazione nuovi Location / Mapset nello stesso Database
• Riproiezione di un vettore
• Importazione e primo filtraggio edificato
• Determinazione buffer e grigliato
  • Buffer intorno al tracciato della linea
  • Grigliato regolare di punti
  • Estrazione dei punti del grigliato all'interno del buffer a 300 m
  • Quotatura grigliato in base a DTM
• Modellazione elettrodotto e calcolo del campo magnetico
  • Modellazione elettrodotto
  • Calcolo del campo magnetico
• Modellazione elettrodotto e calcolo del campo magnetico
  • Visualizzazione shp puntuale con "simbolo graduato", "colore continuo", "valore univoco"
  • Individuazione edifici con altezza superiore a 7 m con livelli di esposizione localmente maggiori di 3 microtesla (approccio 1)
  • Suddivisione dell'edificato con altezza superiore a 7 m in porzioni di edificio differenziate in base ai livelli di esposizione, approccio 1
  • Suddivisione dell'edificato con altezza superiore a 7 m in porzioni di edificio differenziate in base ai livelli di esposizione, approccio 2
  • Individuazione edifici con altezza superiore a 7 m con livelli di esposizione localmente maggiori di 3 microtesla (approccio 2)
  • Suddivisione dell'edificato con altezza superiore a 7 m in porzioni di edificio differenziate in base ai livelli di esposizione, approccio 3
• Visualizzazione 3D
• Elaborazione 3D

Questo tutorial è stato creato come supporto alle 10 ore di esercitazione tenute da Nicola Zoppetti nell'ambito del Master di secondo livello "Sistemi Informativi Geografici per il monitoraggio e la gestione del territorio" 2006-2007. Il tutorial è stato ulteriormente sviluppato nell'ambito delle succssive edizioni del master fino a giungere alla forma attuale. Le esercitazioni seguono 2 ore di seminario sul tema dell'inquinamento elettromagnetico a bassa frequenza ed in particolare sulle linee elettriche ad alta tensione.
Nel seminario ci si è concentrati sugli aspetti fisici e tecnici riguardanti gli elettrodotti ed i campi da essi generati.
Nelle esercitazioni invece si affronta un caso studio utilizzando strumenti GIS. L'aspetto del calcolo dei campi viene solo accennato, trattandosi di un tema che richiederebbe più ore per essere affrontato. Ci si concentra invece sulla preparazione dei dati necessari al calcolo e sull'elaborazione dei risultati. Per vari motivi si è scelto di sviluppare le esercitazioni vincolandosi al sistema operativo Microsoft Windows.
In ogni caso QGIS è disponibile anche per unix e per MAC-OS. Il carattere della esercitazione vuole essere estremamente applicativo: non ci si soffermerà su aspetti teorici o tecnici, pur fondamentali, riguardanti i GIS in generale o QGIS-GRASS in partcolare. Piuttosto si illustrerà come possono essere raggiunti determinato obiettivi rimandando per gli approfondimenti e per gli aspetti generali o di principio, ad opportuni riferimenti.

Per le esercitazioni si è scelto di utilizzare il software QGIS. Si tratta di un software Open Source, disponibile gratuitamente per vari sistemi operativi. Il software può essere scaricato qui. Il programma è stato originariamente sviluppato per linux e da meno tempo è possibile utilizzarlo anche in windows. Le versioni per windows, essendo meno mature sono anche anche meno stabili di quella per linux.
Si tratta in ogni caso di un programma in continua evoluzione e dalle grandi potenzialità e quindi vale la pena provare ad utilizzarlo. Attualmente QGIS, che è nato come GIS leggero e con funzionalità di base, può essere utilizzato insieme a GRASS, un GIS tra i più complessi ed evoluti, e ciò ha ampliato di molto le possibilità di analisi. Quanto detto per QGIS, vale a maggior ragione anche per GRASS, visto che il suo "porting" per windows è relativamente recente.
Più recentemente QGIS sta evolvendo come GIS a sé, grazie anche alla possibilità di sviluppare dei plugin in modo relativamente semplice; in questa esercitazione tuttavia ci concentreremo sull'utilizzo di QGIS come interfaccia per GRASS. Questo tutorial verrà auspicabilmente aggiornato ed esteso sia per illustrare eventuali nuove funzionalità, sia per affrontare i molti aspetti e le possibili elaborazioni legate al problema della valutazione dell'impatto ambientale del campo elettrico e magnetico generati da elettrodotti che non si è avuto tempo di trattare fino ad ora.

• Sito ufficiale QGIS
• Blog QGIS
• Download ultima versione QGIS
• Download precedenti versioni QGIS
• Wiki QGIS: indice dei titoli

• GFOSS ITALIA: Geospatial Free Open Source Software
• GFOSS: HOWTOs e procedure tipo

• Sito ufficiale GRASS
• Sito utenti GRASS italia
• Lista utenti GRASS italia
• Ottimo tutorial GRASS
• Manuali varie versioni di GRASS

Se non è già presente, installare la versione di QGIS usata nelle esercitazioni del 2015 sul proprio computer semplicemente eseguendo il programma di installazione NB: alla fine dell'installazione eliminare la spunta dalla check-box che limita l'uso del programma per motivi di sicurezza. Se non lo si fa GRASS potrebbe non funzionare.

Creare una nuova cartella di lavoro per le esercitazioni. Ad esempio C:/EsercitazioniQGIS2015/

Se si vuole salvare il tutorial , scaricare il file Tutorial HTML (zip 9 MB) ed estrarre tutti i files nella cartella appena creata (utilizzare l'opzione "Use folder names" di WinZip) Nota: i link ai file .zip non funzioneranno nella versione locale.

Scaricare il file Dati.zip (zip 24 MB) ed estrarre tutti i files nella cartella esercitazioni (utilizzare l'opzione "Use folder names" di WinZip)

Cliccando con il tasto destro del mouse su una zona vuota della barra del menù compare la finestra che permettere di scegliere le barre degli strumenti desiderate.

Per poter seguire meglio il tutorial in aula riportiamoci tutti ad una disposizione delle barre degli strumenti simile a quella in figura.

GRASS, come molti altri sistemi GFOSS, usa la libreria PROJ.4 per gestire i sistemi di riferimento, le proiezioni e le loro trasformazioni. Le PROJ.4 definiscono i sistemi di riferimento in base alle loro differenze rispetto al WGS84 ("il datum del GPS"): impostare un sistema di riferimento significa quindi scegliere i parametri della trasformazione da un sistema al WGS84. In teoria si usa una trasformazione a 7 parametri (tre traslazioni + tre rotazioni + 1 fattore di scala), per come sono realizzati i sistemi di riferimento spesso si hanno errori inaccettabili e si usano griglie di valori, i parametri cambiano cioè da punto a punto. link www alla presentazione pdf "Sistemi di riferimento e proiezioni in GRASS" di Paolo Zatelli

Dopo aver lanciato per la prima volta QGIS, impostare un nuovo sistema di proiezione: menù Impostazioni / SR personalizzato / scheda Definizione

Nome: Gauss-Boaga corretto Parametri: +proj=tmerc +lat_0=0 +lon_0=9 +k=0.999600 +x_0=1500000 +y_0=0 +ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

menù Progetto / Proprietà del progetto /

scheda Generale unità mappa: metri precisione: manuale, 3 posizioni decimali

menù Progetto / Proprietà del progetto / scheda SR abilita la proiezione al volo (check box in alto) user defined ... / Gauss-Boaga corretto

Salvare il progetto menù file / salva creare una nuova cartella C:/EsercitazioniQGIS2015/progettiqgis dove salvare tutti i progetti di QGIS

se non è gi à abilitato (ma dovrebbe esserlo) abilitare il plugin GRASS questo passaggio è utile per rendersi conto dei plugin per qgis disponibili menù Plugin / Gestione Plugins

selezionare GRASS

menù Plugin / Grass / New Mapset

Creare una nuova cartella "DatasetGrass" in C:/EsercitazioniQGIS2015/)

Selezionare la cartella appena creata e confermare

Next

Creare una nuova location "ToscanaGaussBoaga" Next

Selezionare la proiezione definita precedentemente Next

Impostare la regione GRASS sull'Italia Next

Nuovo mapset : Suggerimenti per l'assegnazione del nome: gli si potrebbe assegnare il proprio nome nel caso si utilizzi un Dataset/location accessibile da più utenti ("Nicola" ad esempio) nel caso di un pc a singolo utente potrebbe essere nominato in base alla particolare attività svolta ("MasterGIS" ad esempio). Next

Finish

Grass lavora su dati organizzati in modo rigido (in genere l'utente non può scegliere dove posizionare i layer grass su cui sta lavorando) Il sistema Dataset/Location/Mapset corrisponde ad un sistema di organizzazione dei file su tre livelli.

menù Layer / aggiungi vettore / aggiungi un raster

oppure usare la barra degli strumenti

Cercare in C:\EsercitazioniQGIS2015\Dati\Raster e selezionare 275020.tif 275030.tif

Il sistema chiede di selezionare un sistema di riferimento per i raster: selezionare quello creato in precedenza.

menù Visualizza / zoom completo

Uno dei due layer probabilmente non comparirà nella mappa principale. Risolveremo in seguito questo problema (nelle prove effettuate il layer diventa visibile solo aumentando molto lo zoom pur non essendo impostata la visualizzazione dipendente dalla scala). I layer compaiono in legenda Con il tasto destro sul nome del layer in legenda è possibile accedere a molte funzionalità, tra le quali l'accesso alle proprietà del layer.

menù Layer / aggiungi tutto alla vista panoramica

individuare la zona di interesse (vedere riquadro nella vista panoramica)

menù Visualizza / nuovo segnalibro : "Ponte del ristoro"

D'ora in poi, ogni volta che si vuole visualizzare la zona di interesse sarà sufficiente utilizzare il segnalibro appena creato: menù Visualizza / Mostra segnalibri...

Zoom a

La regione di interesse è coperta dal solo layer 275030 per cui si può eliminare il layer 275020. Eliminare non vuol dire cancellare il file dal disco ma semplicemente non prenderlo in considerazione nel progetto corrente. selezionare il layer nella legenda: menù Layer / elimina layer oppure ctrl+D oppure legenda -> tasto dx sul layer -> rimuovi

menù Plugins / Grass / apri strumenti Grass

Toolbox GRASS -> Moduli GRASS / Gestione File / Importa in Grass / Import raster in GRASS / ... dalla vista QGIS/r.in.gdal.qgis

Selezionare il layer 275030. Nelle opzioni avanzate assicurarsi che sia spuntata l'opzione 'sovrascrive la proiezione' Esegui

i raster grass hanno una risoluzione molto spinta (pixel di 0.5mx0.5m) questo e' da imputare al mancato riconoscimento della risoluzione originale del TIFF (vedi metadati TIFF, risoluzione != TIFFTAG_XRESOLUTION) provare a ricampionare il raster una volta importato in grass, anche senza caricare nella vista il raster a piena risoluzione impostare la risoluzione della regione di lavoro usare il comando grass r.resample

La Toolbox GRASS è composta da diverse schede: Albero moduli:elenco dei moduli raggruppati per funzionalità Lista moduli:elenco dei moduli (non è l'ordine), in questa scheda sono possibili ricerche. Browser:permette di navigare nel dataset GRASS. Schede dei comandi: a loro volta suddivise in 3 sottoschede: opzioni: si specificano i parametri del comando output: si visualizza il log di esecuzione del comando manuale: si visualizza il manuale del comando grass

Toolbox GRASS / Browser / Aggiungere il raster grass alla vista Toolbox GRASS / Browser / Imposta la regione corrente con la mappa selezionata

Rimuovere il layer omonimo (tif) dalla legenda Zoom al layer

Selezionare 275030 nella legenda e Nella vista panoramica

Per visualizzare la regione di lavoro di Grass (Rettangolo rosso) utilizzare l'icona evidenziata nella seguente figura assicurandosi che lo strumento sia selezionato

se non si visualizza un rettangolo rosso intorno alla mappa e' necessario aggiornare le coordinate della regione di lavoro in base a quelle del raster: Toolbox GRASS -> Browser -> Imposta regione attuale con mappa selezionata

Salvare il progetto menù progetto / salva

caricare i layer shp vettoriali delle curve di livello e dei punti quotati menù Layer / Aggiungi vettore / aggiungi vettore

CARTOGRAFIA_BASE/275030_PuntiQuotati.shp CARTOGRAFIA_BASE/275030_CurveLivello.shp

Come nel caso dei raster il sistema chiede di selezionare un sistema di riferimento ...

Attivare lo strumento ordine di visualizzazione Cambiare ordine di visualizzazione nella legenda: trascinare con il puntatore del mouse il raster 275030 in cima alla lista (scompaiono i punti e le curve di livello) trascinare con il puntatore del mouse il raster 275030 in fondo alla lista (ricompaiono i punti e le curve di livello)

Doppio clic nella legenda su 275030 (oppure selezionare 275030 e poi Tasto destro / Proprietà) Trasparenza 50% OK

selezionare nella legenda 275030_CurveLivello Toolbox GRASS /Moduli GRASS /Gestione File / Importa in GRASS / Import vettore in GRASS /v.in.ogr.qgis

OGR vector layer: 275030_CurveLivello.shp name for output...:c275030

Esegui

OGR vector layer: 275030_PuntiQuotati.shp name for output...:p275030 Esegui

Selezionare gli shp nella legenda (anche insieme) layer /elimina layer oppure CTRL+d oppure Tasto destro sulle voci selezionate in legenda -> Rimuovi

aggiungere i layer vettoriali GRASS al prgetto ci sono 3 modi: Plugins / GRASS / aggiungi vettoriale GRASS icona aggiungi vettoriale GRASS Toolbox GRASS / Browser / selezionare icona Aggiungi la mappa selezionata all'area di lavoro

Meglio aggiungere solo i sub-layer 1-line e 1-point onde evitare di dover eliminare successivamente gli altri layer, al momento non significativi.

Toolbox GRASS /Moduli Grass / Vettore / Elabora mappa /Gestisci elementi / v.to.points

input: c275030 output:pc275030 distance 10 m Esegui

Toolbox GRASS / Browser / icona Aggiungi la mappa selezionata all'area di lavoro pc275030 1 point 2 point

Approfondimento: al layer 1 point è associata la stessa tabella di database associata alle curve di partenza al layer 2 point è associata una nuova tabella di database (creata dal comando v.to.points) che descrive la posizione del punto nella curva di origine.

Per visualizzare la tabella degli attributi:

eliminare layer c275030 e pc275030 2 point

unione di due layer puntuali: Toolbox GRASS / Vettore / analisi spaziali / sovrapposizione /v.patch

Fare due prove, una selezionando ed una non selezionando la check box "Copy also attribute table" input=p275030,pc275030 output1=ppc275030noattributes output2=ppc275030

Nel secondo caso si ottiene una mappa puntuale a cui è associata una tabella degli attributi completa.

visualizzare la regione corrente (riquadro rosso)

icona modifica regione grass attuale o del menù plugins/grass/modifica regione grass attuale) impostare la regione di lavoro adattata al layer dei punti considerato impostare una risoluzione di 10m sia per le righe sia per le colonne. La risoluzione ed il corrispondente numero di righe e colonne nella regione di lavoro, determinano la dimensione dei raster generati con gli strumenti di grass. Se si scelgono risoluzioni che comportano immagini con molti pixel, i tempi di calcolo si allungano di molto ed aumentano le risorse di memoria necessarie). Con normali pc si hanno tempi di calcolo ragionevoli con raster di 1000x1000 pixel.

Meglio se impostare le coordinate che definiscono i limiti della regione con dei multipli della risoluzione desiderata

Toolbox GRASS / Moduli / raster / surface management / interpolate surface v.surf.idw

input: ppc275030 1point campo attributi: VALORE nome output: dtm_10m_12pt numero punti: 12

Aggiungere il raster grass creato (icona o menù Plugin / GRASS /aggiungi raster)

Se il raster appare colorato uniformemente come nella precedente immagine fare doppio clic sul nome del raster nella legenda (oppure tasto destro proprietà) ed impostarle come nella seguente immagine (Stira a MinMax e poi Applica).

L'immagine sarà colorata in base alla quota attribuita a ciascun pixel secondo una mappa di colori in scala di grigi.

Ogni comando di GRASS, oltre che dalla toolbox, può essere lanciato utilizzando la shell di grass cioé una finestra di comandi. Si crea un nuovo dtm utilizzando un numero diverso di punti (4 in questo caso) nella GRASS shell v.surf.idw input=ppc275030 column="VALORE" output=dtm_10m_4pt npoints=4

Si crea un nuovo raster in cui sono riportate l differenze tra i due dtm: nella GRASS shell r.mapcalc delta_dtm_idw12pt_idw4pt=dtm_10m_12pt-dtm_10m_4pt

La stessa operazione (r.mapcalc) può essere impostata per via grafica: Toolbox GRASS / Albero moduli / Raster / spatial analysis / map algebra dove nel caso illustrato calcolo il valore assoluto alle differenze (si possono apprezzare i differenti massimi e minimi nel riquadro destro nella scheda browser dei layer nella toolbox grass). Il calcolatore grafico permette anche di salvare l'operazione e di richiamarla in seguito. In caso contrario il set di operazioni va perso alla chiusura della toolbox.

Nel browser si possono notare le differenze tra i due raster ottenuti rispettivamente da riga di comando e con il tool grafico (differenza e valore assoluto della differenza)

Nella finestra proprietà è possibile impostare l'estensione del gradiente di colore in scala di grigi, regolando il contrasto.

NB: TUTTE LE MODIFICHE ALLA MAPPA DI COLORE EFFETTUATE IN GRASS (r.colors sia da shell sia da grass-toolbox) HANNO EFFETTO IN QGIS 2.6.1.SOLO SE EFFETTUATE SU UN RASTER NON CARICATO NELLA VISTA DI QGIS!!!

Aprire la scheda scheda metadati del dtm_10m_12pt Tasto destro / proprietà / Metadati Si nota che le quote sono comprese tra 28 e 153 m s.l.m.

Utilizzare il comando r.colors (lanciarlo da shell) utilizzando "valori percentuali" r.colors map=dtm_10m_12pt color=rules 0% 0 0 0 20% 0 100 0 40% 0 255 0 60% 0 100 100 80% 0 255 255 100% 255 0 0 end

Invece di scrivere sulla shell le regole, è possibile memorizzarle una volta per tutte su un file di testo In particolare creare un file di nome color1.file con il seguente contenuto 0% 0 0 0 20% 0 100 0 40% 0 255 0 60% 0 100 100 80% 0 255 255 100% 255 0 0 end

Da shell, lanciare il comando: cat C:/EsercitazioniQGIS2013/work/color1.file | r.colors map=dtm_10m_4pt color=rules

E' possibile lanciare il comando dall'interfaccia grafica del gestore strumenti di GRASS: Toolbox Grass -> raster -> r.colors.rules in tal caso eliminare end finale dal file delle regole

E' necessario eliminare end finale dal file delle regole altrimenti si ottiene il seguente errore

Dopo aver eliminato end finale dal file delle regole altrimenti il comando va a buon fine

E' inoltre possibile utilizzare il comando r.colors utilizzando "valori assoluti" r.colors map=dtm_10m_12pt color=rules 0 0 0 255 50 0 255 0 70 0 100 0 100 50 100 0 130 100 100 0 200 255 100 0 end

Durante le prove con QGIS 2.6.1 le modifiche effettuate in grass hanno avuto effetto solo dopo aver scaricato e ricaricato il layer dalla vista, portando al seguente risultato.

In QGIS1.7 è stata aggiunta la possibilità di modificare in modo semplice la mappa di colore utilizzando le schede di proprietà del raster tasto dx -> proprietà -> stile tasto dx -> proprietà -> mappa colore

Ci si può poi divertire utilizzando i vari strumenti GRASS dedicati alla "Terrain Analysis".

Raster -> analisi spaziali -> modello di irradiazione solare -> .... esatta posizione del sole

con QGIS2.6.1 + grass si ottiene una mappa di tutti 1 (ERRORE)

Raster -> analisi spaziali -> analisi geomorfologica -> r.slope (PENDENZE)

Raster -> analisi spaziali -> analisi geomorfologica -> r.slope.aspect (ESPOSIZIONE)

Cambio la mappa di colore degli aspect scaricando il raster dalla vista nella shell r.colors map=aspects color=rules 0 0 0 255 45 0 0 255 46 0 255 0 135 0 255 0 136 255 0 0 225 255 0 0 226 255 255 0 315 255 255 0 316 0 0 255 360 0 0 255 end

Impostando opportunamente la trasparenza dei raster (doppio clic sul nome nella legenda o tasto destro + proprietà) è possibile sovrapporre più immagini e ottenere mappe composite (nell'immagine la sovrapposizione tra il dem, la mappa delle pendenze ).

menù Layer / crea vettore /nuovo layer shapefile

tipo: punto attributi: X real Y real Z real h real tipo string codice string

OK salvare lo shapefile nella cartella work con nome sostegni_407

Visualizzare il raster 275030 Aggiungerlo alla panoramica Visualizzare la zona illustrata in figura

attivo la barra di digitalizzazione attivare strumento modifica (icona con pennarello)

icona inserisci punti con il mouse cliccare nel punto in cui sul raster è rappresentato il primo sostegno della linea (all'interno della sottostazione di S.Lorenzo a Greve) inserire gli attributi nella tabella (lasciare vuoti i campi x,y,z).

modificare la simbologia dell' shp in modo che sia più evidente.

utilizzando una scala di visualizzazione che permetta l'inserimento dei sostegni sulla mappa con una precisione sufficiente, alternare l'inserimento dei nuovi sostegni e il pan lungo la linea.

una volta arrivati al fiume arno terminare l'inserimento (17nsostegni) icona toggle editing salvare le modifiche

evidenziare sostegni_407 nella legenda menù Visualizza / zoom sul layer

evidenziare sostegni_407 nella legenda tasto destro / visualizza tabella attributi oppure icona Apri tabella

correggere gli eventuali errori di inserimento Abilita modifica doppio clic sul campo da modificare alla fine di tutte le modifiche Interrompi modifica Salva

importare l'shp sostegni_407 come layer GRASS sost407 apri strumenti GRASS / v.in.ogr eliminare il layer sostegni_407.shp caricare il layer grass sost407 (icona o menù Plugin / GRASS / Aggiungi vettoriale GRASS)

Aggiunta etichette ai sostegni evidenziare sost407 nella legenda tasto destro / proprietà oppure doppio clic sul nome del layer Scheda etichette: campo con etichetta: codice spostamento unità: punti

caricare il raster GRASS dtm_10m_12pt creato in precedenza

Toolbox GRASS moduli / vettore / elabora mappa / v.drape

input: sost407 output: sost407_3d raster: dtm_10m_12pt sampling: a piacere, se c'è tempo fare un po' di prove

Caricare il layer vettoriale GRASS sost407_3d Ancora la tabella degli attibuti ha i campi coordinata vuoti.

SCARICARE IL VETTORE DALLA VISTA (solo qgis2.6.1) Toolbox GRASS Modules / Vectors / Report and statistics / v.to.db

Con questo comando è possibile agiungere alla tabella degli attributi le coordinate dei punti. l'interfaccia del comando v.to.db permette di specificare solo un campo di db in cui inserire le coordinate

si usa la shell v.to.db map=sost407_3d option=coor units=me column=X,Y,Z Ricaricare il layer vettoriale sost407_3d nella vista Si controlla che nella tabella degli attributi siano effettivamente stati aggiornati i campi X,Y,Z

NB: il baco descritto di seguito è stato risolto nelle versioni successive a quella considerata durante la prima stesura del tutorial. Tuttavia si mantiene la spiegazione perchè si ritiene possa essere di qualche utilità Il programma non aggiorna la tabella dbf associata al layer vettoriale ma crea solo una tabella temporanea. si deve sostituire manualmente la tabella associata con quella temporanea. Per farlo: si naviga con esplora risorse fino alla cartella .tmp nel mapset in cui stiamo lavorando. Si ordinano i file secondo la data di ultima modifica Si copia in clipboard l'ultimo dbf creato Si incolla il dbf nella cartella dbf, allo stesso livello di .tmp si rinomina: sostegni_407_3d ---> sostegni_407_3d_old si rinomina 2456 (file temporaneo) ---> sostegni_407_3d Si visualizza la tabella degli attributi in cui i campi X, Y e Z sono correttamente riempiti (Se ciò non accade provare a scaricare e ricaricare il layer sostegni_407_3d affinché l'applicazione si "accorga" delle modifiche alla tabella dbf.

tasto destro nella legenda su sostegni_407_3d esporta come shapefile salvare sostegni_407_3d nella cartella work un'ulteriore opportunità di esportazione è data dallo strumento grass toolbox GRASS / File / Export / Export vectors / v.out.ascii

la personalizzazione del separatore (QGIS1.7) non funziona se si usa un carattere 'spazio' lanciando il comando da shell è possibile personalizzare il carattere separatore (funziona anche con lo spazio) ed altri particolari. salvare il layer grass sostegni_407_3d in file di testo sia in formato point sia in formato standard salvare in due file diversi nella cartella work con i seguenti nomi: s407point.txt s407std.txt

aprire con wordpad il file s407_point.txt sostitire il carattere | con uno spazio per evitare questo passaggio si poteva lanciare da shell il comando di esportazione in formato testo specificando il separatore spazio o tab: v.out.ascii input=sost407_3d format=point output=C:/EsercitazioniQGIS2015/WORK/s407point.txt fs="\t"

aprire con wordpad il file di testo s407std.txt cancellare tutto il contenuto del file dopo la 10 riga VERTI(esclusa) aggiungere la riga L 16 1 copiare dopo l'ultima riga tutto il contenuto di s407_point.txt inserire come ultima riga: 1 1

dalla shell di GRASS lanciare il comando v.in.ascii input=C:/EsercitazioniQGIS2015/WORK/s407std.txt output=tracciato407 --overwrite format=standard

caricare il layer vettoriale tracciato407

In QGIS1.7 l'interfaccia grafica del comando v.in.ascii è stata molto potenziata. Adesso non è necessario usare la shell.

Il precedente approccio può essere automatizzato mediante la generazione di uno script. Non si sta qui ad entrare nel merito di come è possibile generare uno script da utilizzare in grass. Lo script già fatto si trova nella cartella Dati/varie e si chiama v.to.lines Copiare il file v.to.lines ed incollarlo nella cartella C:\Programmi\QGIS Brighton\apps\grass\grass-6.4.2\scripts Copiare il file v.to.lines.bat ed incollarlo nella cartella C:\Programmi\QGIS Brighton\apps\grass\grass-6.4.2\bin Lo script richiede che si specifichi: il nome del layer puntuale di ingresso il nome del layer delle linee in uscita Una volta copiato lo script provare a digitare nella shell: v.to.lines input=sost407_3d output=tracciato407bis

Invece di usare la schell potevo usare l'interfaccia grafica standard che grass crea per il nuovo comando alla quale si accede digitando nella shell il nome del nuovo comando v.to.lines

Per generare lo script v.to.lines ho utilizzato i seguenti link: Messagio 1 nella lista utenti grass Messagio 2 nella lista utenti grass Tutorial awk

Se poi si vuole aggiungere alla toolbox grass l'interfaccia grafica riferirsi al seguente link: QGIS: Grass Tools aggiunta moduli

In particolare: copiare i file v.to.lines.qgs v.to.lines.1.png e v.to.lines.2.png nella cartella C:\Programmi\QGIS Brighton\apps\qgis\grass\modules modificare il file C:\Programmi\QGIS Brighton\apps\qgis\grass\config\default.qgc aggiungendo sotto grass name="v.to.points"/ la riga grass name="v.to.lines"/

Affinché le modifiche abbiano effetto, uscire da QGIS e rilanciarlo

menù Plugin / Grass / New Mapset

Riferirsi allo stesso Database "Datigrass" (e quindi cartella) creato in precedenza Next

Creare una nuova location "ToscanaUTMED50" Next

Selezionare la proiezione ED50/utm ZONE 32N Next

Eventualmente cercare per id della proiezione

Nuovo mapset :"nicola" Next

menù Plugin / Grass / New Mapset

Riferirsi allo stesso Database "Datigrass" (e quindi cartella) creato in precedenza Next

Creare una nuova location "ToscanaUTMWGS84" Next

Selezionare la proiezione WGS84/UTM zone 32N Next

Nuovo mapset :"nicola" Next

Riferirsi allo stesso Database "Datigrass" (e quindi cartella) creato in precedenza Next

Creare una nuova location "WGS84" Next

Selezionare la proiezione WGS84 Next

Nuovo mapset :"nicola" Next

Aprire il primo mapset creato ToscanaGaussBoaga (coordinate Gauss-Boaga) Menù plugins -> Grass -> Apri mapset

importare con v.in.ogr lo shapefile C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\sostegni_40700.shp

il file importato ha la tabella completa e rappresenta tutti i sostegni della linea

menù Plugin / Grass / apri mapset Location: WGS84 Mapset: MasterGIS

menù Plugin / Grass / Grass tools Vettore /elabora mappa / riproietta ... / v.proj

In QGIS1.7 è necessario usare la shell v.proj mapset=nicola location=ToscanaGB input=s407 output=s407wgs84

Alternativamente digitare nella shell solo v.proj e compare l'interfaccia grafica del comando di grass.

I vettori sono riproiettati al volo nel sistema di coordinate del progetto (sempre che sia attivata l'opzione della riproiezione al volo nelle proprietà del progetto).

Cambiare il sistema di coordinate del progetto da Gauss Boaga a WGS84

Esportare con il comando v.out.ascii entrambi i file s407 e verificare che le coordinate siano diverse.

Se nel sistema WGS84 i due vettori non fossero allineati, controllare quale è il sistema di proiezione di base del layer in coordinate Gauss-Boaga. Legenda / tasto destro / proprietà / scheda generale / Cambia

Impostare il sistema di proiezione Gauss-Boaga

I vettori adesso dovrebbero risultare allineati

Ad ulteriore verifica che il file di proiezione prj non permette di definire correttamente gli scostamenti rispetto al sistema geografico WGS84, è possibile esportare come shp il vettore in coordinate Gauss-Boaga (tasto destro sulla legenda e "salva come shapefile ..."). Se si importa lo shape nella vista i punti risultano nuovamente sfalsati rispetto a quelli del vettore di origine.

A differenza del formato prj (ESRI), il più generco formato wkt, permette di specificare correttamente il rapporto tra un generico sistema di proiezione ed il sistema WGS84. Per creare un file wkt relativo ad una certa proiezione: Aprire un mapset di una location relativa alla proiezione di interesse (ad esempio ToscanaGaussBoaga / nicola) nella shell di grass digitare g.proj -p sono stampate sulla shell le caratteristiche della proiezione per stamparle in formato wkt digitare g.proj -w copiare la risposta in un file di testo vuoto Dare al file estensione wkt e lo stesso nome dello shape a cui si vuole associare

si imposta la regione di lavoro intorno all'alettrodotto appena digitalizzato con ampio margine W: 1674600 E: 1679700 S: 4847300 N: 4851900

si carica nella vista lo shapefile dell'edificato che sta nella cartella CARTOGRAFIA_BASE se si pianta tutto passare al punto successivo

Si tratta di uno shapefile molto pesante di cui a noi interessa solo una porzione.

Toolbox GRASS -> v.in.ogr

selezionare l'opzione "Limit import to current region" opzione eliminata in QGIS1.7 quindi usare la shell (vedi sotto) o digitare v.in.ogr nella shell e utilizzare l'interfaccia grafica di grass.

oppure dalla shell di GRASS lanciare il comando v.in.ogr dsn=C:\EsercitazioniQGIS2013\Dati\CARTOGRAFIA_BASE layer=edificato_ab_2001_fi output=edificato spatial=1674600,4847300,1679700,4851900 l'opzione spatial=xmin,ymin,xmax,ymax permette di importare nel layer GRASS edificato soli gli elementi all'interno del rettangolo definito dalle coordinate di due suoi angoli Noi utilizziamo le stesse coordinate con cui abbiamo definito la regione di lavoro. E' necessario utilizzare l'opzione -o per riproiettare correttamente l'shp.

digitando v.in.ogr nella shell e utilizzando l'interfaccia grafica di grass ...

dal browser della toolbox GRASS si carica il layer 1 polygon dalla mappa GRASS edificato

la mappa dell'edificato è adesso molto più leggera e gestibile. Verificare che la tabella degli attributi abbia conservato i campi originari.

Il layer 0 polygon è principalmente collegato a edifici "con buco".

L'obiettivo di questa sezione è quello di ricavare una griglia regolare di punti (quotati in base al dtm) limitandosi ad un'area a distanza limitata dal tracciato della linea elettrica.
Questi punti sono quelli in cui le applicazioni per il calcolo del campo magnetico calcoleranno il campo stesso.

Caricare nella vista il tracciato della linea creato in precedenza nome layer: tracciato407

Creare un buffer di 400m intorno al tracciato della linea Toolbox GRASS / Modules / Vettore / Analisi spaziali /Buffer / v.buffer

con la setessa procedura: Creare un buffer di 300m intorno al tracciato della linea

si può usare la shell: v.buffer input=tracciato407 output=buffer400 distance=400 v.buffer input=tracciato407 output=buffer300 distance=300

Impostare la regione di lavoro con l'estensione e la risoluzione desiderata per la griglia di punti. E' molto importante stringere ancora i limiti della regione intorno al buffer per evitare di creare grigliati con troppi punti e quindi operazioni eccessivamente onerose. Provare i prossimi 2 passaggi impostando una risoluzione a 5m. Se il comando r.tovect fosse troppo lento impostare una risoluzione a 10m.

oppure

dalla shell di GRASS r.mapcalc surface=1

dalla shell di GRASS r.to.vect feature=point in=surface out=ptgriglia5m oppure r.to.vect feature=point in=surface out=ptgriglia10m

Se la visualizzazione del grigliato appena generato fosse troppo lenta e più in generale per gestire la visualizzazione di mappe molto pesanti eliminare la spunta dalla check box "visualizza". Ciò permette di effettuare varie operazioni (ad esempio la modifica della simbologia di un layer o dell'ordine di visualizzazione) senza che la vista sia ridisegnata. cambiare l'ordine di visualizzazione delle mappe aumentare lo zoom (la vista panoramica continua ad essere aggiornata) rimettere la croce nella check box "disegna" e visualizzare i punti della griglia

Come esercizio (x casa): provare a eseguire le ultime due operazioni con interfaccia grafica (determinazione surface e trasformazione in vettoriale)

Toolbox GRASS / Modules / Vettore / Analisi spaziali / estrai elementi da un vettore / v.select.overlap

input: ptgriglia5m input: buff300m output: ptgriglia5mb300m ATTENZIONE, l'operazione è abbastanza lenta.

eliminare la spunta dalla check box "disegna" cambiare l'ordine di visualizzazione delle mappe aumentare lo zoom (la vista panoramica continua ad essere aggiornata) impostare la trasparenza (50-60%) dei buffer rimettere la croce nella check box "disegna" e visualizzare i punti della griglia filtrati all'interno del buffer a 300m

rendere un layer grass puntuale a 3D in base al dtm da shell v.drape input=ptgriglia5mb300m type=point rast=dtm_10m_12pt method=bilinear output=ptgriglia5mb300m3D

oppure con la Toolbox GRASS Toolbox GRASS / Modules / Vettore / Elabora mappa / Converti un vettore 2D ...

oppure con l'interfaccia di Grass delcomando v.drape (digitare v.drape nella shell)

esportare il layer grass come testo da shell v.out.ascii input=ptgriglia5m3D_buff300m format=point output=C:\EsercitazioniQGIS2013\WORK\grigliato_407_300m.txt fs="\t"

oppure usare l'interfaccia Grass al comando v.out.ascii

Questa parte del tutorial è solo esemplificativa-informativa e non si forniranno gli strumenti di modellazione e di calcolo. Vuole essere solo un esempio di come si procede con una delle applicazioni utilizzabili per la modellazione di elettrodotti ed il il calcolo dei campi.

esportare il layer sostegni_407_3d come file di testo formato point nella cartella DatiElettrodotti v.out.ascii input=sostegni_407_3d format=point output=C:/EsercitazioniQGIS2013/Dati/DatiElettrodotti/sostegni_407_3d.txt fs="\t"

In questo file sono riportate solo le coordinate dei sostegni ma non le altre informazioni (tipo e altezza) inserite al momento della digitalizzazione (tabella degli attributi).

Esporto la tabella dbf come file di testo echo "select * from sostegni_407_3d"|db.select fs=" " > C:\EsercitazioniQGIS2013\WORK\sostegni_407_3d.txt

Carico il file appena creato nel programma di modellazione

Assegno il tipo al primo sostegno (il portale)

Assegno il tipo al gruppo di sostegni 2-7. Completo l'assegnazione dei tipi come nei due casi precedenti, con riferimento alla tabella testuale.

Una volta assegnati i tipi, si modella il cablaggio della linea.

Importare il flie di testo appena creato nell'applicazione di calcolo del campo magnetico.

Calcolare il campo magnetico a diverse altezze da terra (quota del terreno definita dalla coordinata Z dei punti del grigliato)

Determinare la fascia di rispetto dei primo livello

Calcolare alcune curve isocampo a 1 m da terra considerando un modello del terreno semplificato

II risultati del calcolo sono salvati sia sotto forma di file di testo sia come shapefile. Caricare ad esempio lo shapefile relativo ad un calcolo su grigliato (ad esempio 1m da terra)

Modificare la simbologia dei punti nel modo indicato in figura ed aumentare lo zoom.

Caricare ad esempio lo shapefile relativo aL un calcolo su grigliato a 4m da terra. Modificare la simbologia dei punti nel modo indicato in figura.

isolamento edificato con altezza superiore a 7 m: da shell v.extract input=edificato output=edificato7m where="heigth>7" --o NB: heigth !!!!

OPPURE Toolbox GRASS / Modules /Vector/Spatial analysis /v.extract.where

importazione grigliato a 7 m, e chiamare il layer Bgriglia_7m

isolamento punti di calcolo caratterizzati da campo maggore di 3 uT da shell v.extract input=Bgriglia_7m output=Bgriglia_7m_3uT where="B>3" --o

OPPURE, come in precedenza, Toolbox GRASS / Modules /Vector/Spatial analysis /v.extract.where

Estrazione dell'edificato Con altezza superiore a 7 m con livelli di esposizione localmente maggiori di 3 microtesla: Toolbox GRASS / Modules /Vector/Extract features from vector / v.select.overlap

estrazione dell'edificato on altezza superiore a 7 m con livelli di esposizione localmente maggiori di 3 microtesla in: edificato7m in: Bgriglia_7m_3uT out:edificato_7m_3uT

il modo più semplice per valutare il numero di abitanti esposti è quello di fare una query e copiare la risposta in un foglio di calcolo, in cui effettuare la somma nella shell: echo select ABITANTI01 from edificato_7m_3uT|db.select

generazione di un raster che rappresenta i livelli di campo magnetico se non lo si è già fatto importare la griglia in GRASS Toolbox GRASS / Modules /import vector/v.in.ogr output: Bgriglia_7m

caricare il layer GRASS Bgriglia_7m utilizzare il comando GRASS Toolbox GRASS / Modules / raster /Interpolate surface / v.surf.idw per creare un raster che rappresenti i livelli di campo magnetico a 7 m dal terreno. Dalla shell di GRASS v.surf.idw input=Bgriglia_7m column="B" output=B_7m npoints=4

Oppure dalla Toolbox di GRASS v.surf.idw input=Bgriglia_7m column="B" output=B_7m npoints=4

Toolbox GRASS / Modules /Import / Map Type conversion /v.to.rast Convert vector to raster using constant trasformare l'edificato a 7 m in raster

oppure da shell v.to.rast input=edificato7m@nicola value=1 output=edificato7m use=val

generare un nuovo raster con r.mapcalc in modo che se ci si trova dentro ad un edificio ed il livello di campo è maggiore di 3 uT allora la cella del raster vale 1. da shell: r.mapcalc "ed_7m_3uT=if(edificato7m==1&&B_7m>3,1,null())"

analogamente per 1 < B < 3 : r.mapcalc "ed_7m_1to3uT=if(edificato7m==1&&B_7m>1&&B_7m<3,1,null())" analogamente per 0.4 < B < 1 : r.mapcalc "ed_7m_04to1uT=if(edificato7m==1&&B_7m>0.4&&B_7m<1,1,null())"

riporto tuutti i comandi per comodità, in modo da fare un unico copia incolla dei 3 comandi nella shell: r.mapcalc "ed_7m_3uT=if(edificato7m==1&&B_7m>3,1,null())" r.mapcalc "ed_7m_1to3uT=if(edificato7m==1&&B_7m>1&&B_7m<3,1,null())" r.mapcalc "ed_7m_04to1uT=if(edificato7m==1&&B_7m>0.4&&B_7m<1,1,null())"

modifico le mappe di colori associate ai 3 raster da shell: r.colors map=ed_7m_3uT color=rules 1 255 0 0 end r.colors map=ed_7m_1to3uT color=rules 1 0 0 255 end r.colors map=ed_7m_04to1uT color=rules 1 0 255 0 end

Carico i tre raster nella vista.

genero i raster in cui differenzio 4 regioni in base ai lielli di esposizione (copiare le seguenti righe e poi Maiusc+ clic tasto sinistro nella shell GRASS) r.mapcalc "B_7m_3uT=if(B_7m>3,1,null())" r.mapcalc "B_7m_1to3uT=if(B_7m>1&&B_7m<3,1,null())" r.mapcalc "B_7m_04to1uT=if(B_7m>0.4&&B_7m<1,1,null())" r.mapcalc "B_7m_0to0.4uT=if(B_7m>0&&B_7m<0.4,1,null())"

carico i 3 raster relativi ad esposizioni maggiori di 0,4 microtesla

trasformo i raster in vettori (poligoni) (copiare le seguenti righe e poi Maiusc+ clic tasto sinistro nella shell GRASS) r.to.vect input=B_7m_3uT output=B_7m_3uT feature=area --o r.to.vect input=B_7m_1to3uT output=B_7m_1to3uT feature=area --o r.to.vect input=B_7m_04to1uT output=B_7m_04to1uT feature=area --o visualizzare i vettori grass appena creati

Utilizzando il comando Overlay creo delle mappe vettoriali che rappresentano le porzioni di edificio esposte ad un certo intervallo di vlori di campo (http://grass.itc.it/grass63/manuals/html63_user/v.overlay.html) (copiare le seguenti righe e poi Maiusc+ clic tasto sinistro nella shell GRASS) v.overlay ainput=B_7m_3uT atype=area alayer=1 binput=edificato7m btype=area blayer=1 output=edificato_7m_3uT operator=and olayer=1,0,0 --o v.overlay ainput=B_7m_1to3uT atype=area alayer=1 binput=edificato7m btype=area blayer=1 output=edificato_7m_1to3uT operator=and olayer=1,0,0 --o v.overlay ainput=B_7m_04to1uT atype=area alayer=1 binput=edificato7m btype=area blayer=1 output=edificato_7m_04to1uT operator=and olayer=1,0,0 --o visualizzare i vettori grass appena creati utilizzando dei colori adatti ad evidenziare le elaborazioni fatte

Utilizzando l'approccio precedente alla determinazione degli edifici più alti di 7 m ed esposti a più di 3 microtesla non è pienamente soddisfacente, perché si possono verificare situazioni come quella rappresentata nella seguente figura.

Pure con una risoluzione spaziale dei punti pari a 5m ci sono edifici che non sono considerati perché non sovrapposti ad alcun punto. Per ovviare a questo inconveniente si usa il seguente approccio.

Albero moduli: Raster/Gestione superfici/Genera curve di livello vettoriali/r.contour2 generare le curve di livello vettoriali relative a 3,1 e0.4 microtesla a partire il raster B_7m e chiamarle rispettivamente iso7m04uT iso7m1uT iso7m3uT

oppure da shell: r.contour input=B_7m levels=0.4 output=iso7m04uT r.contour input=B_7m levels=1 output=iso7m1uT r.contour input=B_7m levels=3 output=iso7m3uT

Albero moduli: Vettore/Elabora mappa/Cambia il tipo di elementi geometrici/v.type.lb passare da linee ad aree (senza centroidi): aree7m04uT aree7m1uT aree7m3uT

oppure da shell: v.type input=iso7m04uT output=aree7m04uT type=line,boundary v.type input=iso7m1uT output=aree7m1uT type=line,boundary v.type input=iso7m3uT output=aree7m3uT type=line,boundary

Albero moduli: Vettore/Elabora mappa/Gestisci elementi/v.centroids input: aree7m3uT output: aree7m3uT_cen (v.centroids input=aree7m3uT@nicola output=aree7m3uT_cen option=add cat=1 step=1)

Albero moduli: Vettore/Analisi spaziali/Estrai elementi da un vettore/v.select.overlap input A: edificato_7m 1_polygon input B: aree7m3uT_cen 0_polygon output:edificato_7m_3uT_bis

Come si può notare il problema della mancata selezione di alcuni edifici è stato superato.

Aggiungere i centroidi anche alle aree relative agli 04 e 1 microtesla generando i layer: aree7m04uT_cen aree7m1uT_cen Albero moduli: Vettore/Elabora mappa/Gestisci elementi/v.centroids

da shell: v.centroids input=aree7m04uT output=aree7m04uT_cen option=add cat=1 step=1 v.centroids input=aree7m1uT output=aree7m1uT_cen option=add cat=1 step=1

Creo le aree caratterizzate da valori di campo compresi tra 0,4 e 1 microtesla e tra 1 e 3 microtesla Albero moduli: Vettore/Analisi spaziale/sovrapposizione/v.overlay.xor input A: aree7m04uT_cen 1_polygon input B: aree7m1uT_cen 1_polygon output: aree7m_04_1_uT 1_polygon

Albero moduli: Vettore/Analisi spaziale/sovrapposizione/v.overlay.xor input A: aree7m1uT_cen 1_polygon input B: aree7m3uT_cen 1_polygon output: aree7m_1_3_uT 1_polygon

per suddividere gli edifici esposti in porzioni in base al livello di campo uso nuovamente il comando v.overlay ma con l'opzione and se uso il comando v.overlay nella shell o in Vettore/Analisi spaziale/sovrapposizione/v.overlay.and specificando inputA: edificato7m inputB: aree7m_1_3_uT output: xxx ottengo un errore dovuto al fatto che v.overlay tronca i nomi delle colonne a 8 caratteri e le due colonne ABITANTI01 ed ABITANTI91 vengono a coincidere e conseguentemente il comando fallisce.

Se vogliamo che le porzioni di edificio conservino gli attributi contenuti nella tabella di database, prima è necessario rinominare una delle colonne Albero moduli / Database /Gestione colonne /v.db.renamecol Rinomino la colonna ABITANTI91 con AB91

Una volta modificato il nome della colonna è possibile estrarre delle porzioni di edificio con Vettore/Analisi spaziale/sovrapposizione/v.overlay.and in modo che mantengano le informazioni della tabella associata

Prendo come riferimento le porzioni esposte a livelli compresi tra 1 e 3 microtesla (e quindi il layer edificato7m_1_3_uT) Aggiungo al db la colonna che contiene l'area della porzione di edificio

da shell v.db.addcol map=edificato7m_1_3_uT layer=1 "columns=areaporz double precision" oppure nella toolbox vettore / Database / Gestione colonne / v.db.addcol

aggiorno il contenuto della colonna con l'area delle porzione di edificio

da shell v.to.db map=edificato7m_1_3_uT layer=1 option=area units=meters columns=areaporz oppure nella toolbox vettore / Report e statistiche / v.to.db

stampo sulla shell il contenuto della colonna a_AREA (l'area dell'intero edificio) e areaporz (l'area della porzione appena calcolata)

da shell echo select a_AREA,areaporz from edificato7m_1_3_uT | db.select

aggiungo la colonna abitporz che conterrà gli abitanti della porzione di edificio da shell (vedi sopra l'aggiunta della colonna areaporz) v.db.addcol map=edificato7m_1_3_uT layer=1 "columns=abitporz double precision"

Calcolo il numero degli abitanti della porzione di edificio come il numero degli abitanti contenuti nell'intero edificio moltiplicato per il rapporto tra l'area della porzione con quella dell'intero edificio da shell v.db.update map=edificato7m_1_3_uT column=abitporz value="areaporz/a_AREA*a_ABITANTI" Oppure vettore / Database / Gestione colonne / v.db.update_op

Infine stampo sulla shell il contenuto della colonna a_ABITANTI (gli abitanti dell'intero edificio) e areaporz (l'area della porzione appena calcolata)

da shell echo select a_AREA,areaporz from edificato7m_1_3_uT | db.select

Correzione per tenere conto del numero di piani degli edifici

Importare come layer GRASS, con v.in.ogr i seguenti shp: C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\sostegni_40700_strutturasostegni.shp C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\catenarie_40700_catenarie.shp C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\sostegni_40900_strutturasostegni.shp C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\shp chiamandoli rispettivamente: strutsost407 cat407 strutsost409 cat409 ed avendo cura di importarli come layer 3D (si veda la seguente immagine

Nella ToolboxGrass, aprire il modulo NVIZ

Utilizzare come dtm uno dei raster creati durante le esercitazioni e scegliere uno dei vettori grass 3D appena importati

L'estensione della superficie visualizzata è quella dell'area di lavoro grass attualmente impostata. Con la freccia verde è possibile modificare posizione ed orientamento del punto di vista.

Menù visualize -> Vector lines / 3D polygons New aggiungere i layer delle catenarie e delle strutture dei sostegni di entrambe le linee

Per visualizzare l'edificato è prima necessario "estruderlo" con il comando v.extrude v.extrude input=edificato output=edificato_vol elevation=dtm_10m_12pt hcolumn=heigth

Torno in NVIZ Menù visualize -> Vector lines / 3D polygons New aggiungere il layer edificato_vol

Importare come layer GRASS, con v.in.ogr i seguenti shp: C:\EsercitazioniQGIS2013\Dati\RisultatiCalcolo\02B**,0_punti.shp chiamandoli: B** ed avendo cura di importarlo come layer 3D

Piuttosto che ripetere 15 volte la stessa procedura usare la shell + uno script: cd C:/EsercitazioniQGIS2013/work Creare nella cartella work un file seriecomandi.txt di testo incollandoci il seguente contenuto: v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B24,0_punti_B.shp output=B24 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B22,0_punti_B.shp output=B22 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B20,0_punti_B.shp output=B20 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B18,0_punti_B.shp output=B18 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B16,0_punti_B.shp output=B16 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B14,0_punti_B.shp output=B14 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B12,0_punti_B.shp output=B12 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B10,0_punti_B.shp output=B10 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B8,0_punti_B.shp output=B08 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B6,0_punti_B.shp output=B06 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B4,0_punti_B.shp output=B04 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B2,0_punti_B.shp output=B02 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/02B0,0_punti_B.shp output=B00 -z snap=-1 min_area=0.0001 --overwrite Da shell: bash ./seriecomandi.txt

Unire tutti i layaer appena creati in un unico layer puntuale: C:\>v.patch -e input=B00,B02,B04,B06,B08,B10,B12,B14,B16,B18,B20,B22,B24 output=Bmultiq --overwrite

nella shell: g.region ed impostare i limiti della regione corrente basandosi sul layer appena creato

Spuntare le opzioni che permettono di stampare i limiti della regione modificata

Controllare l'output ed eventualmente aggiustare la risoluzione (evitare valori inferiori a 2 m altrimenti i tempi di calcolo potrebbero allungarsi)

nella shell: v.to.rast3 e quindi nell'interfaccia grafica impostare i campi come in figura

Albero moduli /visualizzazione 3D/ nviz scegliere uno dei raster dtm creati durante le esercitazioni

Visualizza / Volumi / nuovo

Visualizza / Raster / posizione z = -40 da notare che la regione ha 'bottom = 40' Piu' correttamente un effetto analogo si ottiene spostando il volume (posizione Z) a + 40m

Visualizza / Volumi / Aggiungi / Nuova costante / 10

Visualizza / Volumi / Aggiungi / Nuova costante / 5

Visualizza / Volumi / Aggiungi / Nuova costante / 3

Cambiare gli attributi ed in particolare il colore delle tre superfici

Visualizza / Linee Poligoni 3D vettoriali Nuovo aggiungere il layer edificato_vol

Importare come layer GRASS, con v.in.ogr il seguente shp: C:\EsercitazioniQGIS2013\Dati\DatiElettrodotti\shp\grigliatononregolare_punti_B.shp chiamandolo: grigliato3D ed avendo cura di importarlo come layer 3D (si veda la seguente immagine

Impostare la regione di lavoro in modo che copra i punti importati Da notare che adesso i punti non sono disposti secondo una griglia regolare.

Nella shell, per controllare la regione corrente g.region -p -3

nella shell: v.vol.rst input=grigliato wcolumn="B" elev=B3dinterp --overwrite

Toolbox Grass / Albero moduli /visualizzazione 3D / Nviz

Importare in grass un layer che rappresenta la superficie sulla quale si calcolano 10 microtesla (da shell): v.in.ogr -o "dsn=C:/EsercitazioniQGIS2013/Dati/RisultatiCalcolo/3D/04a_anello_1_tubo.shp" output=tubo10uT -z snap=-1 min_area=0.0001 -r -c

Visualizzare in NVIZ le superfici importate e verificare se ci sono disallineamenti. (Approfondire origine disallineamenti in prossime edizioni).

Per generare un raster che rappresenta il campo a 10m e a 20 m dal suolo, tenendo conto dell'andamento del terreno: shell: r.mapcalc 10mdalsuolo=dtm_10m_12pt+10 shell: r.mapcalc 10mdalsuolo=dtm_20m_12pt+10 shell: r3.cross.rast utilizzare l'interfaccia grafica del comando in Grass come nelle seguenti immagini

Modificare la mappa di colore del raster: