In this article, Luca Congedo presents a tutorial for the land cover classification of cropland. In particular, he uses Landsat imagery acquired over the State of Kansas, near the city of Ulysses, using the new version 2.3.2 of the Semi-Automatic Classification Plugin for QGIS.
Before the classification process, we are going to convert the DN values to reflectance. The reflectance is defined as the ratio between reflected and incident energy on a surface (see here for more information about the conversion to reflectance of Landsat images).
First, download the image from here. The image is a subset of the entire scene, including the following Landsat bands (each band is a single 16 bit raster):
2 - Blue;
- 3 - Green;
- 4 - Red;
- 5 - Near-Infrared;
- 6 - Short Wavelength Infrared 1;
- 7 - Short Wavelength Infrared 2.
1. Conversion of raster bands from DN to reflectance
|Original DN spectral signature|
|TOA spectral signature|
DOS1 corrected spectral signature
- Open QGIS and start the Semi-Automatic Classification Plugin; in the main interface select the tab Pre processing > Landsat;
- Select the directory that contains the Landsat bands (and also the required metafile MTL.txt), and select the output directory where converted bands are saved;
- Check the option Apply DOS1 atmospheric correction, and click Perform conversion to convert Landsat bands to reflectance;
- At the end of the process, converted bands are loaded in QGIS.
2. Definition of the classification inputs
In order to create the color composite in QGIS, we create a virtual raster. Then, we define the band set that we are going to classify. Finally, we create a training shapefile which stores the ROIs (i.e. Regions Of Interest) that define the land cover classes.
- From the menu Raster select Miscellaneous > Build Virtual Raster (catalog); click the button Select... and select the bands 3, 4, and 5; select the output file (for instance rgb.vrt); check Separate (bands will be separated) and click OK;
- In the plugin dock ROI creation click the button Band set beside Select an image; the tab Band set will open; click the button Select All, then Add rasters to set (order the band names in ascending order, from top to bottom using ↑ and ↓ arrows);
- In order to create the training shapefile (a polygon layer having several attribute fields that will store ROIs), in the dock ROI creation click the button New shapefile, and select where to save the shapefile (for instance ROI.shp).
3. Creation of the ROIs
We need to create several ROI, considering the spectral variability of land cover classes. It is important that ROIs represent homogeneous areas of the image, therefore we are going todraw the ROIs using a region growing process (i.e. a segmentation of the image, grouping similar pixels).
The following are the land cover classes that we are going to identify in the image:
- crop (e.g. fields with green vegetation)
- low vegetation (e.g. fields without green vegetation, or shrubland)
- built-up (e.g. artificial areas, buildings and asphalt)
- farms (e.g. farm areas)
- bare soil (e.g. soil without vegetation)
- water (e.g. surface water)
- In order to create a ROI, in the dock ROI creation click the button + beside Create a ROI and then click any pixel of the image; zoom in the map and click on a red pixel of a circular field; after a few seconds the ROI polygon will appear over the image (a semitransparent orange polygon);
- Under ROI definition type a brief description of the ROI inside the field Information (e.g. crop; notice that descriptions are not used during the classification process, but are useful for distinguishing ROIs); the field ID (i.e. identifier of the class) is used as reference for the land cover classification, therefore it is important that each category has a unique ID value (ROIs sharing the same ID are evaluated as a single ROI); since we need more than one ROI per class, set the ID = 11, so that the first figure represents the class, and second one the ROI number (other crop ROIs will have values from 11 to 19; in the same way we set ID from 21 to 29 for low vegetation, and so on);
- In order to save the ROI to the training shapefile click the button Save ROI to shapefile;
- In the plugin main interface, select the tab ROI tools > Spectral signature, which displays the plots of selected ROIs, and select the item crop; if the checkox Plot σ is checked, then the plot will display the standard deviation of each ROI.
Repeat the above steps for every class (in order to get the desired ROI, change the values for Minimum ROI size and Range radius); after the collection of several ROIs, it is useful to visualize the ROI spectral signatures in order to avoid the collection of ROIs that are spectrally similar. Following, some examples of ROIs for the land cover classes.
You can download the final ROI shapefile (in which I saved 14 ROIs) from here.
4. Classification of the study areaIt is useful to perform classification previews, in order to assess the quality of classification (for example a common problem is soil classified as built-up and vice versa). If the result of the preview is not good, then we need to go back to the phase of ROI creation, delete the faulty ROIs and/or create new ones.
The classification can be performed through several algorithms; this time we are going to use the Spectral Angle Mapping classification.
- In the Classification dock, under Classification preview set Size = 500 (i.e. the side of the classification preview in pixel unit); click the button + and then click on the image; after a few seconds, the classification preview will be loaded (with default colors);
For a better interpretation of classes, it is useful to assign a classification color to each class.
Also, we can save a visualization style in QGIS, and set this .qml file as the default style for every classification or preview.
- In the panel Layers, left click on the classification preview and select Properties > Style; select the render type Singleband pseudocolor and change the colors and labels of classes, according to the training ROIs; then, click Save style ... to save the .qml file (e.g. classification.qml);
- Under Classification style click the button Select qml to select the file classification.qml; thus the next classification or preview will be loaded with this style;
- click the button Perform classification and select where to save the output (e.g. classification.tif).
Here you can download the final land cover classification.
5. Calculation of classification accuracyIt is useful (and necessary) to assess the accuracy of land cover classification, in order to understand the reliability thereof, and to identify map errors.
Please, notice that this classification is only for demonstration purposes, because an accurate land cover classification would require ancillary data and field survey.
However, in order to assess the classification accuracy, we can compare the land cover classification to the ROIs we have created (theoretically, image pixels belonging to a certain ROI should be classified as that ROI's class).
- Select the tab Post processing > Accuracy of the plugin main interface; select the classification.tif beside Select a classification to assess and select the ROI shapefile beside Select the reference shapefile; then click the button Calculate error matrix and the matrix will be displayed; you can save the error matrix by clicking the button Save error matrix to file;
Generally, classification accuracy is considered good if it is > 80%. For this classification, we can see that the overall accuracy of the classification is about 91.4% (the error matrix is reported below). Also, it is useful to consider the error for single classes (i.e. omission and commission).
In order to estimate the cropland in the area, we are going to reclassify the classification raster in fewer classes, according to the main land cover classes.
- From the Processing toolbox navigate to SAGA > Grid - Tools > Reclassify Grid Values; in the tool window, under Grid selectclassification.tif; under method select  simple table;
- Under Lookup table click the button ... to open the window Fixed Table; the lookup table has 3 columns: minimum, maximum, and new. Minimum and maximum define the reclassification range for the new value, according to the following table.
- Click "OK" and after a few seconds the reclassified classification will be loaded.
- Select the tab Post processing > Classification report of the plugin main interface; select Reclassified grid beside Select a classificationand click Calculate classification report; after a few seconds the report will be displayed (see the following table).
|Class||PixelSum||Percentage %||Area [metre^2]|
From these results, we can see that in the study area (at the moment of image acquisition) there are about 738 square kilometers of crop class (class 1, about 10% of the study area) and 6231 square kilometers of low vegetation class (class 2, about 87%). Also, we can notice that about 14 square kilometers are occupied by farms (class 4) and about 18% is built-up (class 3).
Of course, these numbers are the result of a tutorial for demonstration purposes; for a better estimation of class areas (and possibly the identification of the crop type), we need field data which can improve the creation of ROIs and help for the assessment of classification accuracy.
However, the aim of this tutorial is to show how land cover classifications can be performed simply and rapidly using open source programs, which are useful for several activities, such as land cover monitoring or precision farming.
The following is the video of this tutorial.
Landsat image conversion to reflectance and DOS1 atmospheric correction
Lλ = ML * Qcal + AL
- ML = Band-specific multiplicative rescaling factor from Landsat metadata (RADIANCE_MULT_BAND_x, where x is the band number)
- AL = Band-specific additive rescaling factor from Landsat metadata (RADIANCE_ADD_BAND_x, where x is the band number)
- Qcal = Quantized and calibrated standard product pixel values (DN)
“For relatively clear Landsat scenes, a reduction in between-scene variability can be achieved through a normalization for solar irradiance by converting spectral radiance, as calculated above, to planetary reflectance or albedo. This combined surface and atmospheric reflectance of the Earth is computed with the following formula” (NASA, 2011, p. 119):
ρp = (π * Lλ * d^2 )/ (ESUNλ *cosθs)
- ρp = Unitless planetary reflectance, which is “the ratio of reflected versus total power energy” (NASA, 2011, p. 47)
- Lλ = Spectral radiance at the sensor's aperture (at-satellite radiance)
- d = Earth-Sun distance in astronomical units (provided with Landsat 8 metafile, and an excel file is available from
- ESUNλ = Mean solar exo-atmospheric irradiances
- θs = Solar zenith angle in degrees, which is equal to θs = 90° - θe where θe is the Sun elevation
It is worth pointing out that Landsat 8 images are provided with band-specific rescaling factors that allow for the direct conversion from DN to TOA reflectance. However, the effects of the atmosphere (i.e. a disturbance on the reflectance that varies with the wavelenght) should be considered in order to measure the refletance at the ground. As described by Moran et al. (1992), the land surface reflectance (ρ) is:
ρ = [π * (Lλ - Lp) * d^2]/ [Tv * ( (ESUNλ * cosθs * Tz ) + Edown )]
- Lp is the path radiance
- Tv is the atmospheric transmittance in the viewing direction
- Tz is the atmospheric transmittance in the illumination direction
- Edown is the downwelling diffuse irradiance
Therefore, we need several atmospheric measurements in order to calculate ρ (physically-based corrections). Alternatively, it is possible to use image-based techniques for the calculation of these parameters, without in-situ measurements during image acquisition.
The Dark Object Subtraction (DOS) is a family of image-based atmospheric corrections.
Chavez (1996) explains that “the basic assumption is that within the image some pixels are in complete shadow and their radiances received at the satellite are due to atmospheric scattering (path radiance). This assumption is combined with the fact that very few targets on the Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent”. It is worth pointing out that the accuracy of image-based techniques is generally lower than physically-based corrections, but they are very useful when no atmospheric measurements are available as they can improve the estimation of land surface reflectance.
The path radiance is given by (Sobrino, et al., 2004):
Lp = Lmin – LDO1%
- Lmin = “radiance that corresponds to a digital count value for which the sum of all the pixels with digital counts lower or equal to this value is equal to the 0,01% of all the pixels from the image considered” (Sobrino, et al., 2004, p. 437), therefore the radiance obtained with that digital count value (DNmin)
- LDO1% = radiance of Dark Object, assumed to have a reflectance value of 0,01
Therfore for Landsat images:
Lmin = ML * DNmin + AL
LDO1% = 0,01* [(ESUNλ * cosθs * Tz ) + Edown] * Tv / (π * d^2)
Lp = ML* DNmin + AL – 0,01* [(ESUNλ * cosθs * Tz ) + Edown] * Tv / (π * d^2)
There are several DOS techniques (e.g. DOS1, DOS2, DOS3, DOS4), based on different assumption about Tv, Tz , and Edown .
The simplest technique is the DOS1, where the following assumptions are made (Moran et al., 1992):
- Tv = 1
- Tz = 1
- Edown = 0
Therefore the path radiance is:
Lp = ML* DNmin + AL – 0,01* ESUNλ * cosθs / (π * d^2)
ρ = [π * (Lλ - Lp) * d^2]/ (ESUNλ * cosθs)
ESUN [W/(m2 * μm)] values for Landsat sensors are provided in the following table.
BandLandsat 4*Landsat 5**Landsat 7**11957198319972182517691812315571536153341033103110395214.9220230.8780.7283.4484.90
* from Chander & Markham (2003)
** from Finn, et al. (2012)
For Landsat 8, ESUN can be calculated as (from http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html):
ESUN = (π *d2) * RADIANCE_MAXIMUM / REFLECTANCE_MAXIMUM
where RADIANCE_MAXIMUM and REFLECTANCE_MAXIMUM are provided by image metadata.
-Chander, G. & Markham, B. 2003. Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges Geoscience and Remote Sensing, IEEE Transactions on, 41, 2674 – 2677
-Chavez, P. S. 1996. Image-Based Atmospheric Corrections - Revisited and Improved Photogrammetric Engineering and Remote Sensing, [Falls Church, Va.] American Society of Photogrammetry, 62, 1025-1036
-Finn, M.P., Reed, M.D, and Yamamoto, K.H. 2012. A Straight Forward Guide for Processing Radiance and Reflectance for EO-1 ALI, Landsat 5 TM, Landsat 7 ETM+, and ASTER. Unpublished Report from USGS/Center of Excellence for Geospatial Information Science, 8 p
-Moran, M.; Jackson, R.; Slater, P. & Teillet, P. 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output Remote Sensing of Environment, 41, 169-184
-NASA (Ed.) 2011. Landsat 7 Science Data Users Handbook Landsat Project Science Office at NASA's Goddard Space Flight Center in Greenbelt, 186
-Sobrino, J.; Jiménez-Muñoz, J. C. & Paolini, L. 2004. Land surface temperature retrieval from LANDSAT TM 5 Remote Sensing of Environment, Elsevier, 90, 434-440
Reprinted under Creative Commons License 3.0. Originally published from the "From GIS to Remote Sensing" blog.