Unsupervised Classification and Polygonization With QGis

Video Walkthrough

If you would like a video walkthrough please follow this link to youtube!

https://youtu.be/H06wWFFNtgU

Introduction

There are lots of reasons for wanting to make a land cover classification. Whether that is to make beautiful maps, a temporal study of change or a more complicated statistical analysis. The goal of this tutorial is to be a jumping off point so to speak into the profound and complicated world of classifications.

In this tutorial we will be looking at how to perform an unsupervised classification; that is, a classification where the computer program chooses the classes without a training data set

Programs and Plugins Required

Qgis (3.22 was used in this tutorial) https://qgis.org/en/site/forusers/download.html
7Zip https://www.7-zip.org/
SCP Plugin (download through QGis)

Acquiring Data

For this tutorial pretty much any remote sensing data can be used. Depending on what data it it it may require pre-processing in order to produce a usable product, and that is outside the scope of this tutorial.

A good resource for free remote sensing data us the USGS Earth explorer. It involves making an account but it is such a good resource for remote sensing data that it is completely worth the extra step.

Sign up here https://ers.cr.usgs.gov/login

Now that you have made an account or signed into your existing account it is time to find the data to use.

You can use pretty much whatever data you see fit but Landsat8 is a very common choice and so that is what's being used in this tutorial. That being said the first step is to select which data set you want to use. To get to Landsat 8 image used in this tutorial select Landsat ->Landsat Collection 2 Level-1 ->Landsat 8 OLI/TIRS C2 L1

The next step is to pick the additional criteria. In this tutorial the image used was LC08_L1TP_046007_20200722_20200911_02_T1 this code can be copy and pasted into the "Landsat Product Identifier L1" text box.

Next step is to click results.

The results should bring up one image. Click on the footprint icon to show where it is on the world map (should be in northern Canada).

To download the image click on the download icon, the one with the little green arrow four icons left of the footprint.

Select the product options and right at the top click on the download button for the entire Landsat 8 Collection 2 Level-1 Product Bundle. It should be a big file (1.13 Gigabytes)

Unzipping Data

The data downloaded from USGS Earth Explorer is in .tar format and most windows machines are unable to extract that type of file without an additional program. You might have to download 7Zip (see Programs and Plugins Required)

To unzip the data simply right click on the file you wish to extract, select 7Zip and select extract here. This will not create a new folder it will simply extract the file into the folder that it is already in.
Once the extraction is done feel free to delete the .tar file as it is no longer needed.

Importing and Displaying Data

To import data, launch QGis and create a new project. once the project it loaded go to the rater tab at the top:
Raster-> miscellaneous -> Build Virtual Raster. This will open up the build virtual raster tool. Click the thee dots beside the input layers tab.
The input layers are going to be everything that was extracted from the .tar file. To do this use the:
Add Directory button
navigate to the folder where the .tar file was extracted and click:
Select Folder button
This will bring up everything from that file but not everything is needed. De-select every layer that doesn't end in B followed by a number.
The only layers that should be selected end with:
B1.TIF
B2.TIF
B3.TIF
B4.TIF
B5.TIF
B6.TIF
B7.TIF
B8.TIF
B9.TIF
B10.TIF
B11.TIF
Each file represents a different wavelength that the sensor collected. Landsat 8 uses a multispectral sensor so it collects the entire light spectrum across 11 wave lengths from blue to infrared. For instance B2 is blue, B3 is green, B4 is red NIR is B5 and thermal infrared is B11. Check out the USGS website for more information:https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites?qt-news_science_products=0#qt-news_science_products
With all the correct bands selected hit the ok button.
Change the resolution to highest
Check the box for place each band in separate band file.
Click run
The image should be displayed in the QGis project. The colouring might look weird and that is because it is currently a false colour composite. To change it to a natural colour composite
Right click on the Virtual layer-> select properties-> select symbology
Select Band 04 for red
Select Band 03 for green
Select band 02 for blue
Click apply and then ok to close the window.
Now it c=should be a natural colour composite.

Trim to Mask Area

It is possible to do the classification for the entire Landsat image but it could take a long time so for this tutorial we are going to trim the raster to a mask layer. This will allow us to classify a smaller subset of the region.

To do this a new layer needs to be made:
Layer tab-> Create layer-> New shape file layer
The geometry type can be a polygon, but in this tutorial we are going to simulate a study site and create a buffer around it.
Geometry type should be point. Click ok.
To create the point click:
On the new layer-> toggle editing (pencil icon)-> Add point feature (button two down from toggle editing)-> click where you want the mask area to be centred -> set the point id to 1 -> click ok ->Click save edits (between toggle editing and add point feature) -> click toggle editing again.

Make the buffer by going to:
Vector-> Geoprocessing tools -> buffer
The input layer will be the point layer you just created and select the distance you want to classify, this will be the radius from the point. I used 10km, if km is not an option and it is in degrees this can be changed by changing the projection. Click run. this will draw a circle on your map.

To clip the raster to the new buffer layer go to:
Raster-> Extraction-> Clip raster by mask layer
Input layer will be the virtual raster, mask layer will be the buffer. Click run.
There we go the raster has been clipped.

Unsupervised Classification and Filtering

Time to classify the image. You need the SCP tool. To download it go to:
Plugins-> manage nd install plugins-> search up scp-> select the semi-automatic classification plugin-> install-> close.
to use the SCP tool: SCP->Band set. This will open the tool. Under the Multiband image lost option right at the top go to the drop down menu and select the Clipped (mask) band created. Hit the refresh button on the right (second from the top), select all the bands and hit the Add bands to bandset button. Those steps will bring all the bands into the tool.
Next step is to go to the clustering tab in to tool.
The inputs will be:
method: K-means
Number of classes: 5
Use random seed signatures
All other settings can be left default for now. Hit run, it will prompt you to save it so navigate to the correct folder and hit save. The classification might take a few minutes be patient.

There might be some noise speckling in the classification and classes where there shouldn't be, this is likely due to noise in the imagery. To remove the noise open up the SCP tool again:
post processing tab-> Classification sieve. Set the threshold to 25 (makes the smallest area classified 25 pixels) and set the pixel connection to 8 (this means pixels touching on the diagonal are considered connected and in one cluster). Hit Run. You will be prompted to save the output, name it and hit save. The effect this will have is a more smoothed out classification.

Polygonization (Raster to Vector Conversion)

Rasters can sometimes be difficult to use when calculating land cover statistics about a region so a good option is to transform the raster layer into a vector layer. This can be done with the built in polygonization tool. This transforms every cluster in the classification into a separate polygon. To do this go to:
Raster-> Conversions-> Polygonize. The input layer will be the filtered layer created in the pervious step. Hit run.
The new polygon layer may end up one colour. To represent it as the different land use types:
right click on the layer-> Symbology-> Categorized->Set the value to DN (fid is just the id of each individual polygon DN is the value taken from the raster)-> Classify-> ok.

Field Calculator

To get some more useful information from this new poligonized classification right click on the layer and open up the attribute table. The DN field is corresponds to the number of the land cover from the raster version of the classification but that is difficult to understand at a glance so we are going to add a new field that displays the name of the land cover type.
Open field calculator (from the attribute table screen)->Create new field-> name it land cover-> Output field type, text (string). In the expression box type in: CASE WHEN "DN" IS 1 THEN 'Land' WHEN "DN" IS 2 THEN 'Bound' WHEN "DN" IS 4 THEN 'Water' ELSE 'Ice' END
This expression will correspond all the DN values to a name that the field calculator will put in a new field. Change the numbers as necessary for you classification as it may not end up being identical to mine. Click ok.
this will add a new field with a land cover name for every polygon.

The SCP tool classifies a box round the entire region and now is a good time to get rid of that:
Click on the title for the field we just made, this will sort by name-> click on toggle editing in the attribute table-> select all the rows you want to delete (in my case everything named bound)-> hit the delete features button.

Resources

Unsupervised Classification using SCP: https://fromgistors.blogspot.com/2021/01/unsupervised-classification-using-scp7.html
USGS Earth Explorer: https://earthexplorer.usgs.gov/
Loading Landsat 8 into Qgis: https://www.youtube.com/watch?v=mBk2VIMawRE
Raster Conversion: https://docs.qgis.org/3.22/en/docs/user_manual/processing_algs/gdal/rasterconversion.html?highlight=raster%20vector
USGS band designations: https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites?qt-news_science_products=0#qt-news_science_products