Questo tutorial illustra la classificazione della copertura del suolo di vari immagini Landsat per creare una classificazione di un'area vasta utilizzando il Semi-Automatic Classification Plugin (SCP). Per i tutorial base vedi Tutorial 1: Your First Land Cover Classification e Tutorial 2: Land Cover Classification of Landsat Images.
L'area di studio di questo tutorial è il Costa Rica , un Paese del Centro America che ha un'estensione di circa 51,000 chilometri quadrati. In particolare, andremo a classificare immagini Landsat 8 e Landsat 8, mascherando le nuvole e creando un mosaico delle classificazioni. Identificheremo le seguenti classi di copertura del suolo;
- Costruito (Built-up);
- Vegetazione (Vegetation);
- Suolo (Soil);
- Acqua (Water).
Di seguito il video del tutorial. Si raccomanda di seguire anche i passaggi descritti di seguito.
Installazione del Plugin
Per prima cosa, installiamo il SCP. Per informazioni sull'installazione di QGIS in vari sistemi vedere Plugin Installation.
Aprire QGIS; dal menu principale, selezionare
Plugins
> Manage and Install Plugins
;
Dal menu
All
, selezionare Semi-Automatic Classification Plugin e click sul bottone Install plugin
;
Il SCP dovrebbe essere attivato automaticamente; comunque, verificare che Semi-Automatic Classification Plugin sia attivo nel menu
Installed
(potrebbe essere necessario riavviare QGIS);
Andremo ad impostare alcune opzioni di SCP per ottimizzare i processi seguenti.
Aprire Settings con un click sul bottone e selezionare Settings: Processing. Ora, in Classification process attivare le opzioni
Use virtual rasters for temp files
e Raster compression
, per salvare spazio su disco.
In RAM impostare la RAM disponibile (in MB). Inserire metà della memoria RAM del sistema; per esempio se il sistema ha 2GB di RAM inserire 1024.
Per facilitare la fotointerpretazione nei passaggi successivi, si consiglia l'uso del plugin OpenLayers che permette la visualizzazione di varie mappe in QGIS.
Download e Pre processing di immagini Landsat
Scaricheremo immagini Landsat 7 e 8 usando lo strumento Download Landsat. Le immagini Landsat sono disponibili dal U.S. Geological Survey, e queste bande sono scaricate tramite Google Earth Engine Amazon Web Services. Inoltre, convertiremo le immagini in riflettanza, e applicheremo la correzione atmosferica DOS1 (vedi Landsat image conversion to reflectance and DOS1 atmospheric correction).
Per prima cosa, scarichiamo il database Landsat in SCP. Apri il tab Download Landsat e click sul bottone nel SCP menu o Toolbar.
Click sul bottone
Select database directory
per definire dove salvare il database. Creare una nuova cartella (es. LandsatDB
) nel percorso utente. Click sul bottone Update database
e click Yes
nella domanda seguente.
Ora possiamo definire le coordinate dell'area di studio, click
Find images
e cercare le immagini Landsat. Ogni immagine Landsat ha un ID unico (identificativo). In questo tutorial useremo due immagini Landsat 8 acquisite nel Febbraio 2014 (ID: LC80150532014050LGN00 e LC80160532014057LGN00) e un'immagine Landsat 7 acquisita nel Marzo 2014 (ID: LE70150532014090EDC00); ovviamente, ulteriori immagini sono necessarie per la classificazione dell'intero Paese.
Con SCP è possibile trovare un'immagine tramite l'ID usando le opzioni di Ricerca. In particolare, incollare in
Image ID
i seguenti ID e click su Find images
:
Dopo alcuni secondi, le tre immagini saranno elencate in Landsat images.
Prima di scaricare le immagini, dobbiamo definire le opzioni di conversione in riflettanza che sarà eseguita automaticamente per le immagini scaricate.
Aprire il tab Landsat con un click sul bottone nella Toolbar ; abilitare le opzioni
Apply DOS1 atmospheric correction
e Brightness temperature in Celsius
. Inoltre, lasciare abilitato Use NoData value (image has blackborder)
.
Ora, aprire il tab Download Landsat e deselezionare le opzioni
only if preview in Layers
e Load bands in QGIS
(lasciare selezionato Pre process images
per convertire le bande in riflettanza automaticamente). Click sul bottone Download images from list
e selezionare la cartella di output ed iniziare il processo di download (può richiedere del tempo).
Quando il download è completo, varie cartelle saranno create nella cartella di output, con il nome uguale all'ID Landsat, in cui sono contenute le bande originali e le bande convertite (con il suffisso
_converted
).Classificazione delle immagini Landsat
Iniziamo la classificazione dell'immagine Landsat 8
LC80150532014050LGN00
convertita in riflettanza. Aprire la cartella LC80150532014050LGN00_converted
.
In QGIS, aprire le bande seguenti (anche con drag and drop):
- RT_LC80150532014050LGN00_B2.tif = Blue;
- RT_LC80150532014050LGN00_B3.tif = Green;
- RT_LC80150532014050LGN00_B4.tif = Red;
- RT_LC80150532014050LGN00_B5.tif = Near-Infrared;
- RT_LC80150532014050LGN00_B6.tif = Short Wavelength Infrared 1;
- RT_LC80150532014050LGN00_B7.tif = Short Wavelength Infrared 2.
Aprire il tab Band set e click sul bottone nel SCP menu o la Toolbar. Click sul bottone
Select All
, Add rasters to set
, e Sort by name
per ordinare le bande in automatico. Infine, selezionare Landsat 8 OLI
dalla combo box Quick wavelength settings
, per impostare in automatico il centro della lunghezza d'onda per ogni banda.
Nella lista
RGB
selezionare l'elemento 4-3-2
per visualizzare il Color Composite con Near-Infrared, Red, e Green.
Ora dobbiamo creare il
Training shapefile
e Signature list file
per creare le Training Areas (ROIs) e calcolare le Firme Spettrali (per definizioni base vedere Tutorial 1: Your First Land Cover Classification).
Nella ROI Creation dock click sul bottone
New shp
e definire un nome (es. ROI.shp
) per la creazione di un Training shapefile
che conterrà i poligoni delle ROI. Lo shapefile create è aggiunto a QGIS ed il nome è visualizzato in Training shapefile .
Inoltre, click sul bottone
Save
nella Classification dock e definire un nome (es. SIG.xml
) per creare un Signature list file
che conterrà le firme spettrali. Il percorso del Signature list file
è mostrato in Signature list file .
Adesso siamo pronti per la creazione delle ROI. Useremo gli stessi codici per le ROI in tutte le immagini Landsat, come riportato nella tabella seguente.
Macroclass name | Macroclass ID |
---|---|
Built-up | 1 |
Vegetation | 2 |
Soil | 3 |
Water | 4 |
Per le nozioni base su come creare ROI vedere Create the ROIs. Si possono creare ROI disegnano manualmente un poligono usando il bottone o con region growing (segmentazione) premendo il bottone + e cliccando sulla mappa. Usa il bottone in ROI creation per zoomare sul poligono della ROI, e
Show
per mostrare e nascondere le ROI temporanee.
Creare quante più ROI possibile tramite la ROI Creation dock ed assegnare un univoco Class ID (C_ID) ad ogni ROI, e un Macroclass ID (MC_ID) della Macroclasse corrispondente. Se
Displaycursor for
è selezionato in ROI creation, il valore di NDVI del pixel sotto al cursore sarà mostrato nella mappa: utile per identificare la vegetazione (caratterizzata da alti valori di NDVI).
Le ROI sono usate per il calcolo delle firme spettrali che sono utilizzate dagli algoritmi di classificazione per classificare l'intera immagine. In questo tutorial useremo l'algoritmo Maximum Likelihood.
Dopo la creazione di ogni ROI è utile verificare la Distanza Spettrale per verificare la separabilità delle ROI: infatti, ogni ROI dovrebbe essere differente dalle altre, per evitare confusione spettrale ed ottenere risultati migliori nella classificazione.
Nella Signature list evidenziare le ROI e click sul button . Le firme spettrali sono aggiunge al Spectral Signature Plot.
Dopo la creazione di ogni ROI è utile verificare la Distanza Spettrale per verificare la separabilità delle ROI: infatti, ogni ROI dovrebbe essere differente dalle altre, per evitare confusione spettrale ed ottenere risultati migliori nella classificazione.
Nella Signature list evidenziare le ROI e click sul button . Le firme spettrali sono aggiunge al Spectral Signature Plot.
Ora click nel tab Spectral distances. Ogni tabella rappresenta la distanza spettrale di ogni combinazione di ROI.
Come mostrato nella figura seguente, il confronto tra Built-up e Soil mostra un Angolo Spettrale molto basso; questo significa alta similarità se usassimo gli algoritmi Spectra Angle Mapping o Minimum Distance. Il valore Jeffries-Matusita Distance è prossimo a 2; questo significa che le ROI sono separabili per l'algoritmo Maximum Likelihood.
Come mostrato nella figura seguente, il confronto tra Built-up e Soil mostra un Angolo Spettrale molto basso; questo significa alta similarità se usassimo gli algoritmi Spectra Angle Mapping o Minimum Distance. Il valore Jeffries-Matusita Distance è prossimo a 2; questo significa che le ROI sono separabili per l'algoritmo Maximum Likelihood.
Poiché useremo l'algoritmo Maximum Likelihood, è importante che Jeffries-Matusita Distance sia vicino a 2 per ogni combinazione di ROI.
Ora possiamo creare una anteprima di classificazione (vedi Create a Classification Preview per le nozioni base).
In Classification algorithm selezionare l'algoritmo di classificazione
Maximum Likelihood
. In Classification preview impostare Size
= 500 , click sul bottone +
e poi click sinistro nella mappa per creare un'anteprima di classificazione. Usare lo strumento Transparency
per mostrare in trasparenza la classificazione sull'immagine.
In Classification algorithm, click sul bottone
Come mostrato nella figura seguente, il raster algoritmo ha una simbologia a scala di grigi, dove aree scure rappresentano pixel che l'algoritmo ha trovato distanti da tutte le firme spettrali, e le aree chiare rappresentano pixel che sono molto simili alle firme spettrali.
In queste aree scure abbiamo un maggiore livello di incertezza. Quindi dobbiamo creare nuove ROI per migliorare l'accuratezza.
+
e poi click destro nella mappa per calcolare il raster algoritmo
. Il raster algoritmo rappresenta il risultato del calcolo degli Algoritmi di Classificazione; è utile per localizzare dove dobbiamo creare nuove ROI.Come mostrato nella figura seguente, il raster algoritmo ha una simbologia a scala di grigi, dove aree scure rappresentano pixel che l'algoritmo ha trovato distanti da tutte le firme spettrali, e le aree chiare rappresentano pixel che sono molto simili alle firme spettrali.
In queste aree scure abbiamo un maggiore livello di incertezza. Quindi dobbiamo creare nuove ROI per migliorare l'accuratezza.
Possiamo notare la presenza di nuvole nell'immagine. Per evitare errori di classificazione dobbiamo mascherare le nuvole.
Ci sono vari metodi per mascherare le nuvole; durante la fase di classificazione, un semplice metodo consiste nella creazione di ROI. Crea una ROI in una nuvola nell'immagine, e assegna un unico Class ID e il Macroclass ID uguale a 0. Infatti, MC ID = 0 è utilizzato dal plugin per la categoria
Ci sono vari metodi per mascherare le nuvole; durante la fase di classificazione, un semplice metodo consiste nella creazione di ROI. Crea una ROI in una nuvola nell'immagine, e assegna un unico Class ID e il Macroclass ID uguale a 0. Infatti, MC ID = 0 è utilizzato dal plugin per la categoria
Unclassified (Non classificato)
, cioè i pixel delle nuvole saranno non classificati (mascherati).
Nell'immagine seguente possiamo vedere che le nuvole sono mascherate. Tuttavia i pixel vicino al bordo delle nuvole sono classificati non correttamente come costruito (Built-up).
Nel prossimo paragrafo saranno descritti metodi più efficaci per mascherare le nuvole dopo il processo di classificazione.
Quando il risultato delle anteprime è soddisfacente, possiamo eseguire la classificazione dell'intera area. In Classification algorithm, attiva il checkbox
Nel prossimo paragrafo saranno descritti metodi più efficaci per mascherare le nuvole dopo il processo di classificazione.
Quando il risultato delle anteprime è soddisfacente, possiamo eseguire la classificazione dell'intera area. In Classification algorithm, attiva il checkbox
Use Macroclass ID
. In Classification output click sul bottone Perform classification
e definisci il nome della classificazione di output (es. classification_1.tif
).
Possiamo vedere che parte delle nuvole sono nere (non classificate), comunque molti pixel delle nuvole sono classificati come costruito. Inoltre, il bordo dell'immagini è classificato come costruito.
Andremo a correggere questi errori e rifinire la classificazione nei prossimi passaggi.
Adesso, in QGIS apri le seguenti bande Landsat 8 che sono nella cartella
Andremo a correggere questi errori e rifinire la classificazione nei prossimi passaggi.
Adesso, in QGIS apri le seguenti bande Landsat 8 che sono nella cartella
LC80160532014057LGN00_converted
.
- RT_LC80160532014057LGN00_B2.tif = Blue;
- RT_LC80160532014057LGN00_B3.tif = Green;
- RT_LC80160532014057LGN00_B4.tif = Red;
- RT_LC80160532014057LGN00_B5.tif = Near-Infrared;
- RT_LC80160532014057LGN00_B6.tif = Short Wavelength Infrared 1;
- RT_LC80160532014057LGN00_B7.tif = Short Wavelength Infrared 2.
Ripetere i passaggi precedenti per la creazione del Band set, del
Training shapefile
e Signature list file
.
Creare una classificazione di copertura del suolo ripetendo i passaggi precedenti.
In un nuovo progetto di QGIS, aprire le bande del Landsat 7 nella cartella
In un nuovo progetto di QGIS, aprire le bande del Landsat 7 nella cartella
LE70150532014090EDC00_converted
:
- RT_LE70150532014090EDC00_B1.tif = Blue;
- RT_LE70150532014090EDC00_B2.tif = Green;
- RT_LE70150532014090EDC00_B3.tif = Red;
- RT_LE70150532014090EDC00_B4.tif = Near-Infrared;
- RT_LE70150532014090EDC00_B5.tif = Short Wavelength Infrared 1;
- RT_LE70150532014090EDC00_B7.tif = Short Wavelength Infrared 2.
Si può vedere che l'immagine copre la stessa area dell'immagine Landsat 8
LC80150532014050LGN00
. Infatti, useremo la classificazione di questa immagine Landsat 7 image per riempire i pixel non classificati (Unclassified) dell'immagine Landsat 8.
Ancora, creare una classificazione di copertura del suolo seguendo i passaggi precedenti.
Ora abbiamo 3 classificazioni di copertura del suolo che possiamo migliorare i vari modi.
Miglioramento delle Classificazioni Usando l'NDVI
Calcoleremo l’NDVI per migliorare la classificazione usando Band calc (vedi Tutorial: Using the tool Band calc). In particolare, i pixel in cui il valore di NDVI è superiore ad una certa soglia sarà classificato come vegetazione (codice 2). Al di sotto di questa soglia, la classificazione Maximum Likelihood rimarrà invariata.
Ovviamente, questo è un esempio di integrazione di dati ancillari; potremo usare altri dati come altri indici vegetazionali o i risultati di altre classificazioni (es. usando Spectral Angle Mapping).
Ovviamente, questo è un esempio di integrazione di dati ancillari; potremo usare altri dati come altri indici vegetazionali o i risultati di altre classificazioni (es. usando Spectral Angle Mapping).
Ora, in QGIS carica le bande dell’immagine Landsat 8 LC80150532014050LGN00 e la rispettiva classificazione di copertura del suolo. Apri Band calc e click sul bottone Refresh list. In Band calc, calcola l’NDVI copiando la Expression seguente:
In questo caso abbiamo usando una soglia di NDVI pari a 0.6. Comunque, il valore del limite deve essere scelto per ogni immagine, perché l’NDVI varia da immagine a immagine.
Ora eseguiamo lo stesso miglioramento per le altre classificazioni di copertura del suolo. Per l’immagine Landsat 8 LC80160532014057LGN00 calcolare l’NDVI con l’espressione seguente:
Ora eseguiamo lo stesso miglioramento per le altre classificazioni di copertura del suolo. Per l’immagine Landsat 8 LC80160532014057LGN00 calcolare l’NDVI con l’espressione seguente:
e l’espressione seguente per migliorare la classificazione:
Per l’immagine Landsat 7 LE70150532014090EDC00 calcolare l’NDVI con l’espressione seguente:
e l’espressione seguente per migliorare la classificazione:
Ora che la classificazione della vegetazione è stata migliorata per le tre immagini, andremo a mascherare le nuvole e i pixel di bordo per evitare errori di classificazione.
Maschera delle Nuvole
Le immagini Landsat 8 includono la band Quality Assessment (QA) che è utile per identificare le nuvole. I valori dei pixel delle bande QA sono codici che rappresentano combinazioni delle condizioni della superficie e dell’atmosfera. Questi valori indicano con elevata confidenza pixel affetti da cirri o nuvole (per la descrizione di questi codici vedere la tabella a http://landsat.usgs.gov/L8QualityAssessmentBand.php ).
In QGIS, aprire la banda LC80150532014050LGN00_BQA che è all'interno della cartella LC80150532014050LGN00 delle immagini Landsat scaricate e la classification_1_NDVI. Copiare l’ Expression seguente in Band calc:
Le nuvole sono quasi completamente mascherate (Unclassified); tuttavia, alcuni pixel sono ancora classificati come costruito (in rosso). Possiamo fare la stessa cosa per l'immagine
LC80160532014057LGN00
usando l'Expression in Band calc:L'immagine Landsat 7 non ha una banda QA. Un altro metodo per mascherare le nuvole usa la banda del blu e dell'infrarosso termico (convertito in temperatura), basandosi sul fatto che le nuvole sono generalmente più fredde del suolo ed hanno valori di riflettanza alti nella banda del blu.
Landsat 7 è anche affetto da righe nere (cioè SLC-off) che andremo a mascherare ugualmente. Creeremo un'espressione che identifichi i valori dei pixel al di sotto di una certa soglia di temperatura per la banda termica (banda 6 per il Landsat 7), e oltre una certa soglia di riflettanza per la banda del blu (banda 1).
In QGIS caricare tutte le bande Landsat dalla cartella
LE70150532014090EDC00_converted
. Usare la seguente espressione in Band calc:La prima parte (
("RT_LE70150532014090EDC00_B6_VCID_1"<23) & ("RT_LE70150532014090EDC00_B1">0.1)
) significa che maschereremo i pixel con temperatura inferiore a 23°C e riflettanza della banda blu maggiore di 0.1 . Questi valori soglia devono essere identificati nell'immagine, usando lo strumento Identify
di QGIS per i pixel delle nuvole in banda 1 e banda 6.Il carattere
|
significa or
, cioè tutte le altre espressioni a seguire (es. "RT_LE70150532014090EDC00_B1" == 0
) identificato valori dei pixel uguali a 0 (che sono NoData), per ogni banda Landast, al fine di mascherare le linee nere dovute al SLC-off ed il bordo nero dell'immagine.
Potremmo usare lo stesso metodo di maschera delle nuvole anche per le immagini Landsat 8. Per l'immagine
LC80150532014050LGN00
carica le bande RT_LC80150532014050LGN00_B10
e RT_LC80150532014050LGN00_B2
, e usa la seguente Expression in Band calc:La condizione
"RT_LC80160532014057LGN00_B2" == 0
permette di mascherare il bordo nero dell'immagine.
Come si nota, ci sono ancora dei buchi (pixel non classificati) nella classificazione; avremmo bisogno di classificare altre immagini Landsat per riempire questi buchi.
Dopo la mascheratura di queste tre classificazioni, possiamo creare un mosaico, cioè la classificazione dell'intera area di studio.
Parte dei buchi non classificati sono stati riempiti con la classificazione della Landsat 7. Ovviamente sarebbero necessarie molte più immagini per riempire tutti i buchi.
Dopo la mascheratura di queste tre classificazioni, possiamo creare un mosaico, cioè la classificazione dell'intera area di studio.
Parte dei buchi non classificati sono stati riempiti con la classificazione della Landsat 7. Ovviamente sarebbero necessarie molte più immagini per riempire tutti i buchi.
Mosaico delle Classificazioni
Per creare un mosaico delle classificazioni, scriveremo un'espressione che riempirà i pixel non classificati dell'immagine Landsat 8 (ID
In QGIS apri le tre classificazioni. Copia l'espressione seguente in Band calc:
LC80150532014050LGN00
) con la classificazione dell'immagine Landsat 7 (ID LE70150532014090EDC00
). Inoltre, uniremo queste classificazioni con la terza (l'immagine Landsat 8 con ID LC80160532014057LGN00
).In QGIS apri le tre classificazioni. Copia l'espressione seguente in Band calc:
Disabilita il checkbox
Intersection
in Output raster e click Calculate
. Il risultato (es. classification_mosaic
) è mostrato nell'immagine seguente.
Nei prossimi passaggi andremo a verificare l'accuratezza e stimare le superfici di copertura del suolo.
Verifica dell'Accuratezza
La Verifica di Accuratezza è un passaggio importante di una classificazione di copertura del suolo. In questo tutorial andremo ad usare il
Training shapefile
come riferimento per valutare l'accuratezza della classificazione. Comunque, ci sono altri metodi che migliorano l'affidabilità della verifica (vedi http://fromgistors.blogspot.com/2014/09/accuracy-assessment-using-random-points.html ).In QGIS, carica il mosaico della classificazione e il
Training shapefile
usato per l'immagine LC80150532014050LGN00
(ovviamente potremmo usare anche gli altri shapefile). In SCP apri il tab Accuracy e click sul bottone Refresh list
. Seleziona classification_mosaic
come "classification to assess" e il Training shapefile
come "reference shapefile". Inoltre, seleziona MC_ID
come Shapefile field
. Click Calculate error matrix
e scegli la destinazione del raster di output (es. accuracy.tif
).Il processo produce una matrice di errore e un raster di errore che sono utili per valutare la qualità della classificazione.
Clip della Classificazione
Prima di calcolare l'area di ogni classe di copertura del suolo dobbiamo ritagliare la classificazione all'estensione dell'area di studio, che è il Costa Rica.
Scarica lo Shapefile del Costa Rica da http://data.fao.org/map?entryId=c7a0f990-88fd-11da-a88f-000d939bc5d8&tab=metadata (click su Download) fornito dalla FAO .
Estrai il file scaricato (
Scarica lo Shapefile del Costa Rica da http://data.fao.org/map?entryId=c7a0f990-88fd-11da-a88f-000d939bc5d8&tab=metadata (click su Download) fornito dalla FAO .
Estrai il file scaricato (
1173.zip
) e carica lo shapefile costa rica.shp
in QGIS (selezionando WGS84 come proiezione).
In questo caso, dobbiamo definire la proiezione di questo shapefile. In QGIS, apri il comando
Vector > Data management tool > Define current projection
; seleziona lo shapefile costa rica
come Input vector layer
e scegli EPSG:4326 - WGS 84
come riferimento spaziale, e click OK
.
Ora possiamo ritagliare il file
classification_mosaic.tif
. Carica la classificazione in QGIS. Apri il comando Raster > Extraction > Clipper
. Seleziona classification_mosaic
come input raster; imposta il file output (es. classification_clip.tif
), e imposta No data value
uguale a 0. In Clipping mode
abilita Mask layer
e seleziona costa rica
, poi click su OK
.
Abbiamo ritagliato la classificazione per l'estensione del Costa Rica (come si nota, avremmo bisogno di altre immagini Landsat per coprire tutta l'estensione del Costa Rica), e possiamo calcolare il report della classificazione.
Report della Classificazione
In SCP apri il tab Classification report e click sul bottone
Abilita
Si può salvare il report con un click su
Refresh list
.Abilita
Use NoData value
impostando il valore uguale a 0 e click sul bottone Calculate classification report
. Il report della classificazione mostra il conteggio dei pixel, l'area e la percentuale di ogni classe di copertura del suolo.Si può salvare il report con un click su
Save report to file
.
Abbiamo completato questo tutorial sulla classificazione della copertura del suolo di una vasta area, utilizzando più immagini Landsat e la creazione di un mosaico delle classificazioni. È importante sottolineare che i risultati della classificazione dipendono dalla stagione delle immagini. Pertanto, le immagini di input dovrebbero essere acquisite nello stesso periodo, al fine di evitare differenze dovute per esempio allo stato fenologico della vegetazione o al cambiamento della copertura del suolo avvenuto nel tempo.
Per informazioni o domande iscrivetevi al gruppo Facebook o alla Community Google+ del Semi-Automatic Classification Plugin.
Per informazioni o domande iscrivetevi al gruppo Facebook o alla Community Google+ del Semi-Automatic Classification Plugin.