QGIS: confronto tra calcolatori raster

In questa pagina sono confrontate alcune funzioni che permettono di compiere operazioni matematiche/logiche su dati raster in QGIS.
Tali funzioni/comandi sono:
- Calcolatore raster integrato in QGIS
- Raster calculator di GDAL
- r.mapcalculator di GRASS
- Raster calculator di SAGA
- Band Math di Orfeo Toolbox
A parte il calcolatore raster integrato in QGIS, gli altri sono inclusi negli Strumenti di Processing per i quali occorre installare le relative applicazioni e abilitare (ed eventualmente settare) il plugin "Processing".

Il calcolatore raster di QGIS è facilmente apribile dal menù "Raster > Calcolatore raster...". Esso ci mostra l'elenco dei layer raster (e delle loro bande) caricati in QGIS con la dicitura nome_layer@n°_banda; in tal modo possiamo includere nel calcolo le sole bande che ci interessano. Tra le opzioni del "Risultato del layer", oltre a percorso/nome del file di output, possiamo scegliere il formato (default: GeoTIFF), cambiere l'estensione dell'area sulla quale verrà eseguito il calcolo (e la dimensione in pixel del file risultante) e cambiare il sistema di riferimento.
Al termine del calcolo viente prodotto un raster con data type Float32.
La cosa interessante di questo calcolatore raster è che permette di operare su bande con estensioni e risoluzioni spaziali diverse (ma nello stesso sistema di riferimento). L'estensione del risultato sarà quella specificata (X min, X max, Y min, Y max) mentre la risoluzione verrà calcolata automaticamente per adeguarsi al numero di righe e colonne che impostiamo; personalmente avrei preferito la possibilità di fare il contrario, ovvero specificare la risoluzione desiderata in modo che il programma producesse un numero adeguato di righe e colonne. Il cambiamento della risoluzione richiede un ricampionamento; purtroppo non è possibile scegliere il metodo di ricampionamento che è fissato a nearest neighbor (vicino più prossimo).
Sr in output: questa opzione NON permette di trasformare l'output nel CRS desiderato ma serve solo per specificare a quale sistema di riferimento appartengono le coordinate sopra inserite.

Strumenti di Processing: GDAL/OGR > [GDAL] Miscellaneous > Raster calculator.
Il calcolatore raster incluso negli strumenti di GDAL richiama il comando gdal_calc.py (vedi www.gdal.org/gdal_calc.html).
La finestra di dialogo permette di specificare fino a 6 raster da inserire nel calcolo (i raster devono avere la stessa dimensione in pixel). Nel caso di file multibanda è possibile scegliere la banda da considerare.
Tra le opzioni del file in uscita è possibile indicare il data type (di defalut è Float32, ovvero virgola mobile 32bit/4byte che permette la memorizzazione di numeri decimali) e il valore da considerare nullo (NoDataValue).
Inoltre è possibile aggiungere altri parametri aggiuntivi del comando gdal_calc.py ovvero:

Strumenti di Processing: GRASS > Raster (r.*) > r.mapcalculator.
(> https://grass.osgeo.org/grass64/manuals/r.mapcalculator.html).
Permette di inserire fino a 6 raster in uno stesso calcolo; purtroppo non si possono scegliere le bande.
Le opzioni relative alla regione di GRASS permettono di specificare l'estensione geografica e la dimensione del pixel su cui effettuare il calcolo (quindi estensione e risoluzione del raster risultante).
ATTENZIONE: i raster di input devono avere lo stesso data type il quale sarà usato anche per il risultato finale ed i calcoli intermedi. Se la formula contiene divisioni o altri operatori che possono dare risultati decimali è raccomandabile convertire prima i raster di input a float32 (numero in virgola mobile di 32 bit).
Nella figura sottostante è riportata la formula d'esempio "(A-B)/(A+B)*1000". Essa fornisce risultati sbagliati se i raster A e B sono in formato di numero intero (data type: integer) poichè la parte iniziale "(A-B)/(A+B)" può fornire solo risultati interi -1, 0 o +1 che poi vengono moltiplicati per 1000. Quindi occorre convertire prima i raster in virgola mobile (float) oppure cambiare la formula in "(A-B)*1000/(A+B)".

Strumenti di Processing: SAGA > Grid - Calculus > Raster calculator.
Nelle versioni di SAGA dalla 2.1.3 in poi (non ancora incluse in OSGeo4W) è chiamato "Grid calculator".
(> www.saga-gis.org/saga_module_doc/2.1.3/grid_calculus_1.html).
Come molti altri moduli di SAGA io calcolatore accetta solo raster con un'unica banda.
Nella prima riga della finestra di dialogo si sceglie il raster principale, con la seconda si possono indicare i raster aggiuntivi. Non sembra esserci un limite al numero di raster.
Si può operare anche su raster con dimensioni (estensione e risoluzione) diverse; il risultato prenderà le proprietà del raster principale indicato nella prima riga, ma avrà valore nullo nei punti non coperti anche dai raster aggiuntivi.
Infine è possibile specificare il data type. ATTENZIONE: la scelta del data type influenza non solo il risultato finale ma anche le operazioni di calcolo. Ad esempio impostando a integer, anche i valori dei calcoli intermedi vengono arrotondati a numeri interi per cui il risultato della formula "(a-b)/(a+b)*1000" è diverso da "(a-b)*1000/(a+b)". Se la formula contiene divisioni o altri operatori che possono dare risultati decimali è raccomandabile impostare il data type su floating point. Per calcoli con operazioni intermedie che prevedono valori assoluti molto elevati (superiori al limite dei 2 byte) occorre usare data type di 4 o 8 byte.

Strumenti di Processing: Orfeo Toolbox > Miscellaneous > Band Math.
(> www.orfeo-toolbox.org//Applications/BandMath.html).
Le bande raster su cui operare sono indicate con imXbY dove X indica il numero del layer raster e Y il numero della banda.
I raster devono avere la stessa dimensione in pixel.
Data type fisso a float32 (numero in virgola mobile di 32 bit).

La tabella sottostante riporta la presenza o meno di alcune opzioni nei diversi calcolatori; sono state considerate:
- scelta bande: possibilità di usare raster multibanda scegliendo la banda d'interesse;
- limite numero di raster: eventuale limite al numero di raster che possono essere inclusi nel calcolo;
- scelta data type: possibilità di stabilire il data type del raster risultante;
- modifica area di calcolo: possibilità di impostare estensione e risoluzione dell'area di calcolo e quindi del raster risultante;
- tempo ndvi: tempo impiegato nel calcolare un indice ndvi su bande Landsat 8 con la formula "(NIR-RED)/(NIR+RED)*1000" (ovviamente il tempo cambia in base alle prestazioni dell'hardware; tali dati sono utili come informazioni sulla velocità relativa).

Da ricordare che i calcolatori di GDAL, GRASS, SAGA e Orfeo Toolbox, come tutti gli strumenti di processing, permettono memorizzare il raster risultante come file temporaneo che in seguito può essere salvato con formato e impostazioni definite dall'utente.