Calculating Raster Area (QGIS3)

Many applications require quantifying the landuse patterns in a region. Land Use Land Cover (LULC) datasets come as raster files where each pixel is assigned a class value. GIS analysts often needs to generate reports based on this data by computing the area per class in a given region. QGIS comes with many built-in tools to calculate and summarize raster area.

Σημείωση

Historically the suggested approach for calculating areas for rasters was to convert the raster to a vector layer and use vector area calculation techniques. This approach is quite computationally intensive and error-prone. The recommended approach is to use the processing tool Raster layer unique values report which can directly compute the pixel areas. If you have a layer with many polygons and need areas for each of them, you can use the Zonal histogram tool to get pixel counts for each class and then multiply it with the area of each pixel.

Overview of the task

We will be using a raster layer with 11 land cover classes and calculate the area of each class within a national park. We will also post-process the results to create a spreadsheet with class names and areas.

Other skills you will learn

  • How to apply symbology to a layer from a style file in the .qml format.

  • How to write expressions with multiple if/else conditions using the CASE statement.

  • How to export a QGIS table as an Excel spreadsheet.

Get the data

We will be downloading the following datasets

  1. World Database on Protected Areas (WDPA): We will download the shapefile for the boundary of the Kaziranga National Park in India.

  2. ESA WorldCover 2020: The European Space Agency (ESA) WorldCover 10 m 2020 product provides a global land cover map for 2020 at 10 m resolution. We will download the tile covering our region of interest.

Park Boundary

  1. Go to the Protected Planet website, and click on the search toolbox. Search for Kaziranga National Park.

../../_images/data01.png

  1. The protected vector boundary will be shown as a search result. Click on it to view the page for the Kaziranga National Park.

../../_images/data02.png

  1. This page will contain additional information like total area, created year, etc. Click on the Download and click the SHP to download the data in Shapefile format.

../../_images/data03.png

  1. Two options for download will be prompted. Click continue under Non Commercial Use. Now a zip file containing the national park boundary will be downloaded.

../../_images/data04.png

Landcover Data

  1. Go to the ESA WorldCover website and click on the DATA ACCESS menu.

../../_images/data112.png

  1. Scroll down to the DATA DOWNLOAD section and click on the link to open the WorldCover viewer

../../_images/data121.png

  1. You need to create a free account to download the data. Click on the Register link on the top-right corner. Follow the instruction given to create a new account.

../../_images/data13.png

  1. After creating the account, log in using the credentials. Our area of interest for this tutorial is the Kaziranga National Park. Zoom to the North-East India region. Once you are zoomed in enough, the landcover image tiles bounding box will start to appear.

../../_images/data14.png

  1. Search and locate the N24E093 tile region.

../../_images/data15.png

  1. Select the tile and click NEXT.

../../_images/data16.png

  1. Click on the DOWNLOAD to download a zip file containing the landcover information in raster format. Make sure to select WorldCover Version 1 data for download.

../../_images/data17.png

  1. We will also download a symbology file provided by ESA. Visit the ESA WorldCover Data Access page., Scroll-down to the Symbology section. Click on the QGIS to download the ESAWorldCover_ColorLegend.qml file which can be used to style the raster layer with approproate colors and class labels.

../../_images/data18.png

For convenience, you may directly download a copy of the datasets from the links below:

Data Source [WDPA] [WorldCover]

Procedure

  1. Unzip all the downloaded files. In the Browser, locate the folder containing the boundary file WDPA_WDOECM_Apr2022_Publicc_10744_shp-polygons.shp and drag-and-drop it into the QGIS canvas.

../../_images/013.png

  1. Now locate the worldcover raster tile ESA_WorldCover_10m_2020_v100_N24_E093_Map.tif and drop it into the QGIS canvas.

../../_images/023.png

  1. You will now have both the vector boundary and landcover raster layer loaded in the Layers panel. Let’s clip the landcover raster to the national park boundary. Go to Processing ‣ Toolbox to open Processing toolbox. Search for and locate the GDAL ‣ Raster extraction ‣ Clip raster by mask layer algorithm. Double-click to launch it.

../../_images/033.png

  1. In the Clip Raster by Mask Layer dialog, choose the ESA_WorldCover_10m_2020_v100_N24_E093_Map layer as the Input layer and WDPA_WDOECM_Apr2022_Publicc_10744_shp-polygons layer as Mask Layer. Enter -9999 in Assign a specified nodata value to output bands section.

../../_images/043.png

  1. Now open the Advanced Parameters section and choose High Compression in Profile. Now under Clipped (mask), click on the ... and select Save To File…. Enter the file name as worldcover_clipped.tif. Click Run.

../../_images/053.png

  1. Now the worldcover_clipped layer will be loaded in the QGIS canvas. Right-click the ESA_WorldCover_10m_2020_v100_N24_E093_Map layer and select Remove Layer…

../../_images/063.png

  1. Both our layers come in the Geographic CRS EPSG:4326. This CRS has units of degrees and is not suitable for calculating area. We must first reproject the layers to a Projected CRS. for regional analysis like these, UTM is a good choice for a projected CRS. We will reproject the layers to the CRS for the local UTM zone. Open the Processing toolbox and search for Vector general ‣ Reproject layer algorithm. Double-click to launch it.

../../_images/073.png

  1. In the Reproject Layer dialog, choose the WDPA_WDOECM_Apr2022_Publicc_10744_shp-polygons layer as the Input layer, click on the Select CRS button under Target CRS.

../../_images/083.png

  1. Our area of interest falls in the UTM Zone 46N. Search for 46 N and select the WGS 84 / UTM zone 46N CRS.

../../_images/093.png

Σημείωση

To find out which UTM zone for your region, refer to the What UTM Zone am I in website.

  1. In the Reprojected section, click ... and select Save To File…. Enter the name as nationalpark_reprojected.gpkg. Click Run.

../../_images/1010.png

  1. Now the nationalpark_reprojected layer will be loaded in canvas. Right-click the WDPA_WDOECM_Apr2022_Publicc_10744_shp-polygons layer and select Remove Layer… to remove it. Now we will reproject the raster layer. In the Processing Toolbox, search and locate GDAL ‣ Raster projections ‣ Warp (reproject)

../../_images/1113.png

  1. In the Warp (Reproject) dialog choose worldcover_clipped as the Input layer, select WGS 84 / UTM zone 46N CRS in Target CRS. Open the Advanced Parameters and choose High Compression in Profile.

../../_images/1211.png

  1. Now under Reprojected, click on ... and select Save To File…. Enter the name as worldcover_reprojected.tif. Click Run.

../../_images/139.png

  1. Now the worldcover_reprojected layer will be loaded in the canvas, remove the worldcover_clipped layer. Let’s set the project layer to the UTM zone. Click on any layer and choose Layer CRS ‣ Set Project CRS from Layer.

../../_images/148.png

  1. Now the project CRS will be updated. Let’s set the symbology of the raster layer as per the class names and colors of the ESA WorldCover dataset. Right-click on the worldcover_reprojected layer and click Properties…

../../_images/158.png

  1. In the Layer Properties dialog, choose Symbology. You can notice the Layer colors are visualized in a white-black tone. To fix this, click on the Style ‣ Load Style…. Browse and select the ESAWorldCover_ColorLegend.qml file.

../../_images/167.png

  1. Now you can see the updated symbol colors and class descriptions. Click OK.

../../_images/178.png

  1. Expand the worldcover_reprojected layer in the Layers panel to see the legend with correct class descriptions.

../../_images/186.png

  1. Now let’s calculate the area for each class. In the processing toolbox, search and locate the Raster layer unique values report tool. Double-click to open it.

../../_images/196.png

  1. In the Raster Layer Unique Values Report dialog, choose the Input layer as worldcover_reprojected. Under the Unique values table click on ... and choose Save to File…. Enter the name as class_areas.gpkg. Click Run.

../../_images/206.png

  1. Now the class_areas layer will be added to the Layers panel. Right-click on the layer and click Open Attribute Table. The column m2 contains the area for each class in square meters.

../../_images/2110.png

  1. Let’s convert the area to square kilometers. In the Processing Toolbox, search and select Vector table ‣ Field Calculator.

../../_images/226.png

  1. In the Field Calculator dialog, select the class_areas layer in the Input Layer. Enter the Field name as area_sqkm. In the Result field type choose Float. In the Expression window, enter the below expression. This will convert the sqmt to sqkm and round the result to 2 decimal places. Under the Calculated click on ... and choose Save To File… . Enter the name as class_area_sqkm.gpkg. Click Run.

round("m2"/ 1e6, 2)
../../_images/235.png

  1. Now the class_area_sqkm layer will be loaded in canvas. Open the Attribute table and examine the newly added area_sqkm column. You will notice that the Value column contains numbers for each class. To make the results easier to interpret, Let’s also add the description for each class number. The class descriptions are available in the ESA Product User Manual.

../../_images/246.png

  1. Open Field Calculator, and select the class_areas_sqkm layer in Input Layer. Enter the Field name as landcover, in the Result field type, choose String. In the Expression window enter the below expression. This expression uses the CASE statement to assign a value based on multiple conditions. Under the Calculated click on ... and choose Save To File… . Enter the name as class_area_with_landcover.gpkg. Click Run.

CASE
WHEN "value" = 10 THEN 'Tree cover'
WHEN "value" = 20 THEN 'Shrubland'
WHEN "value" = 30 THEN 'Grassland'
WHEN "value" = 40 THEN 'Cropland'
WHEN "value" = 50 THEN 'Built-up'
WHEN "value" = 60 THEN 'Bare / sparse vegetation'
WHEN "value" = 70 THEN 'Snow and Ice'
WHEN "value" = 80 THEN 'Permanent water bodies'
WHEN "value" = 90 THEN 'Herbaceous wetland'
WHEN "value" = 95 THEN 'Moss and lichen'
WHEN "value" = 100 THEN 'Mangroves'
END
../../_images/256.png

  1. Now the class_area_with_landcover layer will be loaded in canvas. Open the Attribute table. The landcover column will contain the landcover name against each landcover value.

../../_images/265.png

  1. Let’s export this result as an excel file. Before export we will also organize the table and remove unwanted fields. In the Processing Toolbox, search and select Vector table ‣ Refactor fields.

../../_images/276.png

  1. In the Refactor Fields dialog, select the class_area_with_landcover layer in the Input Layer. Select all columns except area_sqkm and landcover, then click Delete selected field.

../../_images/286.png

  1. You can also change the order of fields in the table using the Move Selected Fields button. Once you are done with the edits, click on the ... button next to Refactored and choose Save To File…. Select XLSX Files (*.xlsx) as the format. Enter the file name as park_area_by_landcover.xlsx and click Save. In the Refactor Fields dialog, click Run to apply your changes.

../../_images/295.png

  1. The result will be a spreadheet with landcover and area_sqkm columns.

../../_images/305.png

If you want to report any issues with this tutorial, please comment below. (requires GitHub account)