This tutorial aims to analyze land cover change using the Semi-Automatic Classification Plugin (SCP) Postprocessing tools. Basically, we are going to assess land cover change from two raster classifications, and relate the changes to a land use vector file. An overview of several postprocessing tools is also provided.
The following is the video tutorial, and the following text illustrates the phases in detail.
The tools can be applied to any land cover classification, but we are going to use Copernicus data, which are freely available (as established by the EU Regulation No 1159/2013 of 12 July 2013) and cover the European countries. Of course, this tutorial is designed for demonstration purposes and it is not endorsed by the European Union. The original Copernicus data (produced with funding by the European Union) are downloaded from https://land.copernicus.eu/ and remain the sole property of the European Union.
Following, a brief description of the data we are going to use.
The Copernicus High Resolution Layers are raster classifications with 20m spatial resolution. Several land cover classes are available, but in this tutorial we are going to use the Imperviousness Density for 2012 and 2015. These data classify the degree of imperviousness (0-100% of impermeable cover of soil), which is the artificially sealed area. The Imperviousness Density was produced using automatic derivation based on calibrated Normalized Difference Vegetation Index. You can find the detailed product specifications here.
The Copernicus Corine Land Cover is a land use/land cover vector produced by standard methodology of photo-interpretation of satellite images. The vector is classified in 44 classes divided in 3 hierarchical levels with minimum mapping unit of 25 hectares. In this tutorial we are considering only the first level of Corine Land Cover 2012, divided in these classes:
- artificial surfaces
- agricultural areas
- forests and semi-natural areas
- wetlands
- water bodies
1. Refine the classifications with direct editing
You can download the data for this tutorial from this archive , or use your own data (two classification rasters and a land use vector).
For this tutorial, the original Copernicus data were modified by clipping the rasters to a small area over Florence (Italy).
Start QGIS and load the two rasters
IMD_2012.tif
and IMD_2015.tif
that are Copernicus Imperviousness Density for 2012 and 2015 respectively. As you can see, the rasters have values from 0 to 100, representing the degree of imperviousness.
It is useful to refine the classification by photo-interpretation, especially for data produced by semi-automatic processing.
We can use high resolution images or other services such as OpenStreetMap. For example you can follow this tutorial Download the Data to download satellite images, or you can download a subset of a Landsat 8 image, already converted to reflectance, from this link (about 27 MB, data available from the U.S. Geological Survey), unzip the downloaded file, and load the bands in QGIS.
First, we need to define a Band set containing a classification raster (this is required for drawing ROIs).
Open the tab Band set clicking the button in the SCP menu or the SCP dock. Click the button to refresh the layer list, and select the to add selected raster to the Band set 1.
IMD_2012
raster (just this raster is sufficient); then click
Optionally, we can create a band set for the satellite image to display a color composite; open the tab Band set and select all the Landsat bands in the list; click to add a new band set, then click to add selected rasters to the Band set 2.
In QGIS zoom to an area where we want to correct the classification. In this case we are going to manually remove a few pixels pretending they are classification errors.
We need to manually create a ROI, but first check that the Band set 1 is active. Now click the button in the Working toolbar. Left click on the map to define the ROI vertices and right click to define the last vertex closing the polygon. An orange semi-transparent polygon is displayed over the image, which is a temporary polygon (in this case we don’t need to define the Training input).
Now open the tool Edit raster opening the SCP menu and the submenu Use constant value enter 100 (we want to correct impervious pixels; in case we would like to correct not impervious pixels we would enter the value 0). The other options are fine. Therefore, click RUN to edit the raster.
Postprocessing
. Select the Input raster, for instance IMD_2012
. According to the legend of Imperviousness Density, in Attention: the input raster is directly edited; it is recommended to create a backup copy of the input raster before using this tool in order to prevent data loss.
Of course we could repeat these steps to edit any area of the raster.
TIP : Sometimes changes are not immediately visibile because the raster is not refreshed; try to zoom out and zoom in to refresh the view.
2. Classification report
It could be interesting to know the area of each land cover class. In order to get the area statistics, open the SCP menu and click the tab Classification report under the submenu
Postprocessing
.
Click the button to refresh the layer list, and select the ; next click RUN to start the calculation; the output report is saved in a text file and displayed in the tab Output.
IMD_2012
raster in Select the classification
We can repeat the same steps for the
IMD_2015
raster. Over the 86% of the area is not impervious.3. Reclassification
Before calculating land cover change it is convenient to reclassify the imperviousness degree into two classes: built-up and not built-up. A possible threshold for the distinction between built-up and not built-up is 30% (for further information read this document ). We can reclassify the raster using the SCP tool, obtaining the simple classification 1 = built-up and 0 = not built-up.
Open the tool Reclassification. In Select the classification select the raster twice to add two rows to the table. We need to enter the expressions illustrated in the following table.
IMD_2012
. Click the button Reclassification table
Old value | New value |
---|---|
raster < 30 | 0 |
raster >= 30 | 1 |
Uncheck the options Use code from Signature list and click RUN to start the reclassification. A new raster will be created (e.g.
BU_2012
).
Now select the to reclassify the 2015 raster (e.g.
IMD_2015
(the reclassification table is the same as before) and click RUN BU_2015
). Now the two reclassified rasters are loaded in the map and we can assing an appropriate symbology.4. Remove isolated pixels
We are going to compute the land cover change, but first we may want to remove isolated pixels in order to improve the analysis. In fact, single pixels may not represent real changes between the two classifications, because of geometrical shifts or isolated classification errors. Of course, this step is not always required, and it should be avoided if the purpose of the analysis is to find also the smallest changes.
We are going to use classification sieve for removing single pixels. Open the tool Classification sieve.
In Size threshold leave 2; all patches smaller the the selected number of pixels (i.e. single pixels) will be replaced by the value of the largest neighbour patch. Of course we could increase this value if we want to remove larger patches.
In Select the classification select the raster
BU_2012
. The option 4 in Pixel connection determines how pixels are considered connected, that is in a 3x3 window diagonal pixels are not considered connected. If we select the option 8 also diagonal pixels are considered connected.
Of course, we should repeat these steps also for raster
BU_2015
to create the new raster BU_2015_sieve
.5. Assess land cover change
Now we can use the tool to assess land cover change between the two classifications 2012 and 2015. Open the tool Land cover change.
This tool is quite straightforward. Click the button to refresh the layer list. In Select the reference classification select the
BU_2012_sieve
raster, that is the first classification. In Select the new classification select the BU_2015_sieve
raster, that is the latest classification.
Uncheck the option Report unchanged pixels, because we want to report only the pixels where the classification changed between 2012 and 2015. Now click RUN to create the new land cover change raster (e.g.
change
). Also, a text file is created (i.e. a file .csv separated by tab) containing the land cover change statistics.
The values of the land cover change raster represent a combination between reference and new classification, as described in the text file. In this case, only the value 1 is present that is the condition where
BU_2012_sieve
= 0 and BU_2015_sieve
= 1.
From the report we ca read that 520 pixels changed from 0 to 1, while no pixel changed from 1 to 0 between years 2012 and 2015.
6. Analyze the context of land cover changes
Now, it could be interesting to compare land cover change to other data such as land use, in order to analyze the context of new built-up areas. We are going to cross the land cover changes to the vector of Corine Land Cover; this way we can differentiate the new built-up areas according to Corine Land Cover classification system.
The original Corine Land Cover data were modified by clipping to a small area over Florence (Italy) and adding a field
Class_1
filled with the first level of classification.
Load in QGIS the Copernicus Corine Land Cover shapefile
CLC_2012.shp
previously downloaded. You can see the symbology of the first level Corine Land Cover classes that are:- artificial surfaces
- agricultural areas
- forests and semi-natural areas
- wetlands
- water bodies
Open the tool Cross classification. Click the button to refresh the layer list. In Select the classification select the Use NoData value and set the value 0, in order to exclude unchanged pixels (having value 0 in the
change
raster, that is our land cover change. Check change
raster) from the analysis.
In Select the reference vector or raster select the vector
CLC_2012
and in Vector field select the field Class_1
, containing the code of first level classes.
Now click RUN to create a new raster of comparison (e.g.
change_CLC
). The output will report the area of each combination between change
code and CLC_2012
code.
From the cross matrix we can evaluate the area in of built-up changes occurred in the 5 classes of Corine Land Cover classification.
Cross matrix
CLC_2012 | ||||||
---|---|---|---|---|---|---|
1 | 2 | 3 | 4 | 5 | ||
Change | 1 | 157600 | 48400 | 2000 | 0 | 0 |
The tool Cross classification can be very useful also for other analyses that involve the comparison with other data, such as population or flood risk, but this could be the subject of other tutorials.
7. Assess the spectral signature of changes
An optional step could be the assessment of the spectral signature of changes. We can download satellite images (see Tutorial 1: Your First Land Cover Classification for the details) and calculate spectral signatures for monitoring the changes through time.
We are going to use the Landsat 8 image downloaded at the beginning of this tutorial for calculating the spectral signature of changes. First, we need to create a Training input to store the spectral signatures calculated from the classes. In the SCP dock select the tab Training input and click the button to create the Training input (define a name such as
signatures.scp
). The path of the file is displayed and a vector is added to QGIS layers with the same name as the Training input (in order to prevent data loss, you should not edit this layer using QGIS functions).
Now open the tool Class signature opening the SCP menu and the submenu to start the calculation.
Postprocessing
. In Select the classification select the raster change_CLC
, thus we can distinguish the spectral signatures of changes. In Select input band set enter the number of the band set containing the Landsat 8 bands (i.e. 2). Now click RUN
After a while the spectral signatures are loaded in the Training input.
If the changes involved vegetation, we could calculate spectral signatures for images acquired in different seasons and assess the phenological variations of vegetation through spectral signatures. Also, these spectral signatures could be used as training input for further land cover classifications.
In order to display the signature plot, in the ROI Signature list highlight two or more spectral signatures (with click in the table), then click the button . The Spectral Signature Plot is displayed in a new window.
8. Export the changes to vector format
This is an optional step that may be useful for further analyses and integration with other data. We are going to convert the change raster to vector.
Open the tool Classification to vector. In Select the classification select the Use code from Signature list. Now click RUN to create a new vector (e.g.
change_CLC
raster and uncheck the change_vector
).
In the attribute table of this
change_vector
you can see the field C_ID
that represents the code of the change raster as described in Assess land cover change. Of course we could delete the polygons with code 0 (unchanged area), displaying only changes with code 1.For any comment or question, join the Facebook group about the Semi-Automatic Classification Plugin.