Tutorial esercitazioni inquinamento elettromagnetico a bassa frequenza Aggiornato al 4 maggio 2012

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)
• Importazione e primo filtraggio edificato
• 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
• 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
  • 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
• 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 edizioni 2007/2008, 2008/2009 del master fino a giungere alla forma attuale dopo le ultime modifiche effettuate prima dell'edizione 2011/2012. Le esercitazioni seguono 4 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.
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 Opne 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.
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
• Wiki QGIS: indice dei titoli
• Forum di QGIS
• Mailing lists

• 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

• GRASS per windows

Per installare la versione di QGIS usata nelle esercitazioni del 2012 sul proprio computer: Il file zip osgeopkgs2012.zip contiene qgis 1.7 e grass. Queste sono le versioni del software usate nel 2012. Utilizzando questo pacchetto è possibile installare qgis anche su computer non connessi ad internet. Per fare questo è necessario: scompattare il pacchetto osgeopkgs2012.zip in c: usando l'opzione "use folder names" lanciare il file setup17.bat

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

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ù Impostazioni / Proprietà del progetto /

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

menù Impostazioni / Proprietà del progetto / scheda Proiezione 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:/EsercitazioniQGIS2012/Qgis_proj dove salvare tutti i progetti di QGIS

menù Plugin / Gestione Plugins

selezionare GRASS e confermare (OK)

menù Plugin / Grass / New Mapset

Creare una nuova cartella "DatasetGrass" (ad esempio in C:/EsercitazioniQGIS2012/)

Selezionare la cartella appena creata e confermare

Next

Creare una nuova location "LocationToscanaGaussBoaga" 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)<\br> Il sistema Dataset/Location/Mapset corrisponde ad un sistema di organizzazione dei file su tre livelli.

menù Layer / aggiungi un raster

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

menù Visualizza / vista massima

Se non si dovesse vedere nulla fare doppio clic sul nome del raster nella legenda Selezionare la scheda Generale

Cambiare il sistema di riferimento con quello generato precedentemente. NB: i raster sono immediatamente visualizzati se non si abilita la riproiezione al volo nelle proprietà del progetto.

menù Layer / aggiungi tutto alla vista panoramica

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

menù Visualizza / nuovo segnalibro : "Argingrosso"

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

menù Plugins / Grass / apri strumenti Grass

Toolbox GRASS /Modules /File / Import / Import raster /r.in.GDAL

Esegui

Chiudi

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

Selezionare 275030 nella legenda e Zoom al layer

Selezionare 275030 (il .tif) nella legenda e Elimina Layer

Aggiungi raster Grass / 275030

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

Salvare il progetto menù file / salva

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

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

Verranno visualizzate correttamente solo le curve di livello. Questo perché solo delle shp delle curve ha associate le informazioni relative alla proiezione e quindi qgis/grass riesce a riproiettare "al volo" il layer.

Per visualizzare correttamente anche i punti quotati si deve associargli la proiezione, come si è fatto precedentemente con i raster. Tasto destro sulla legenda -> proprietà -> scheda Generale

Oppure Tasto destro sulla legenda -> imposta il SR del layer

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)

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

selezionare nella legenda 275030_CurveLivello Toolbox GRASS /Modules /File / Import / Import vector /v.in.ogr

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 (uno alla volta) layer /elimina layer

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

Toolbox GRASS /Modules /File / Map type conversion / convert a raster to vector / v.to.points Nella versione 1.4 Toolbox GRASS /Modules / vector / geometry management /manage features / v.to.points Nella versione 1.7 Toolbox GRASS /Albero moduli / 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 / vector / spatial analysis / overlay /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 / Modules / raster / surface management / interpolate surface v.surf.idw chiamo il raster dtm_10m_12pt

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).

Con QGIS1.4, se si vogliono visualizzare correttamente i raster appena generati si deve modificarne le propreià, come fatto per il primo (doppio clic sulla legenda-> proprietà ->Scheda Simbologia -> Stira MinMax -> Apply) Ciò non è più necessario con QGIS1.7

Aprire la scheda scheda metadati di uno dei dtm creati 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:/EsecitazioniQGIS2010/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

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". ad esempio è possibile generare la cosiddetta mappa shaded cioè una visualizzazione del modello del terreno in cui compaiono le ombre che ne evidenziano l'orografia. Toolbox Grass / Raster / Spatial analysis /

Si esegue l'analisi per le pendenze (r.slope.aspect.slope)

Si esegue l'analisi per le pendenze (r.slope.aspect.aspect)

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 shaded, quella delle pendenze e quella degli "aspect").

menù Layer / nuovo /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

visualizza / mostra segnalibri / Argingrosso / zoom a seguire il tracciato dell'elettrodotto prima verso est e poi verso sud fino alla sottostazione presso "il ronco".

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 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

è possibile cambiare colore e font al testo delle etichette.

caricare il raster GRASS dtm_10m_12pt creato in precedenza

Toolbox GRASS Modules / Vectors / Develop map / convert 2D vector to 3D by sampling raster

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 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 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

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:/EsercitazioniQGIS2012/WORK/s407_point.txt fs="\t"

aprire con wordpad il file di testo appena creato 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:/EsercitazioniQGIS2012/WORK/s407_std.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\OSgeo4w\apps\grass\grass-6.4.2\scripts Copiare il file v.to.lines.bat ed incollarlo nella cartella C:\Programmi\OSgeo4w\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

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 Quick card reference awk

Se poi si vuole aggiungere ala 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\OSgeo4w\apps\qgis\grass\modules modificare il file C:\Programmi\OSGEO\apps\qgis\grass\config/default.qgc agiungendo 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

Nuovo mapset :"MasterGIS" 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 :"MasterGIS" 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 :"MasterGIS" Next

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

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

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

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:\EsercitazioniQGIS2012\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.

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.

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

dalla shell di GRASS r.mapcalc surface=1

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

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 "disegna". 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: ptgriglia5m_buff300m

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=ptgriglia5m_buff300m type=point rast=dtm_090_10m method=bilinear output=ptgriglia5m3D_buff300m

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

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

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:/EsercitazioniQGIS2012/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:\EsercitazioniQGIS2012\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.

Non si ha il diretto controllo sul colore associato al valore di campo.

Simile a "simbolo graduato" ma in questo caso la classe è costituita dai punti che assumono lo stesso valore del campo di classificazione. Modalità di visualizzazione adatta a casi in cui ci siano pochi valori del campo di classificazione che si ripetono per molti punti.

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

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

importazione grigliato a 7 m

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 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

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 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

generare il nuovo raster 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(edificato_7m==1&&B_7m>3,1,null())" r.mapcalc "ed_7m_1to3uT=if(edificato_7m==1&&B_7m>1&&B_7m<3,1,null())" r.mapcalc "ed_7m_04to1uT=if(edificato_7m==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

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

Importare come layer GRASS, con v.in.ogr i seguenti shp: C:\EsercitazioniQGIS2012\Dati\DatiElettrodotti\shp\sostegni_40700_strutturasostegni.shp C:\EsercitazioniQGIS2012\Dati\DatiElettrodotti\shp\catenarie_40700_catenarie.shp C:\EsercitazioniQGIS2012\Dati\DatiElettrodotti\shp\sostegni_40900_strutturasostegni.shp C:\EsercitazioniQGIS2012\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:\EsercitazioniQGIS2012\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:/EsercitazioniQGIS2012/work Creare nella cartella work un file seriecomandi.txt di testo incollandoci il seguente contenuto: v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B24,0_punti_B.shp output=B24 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B22,0_punti_B.shp output=B22 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B20,0_punti_B.shp output=B20 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B18,0_punti_B.shp output=B18 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B16,0_punti_B.shp output=B16 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B14,0_punti_B.shp output=B14 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B12,0_punti_B.shp output=B12 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B10,0_punti_B.shp output=B10 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B8,0_punti_B.shp output=B08 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B6,0_punti_B.shp output=B06 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B4,0_punti_B.shp output=B04 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/Dati/RisultatiCalcolo/3D/02B2,0_punti_B.shp output=B02 -z snap=-1 min_area=0.0001 --overwrite v.in.ogr -o dsn=C:/EsercitazioniQGIS2012/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:\EsercitazioniQGIS2012\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:/EsercitazioniQGIS2012/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).