Universita` degli Studi della BasilicataDipartimento di Ingegneria e Fisica dell'Ambiente (DIFA)

GRASS GIS: Primi passi e operazioni sui raster

Ing. Margherita Di Leo

Avvio del programma

● Avviare il programma dal launcher;

Avvio del programma

● Avviare il programma dal launcher;

● Sfogliare la cartella grass data;

Avvio del programma

● Avviare il programma dal launcher;

● Sfogliare la cartella grass data;

● Scegliere la location North Carolina;

Avvio del programma

● Avviare il programma dal launcher;

● Sfogliare la cartella grass data;

● Scegliere la location North Carolina;

● Scegliere il mapset user1;

Avvio del programma

● Avviare il programma dal launcher;

● Sfogliare la cartella grass data;

● Scegliere la location North Carolina;

● Scegliere il mapset user1;

● Avvio.

Map Display

Layer Manager

TerminalMap Layers

Command console

Search module

Interfaccia utente (GUI)

Python shell

● La prima volta che si usa GRASS e` importante impostare la regione computazionale.

● Tutte le elaborazioni successive vengono effettuate sulla regione computazionale prescelta.

● La regione (estensione e risoluzione) si imposta con il comando g.region.

Impostare la regione

Impostare la regione

●Digitando g.region nella command console e premendo invio, si manifesta l'interfaccia grafica del comando.

●Nel tab Existing e` possibile impostare la regione sulla base di una regione esistente salvata precedentemente oppure su un file esistente.

Impostare la regione

●Nel tab Bounds e` possibile impostare la regione indicando i limiti precisi.

Impostare la regione

●Nel tab Resolution e` possibile impostare la risoluzione della cella.

Impostare la regione

●Nel tab Effects e` possibile salvare la regione corrente (regione di visualizzazione) con un nome, in modo da poterla successivamente richiamare.

Impostare la regione

●Nel tab Print e` possibile spuntare le varie opzioni di cio` che desideriamo venga stampato a schermo.

Impostare la regione

g.region -ap rast=elevation@PERMANENT

●Per il nostro progetto, impostiamo come regione computazionale la risoluzione ed estensione del raster elevation.

Visualizzare il contenuto del mapset● Vogliamo vedere cosa e` contenuto nel

dataset di esempio. Il comando che si usa allo scopo e` g.list.

● La sintassi e`: g.list type=tipo1[,tipo2,...]

● Ad esempio:g.list type=rast,vect

elenchera` tutti i raster e i vettoriali presenti nel mapset, mentre:g.list type=rast3d

dira` che non e` presente alcun raster3d.

Visualizzare un rasterAdd raster map layer

d.rastd.rast –helpd.rast map=elevation@PERMANENT

Visualizzare un rasterZoom to selected map

Informazioni su una mappa raster

r.infor.info –helpr.info rast=elevation@PERMANENT

r.infor.info –helpr.info rast=elevation@PERMANENT

r.info da informazioni sulla mappa raster: titolo,

tipo di dato, righe, colonne, numero di celle, proiezione, estensione, risoluzione, intervallo dei valori, informazioni

sulla provenienza (come e` stata generata e da

chi).

Informazioni su una mappa raster

Interrogare una mappa rasterQuery raster/vector map(s)

r.what restituisce il valore della mappa nel punto selezionato.

Utile per individuare le coordinate della sezione di chiusura!

Interrogare piu` mappe raster● r.what puo` essere utilizzato da display oppure da linea di

comando, come in generale tutti i comandi di GRASS.

● Utilizzarlo da linea di comando e` utile quando, note le coordinate, vogliamo conoscere i valori di una o piu` mappe raster in quel dato punto.

● Per fare un esempio, consideriamo la mappa raster landclass96.

● Interroghiamo con il comando r.what entrambe le mappe contemporaneamente, per conoscere il valore di entrambe le mappe nella posizione considerata:

r.what input=elevation,landclass96 east_north=636889.830508,222716.101695

● La risposta e`:

636889.830508|222716.101695||98.26356|5

ovvero: east | north || elevation | landclass96

Suggerimenti

● Per interrogare la mappa landclass96, utilizzando la linea di comando, non e` necessario caricare la mappa nel layer tree!

● L'uso della linea di comando velocizza il lavoro;

● Per dubbi sulla sintassi dei comandi da riga di comando e` sempre possibile farli seguire da: --help

Generare report e statistiche

r.report fornisce le seguenti informazioni:

● Location in cui stiamo lavorando;

● Limiti n-s e-o della mappa;

● Risoluzione della mappa n-s e-o;

● Eventuale presenza della mask;

● Descrizione delle categorie presenti nella mappa;

● Estensione delle varie categorie nell'unita` di misura scelta dall'utente.

Generare report e statistiche

r.report map=name[,name,...][units=string[,string,...]]

[output=name][nsteps=value]

● Units: Le unita` possono essere: mi,me,k,a,h,c,p

(mi: miglia, me: metri, k: chilometri, a: acri, h: ettari, c: numero di celle, p: percentuale);

● Output: Nome del file di output;

● Nsteps: Numero di intervalli (default: 255). Ha senso per mappe con valori continui, come elevazione, pendenza, ecc.

Generare report e statistiche

r.report map=geology_30m units=k

Generare report e statisticher.univar serve per riportare le statistiche univariate dell'intervallo di dati della mappa:

● Numero di celle;

● Minimo e massimo valore delle celle;

● Intervallo;

● Media aritmetica;

● Varianza della popolazione;

● Deviazione standard

● Coefficiente di variazione.

● Con la flag -e vengono calcolati anche: primo quartile, secondo quartile (mediana), terzo quartile e percentile dato.

Generare report e statistiche

r.univar map=slope

Riclassificazioner.reclass serve per aggregare delle categorie a partire da una mappa raster. Consideriamo ad esempio la mappa landclass96. Utilizzando r.report su tale mappa otteniamo le categorie riportate a lato.

Procediamo alla riclassificazione creando un file di testo sul nostro Desktop (utilizziamo WordPad).

Il file di testo, denominato “regole.txt”, conterra` le regole per la riclassificazione.

r.report map=landclass96

RiclassificazioneSupponiamo di voler aggregare tutte le categorie che rappresentano la copertura vegetale non di origine agricola.

Per costruire il file delle regole:

vecchio_valore = nuovo_valore etichetta

Il file regole.txt sara` fatto cosi`:

1 = 1 developed

2 = 2 agriculture

3 = 3 vegetazione

4 = 3 vegetazione

5 = 3 vegetazione

6 = 4 water

7 = 5 sediment

Salviamo il file regole.txt sul Desktop.

RiclassificazioneOra siamo pronti a riclassificare la mappa landclass96.

Digitiamo nel command console r.reclass e diamo invio. Si aprira` l'interfaccia grafica.

Inseriamo in input: landclass96 e in output landclass_riclassificata.

RiclassificazioneNel tab “Optional” sfogliamo il percorso al file “regole.txt” che abbiamo salvato sul Desktop.

RiclassificazioneIl risultato e` una mappa raster riclassificata:

Prima: landclass96 Dopo: landclass_riclassificata

Riclassificazione

IMPORTANTE: La riclassificazione NON crea una nuova mappa, ma crea soltanto un file di riclassificazione

che viene conservato in una cartella del mapset di GRASS. Per generare un file raster riclassificato, e` necessario ricorrere alla map algebra. Nel terminale, digitiamo:

r.mapcalc “nuova_landclass = landclass_riclassificata”

● NOTA: Osservando bene “landclass_riclassificata” e “nuova_landclass”, si puo` notare che hanno risoluzione

diversa. Perche`? (Consiglio per la risposta: usare r.info e ricordare cio` che e` stato detto riguardo la regione

computazionale e g.region!)

Esercitazione

Supponiamo di voler creare una mappa di localizzazione probabile di reperti archeologici appartenenti ad una antica popolazione.

● Si ipotizza che tale popolazione costruisse i propri villaggi ad una quota inferiore a 100 m s.l.m., ad una distanza entro i 200 m dai corsi d'acqua, in territori con esposizione da sud-est a sud-ovest e pendenza compresa tra 0 e 25 gradi.

● Solo i terreni con superficie superiore a 1 ettaro saranno poi selezionati, poiché le aree con superfici inferiori non possono permettere lo sviluppo di un villaggio.

Creazione di un buffer

Ricaviamo prima di tutto la mappa dei punti che si trovano ad una distanza entro i 200 m dai corsi d'acqua. Utilizziamo il comando r.buffer.

r.buffer input=streams_derived@PERMANENT output=fiumi_buffer distances=200 units=meters

Creazione di un buffer

Il risultato è una mappa in cui hanno valore 1 le celle contenenti i fiumi, valore 2 le celle distanti fino a 200 metri da essi e * le altre.

[* questo simbolo indica i null data, e si gestisce con il comando r.null]

r.buffer input=streams_derived@PERMANENT output=fiumi_buffer distances=200 units=meters

Map algebra

r.mapcalc "fino_a_100=if(elevation<=100,1,0)"

Vogliamo ricavare la mappa delle aree con quota non superiore a 100 m s.l.m. Nel TERMINALE inseriamo il comando:

Map algebra

r.mapcalc "fino_a_100=if(elevation<=100,1,0)"

Il risultato e` una mappa binaria, con valore 1 dove la quota e` minore o uguale a 100, e valore 0 altrove.

Map algebraCerchiamo ora le aree con esposizione a sud-est, sud e sud-ovest. La convenzione delle direzioni in GRASS e` la seguente:

I valori della mappa dell'esposizione (aspect) che ci interessano sono dunque: >=225 & <=315

Tradotto in map algebra, si scrive:

r.mapcalc "buona_esposizione=if(aspect>=225&&aspect<=315,1,0)"

Map algebra

Scriviamo nel terminale:

r.mapcalc "buona_esposizione=if(aspect>=225&&aspect<=315,1,0)"

Map algebra

r.mapcalc "buona_esposizione=if(aspect>=225&&aspect<=315,1,0)"

Mappa delle esposizioni sud-est, sud, sud-ovest.

Map algebraNOTA: la map algebra ha anche un'interfaccia grafica, per chi la preferisce.

Map algebraCalcoliamo ora la mappa delle aree con pendenza non superiore a 25 gradi, usiamo la mappa slope presente nel dataset, e impostiamo la condizione con la solita map algebra, scrivendo nel terminale:

r.mapcalc "buona_pendenza=if(slope<=25,1,0)"

Map algebra

Mappa delle aree con pendenza non superiore a 25 gradi.

r.mapcalc "buona_pendenza=if(slope<=25,1,0)"

Map algebraCalcoliamo infine la mappa di intersezione di tutte le condizioni richieste(si ricorda che l'intersezione corrisponde all'operatore AND (&), mentre l'unione corrisponde all'operatore OR(|)):

r.mapcalc "intersezione = if( ( fiumi_buffer == 2 && fino_a_100 == 1 && buona_esposizione == 1 && buona_pendenza == 1 ), 1, null() )"

Map algebraLa mappa di intersezione di tutte le condizioni richieste:

r.mapcalc "intersezione = if( ( fiumi_buffer == 2 && fino_a_100 == 1 && buona_esposizione == 1 && buona_pendenza == 1 ), 1, null() )"

Riclassificazione in base all'areaCome ultimo step, dobbiamo eliminare tutte le aree che risultano inferiori ad un ettaro. Per fare cio` utilizziamo il comando r.reclass.area. Digitiamolo della command console, e nell'interfaccia grafica inseriamo i parametri richiesti:

r.reclass.area input=intersezione output=intersezione_max_ettaro greater=1

Riclassificazione in base all'area

Nel tab “Optional” inseriamo la condizione che le aree debbano essere superiori ad 1 ettaro:

r.reclass.area input=intersezione output=intersezione_max_ettaro greater=1

Risultato

Et voila`!

La mappa intersezione_max_ettaro

rappresenta il risultato della nostra esercitazione.

License of this document

             This work is licensed under a Creative Commons License.                         http://creativecommons.org/licenses/by­sa/3.0/

                                      © 2011, Margherita Di Leo, Italy                                      dileomargherita AT gmail DOT com

License details: Attribution­ShareAlike 3.0:You are free:  * to Share — to copy, distribute and transmit the work  * to Remix — to adapt the work Under the following conditions:  * Attribution — You must attribute the work in the manner specified by the author or licensor (but not in any way that suggests that they endorse you or your use of the work).  *Share Alike — If you alter, transform, or build upon this work, you may distribute the resulting work only under the same, similar or a compatible license.

With the understanding that:  * Waiver — Any of the above conditions can be waived if you get permission from the copyright holder.  * Other Rights — In no way are any of the following rights affected by the license:       o Your fair dealing or fair use rights;       o The author's moral rights;       o Rights other persons may have either in the work itself or in how the work is used, such as publicity or privacy rights.