Capacity Building for GIS-based SDG Indicator Analysis with Global High-resolution Land Cover Datasets
ISPRS Educational and Capacity Building Initiative 2022
In the framework of the United Nations Sustainable Development Goals (UN SDGs), the support of geospatial data and technologies has turned out to be critical for both the assessment and the monitoring of key indicators, revealing the trajectory of our planet and society towards sustainability. The increasing availability of global open geospatial datasets - above all the global high-resolution land cover (HRLC) datasets - opens newsworthy opportunities for the computation and comparison of these indicators across different geographical regions as well as multiple spatial and temporal scales. The added value of these datasets is tangible, especially for developing countries, where often such information is only partially available from local authorities. Nevertheless, there are still several barriers to their proficient use due to the lack of data management and processing capacity using proper Geographic Information Systems (GIS) software tools.
In view of the above, this project, supported by the Educational and Capacity Building Initiative 2022 of the International Society of Photogrammetry and Remote Sensing (ISPRS), addresses the creation of open training material covering the complete learning process of discovering, accessing and manipulating global open geospatial datasets for computing SDG indicators, with a focus on those directly connected to marine and terrestrial ecosystems, urban environment, and climate. To ensure the widest possible accessibility, the material primarily leverages the Free and Open Source Software (FOSS) QGIS and it is released under a Creative Commons Attribution 4.0 License (CC BY 4.0).
Background
Introduction to United Nations Sustainable Development Goals
The 2030 Agenda for Sustainable Development, adopted by all United Nations Member States in 2015, provides a “shared blueprint for peace and prosperity for people and the planet, now and into the future”. At its core are the 17 Sustainable Development Goals (SDGs), which are an urgent call for action by both developed and developing countries in a global partnership (Fig. 1.1.1). The 17 SDGs comprise 169 targets which can be tracked and monitored through 231 unique indicators. These indicators aim to balance the economic, social and ecological dimensions of sustainable development, by means of measurable and comparable indexes [1].
Sustainable Development Goals. Source: https://www.un.org/en/sustainable-development-goals
Footnotes
Geospatial data and technologies for SDGs
Despite the development of the Global SDG framework was originally based on traditional statistical data [1], many statistical offices and governmental agencies have recognised the need for geospatial data to augment non-spatial information and provide new and consistent data sources, which are critical to SDGs monitoring. In fact, It has been estimated that approximately 20% of the SDG indicators can be interpreted and measured either through the direct use of geospatial data itself or through integration with complementing statistical data [2]. GIS, Earth observations and global-coverage geospatial data provide an important link to enable consistent comparison among countries, provide granularity and disaggregation of the indicators and communicate their geographic dimensions [3].
Key global geospatial data products that cover major thematic areas of the biosphere and society (e.g. land cover, vegetation productivity, forests, wetlands, surface water, human settlements, etc.) are strongly supporting the methodological development and measurement of a number of SDG indicators. Generally, national or local geospatial data provide (where available) a higher spatial and temporal resolution than coarser global products and - in turn - better capabilities for accurate monitoring of SDG indicators. However, the uneven availability of such data often prevents consistent assessment of SDGs’ progress across countries worldwide. In this context, the use of open global geospatial datasets represents a quick, efficient and cost-effective solution to provide a first comprehensive assessment before national data is produced and analyzed [4].
Footnotes
Data
Global High-resolution Land Cover datasets
A notable class of key global geospatial products, which plays a major role in many SDGs, is the class of high-resolution global land cover (HRLC) maps [1]. The land cover stands for all the geospatial information describing the physical and biological cover of the Earth’s surface. Land cover data products at high spatial resolution have a critical role in many scientific and policy-making applications, including climate modeling, natural resources management, landscape and biodiversity preservation, urbanization monitoring, and spatial demography. In the last two decades, the advancement of satellite land Remote Sensing as well as the availability of global coverage satellite imagery with high temporal resolution - available as free and openly-licensed data - have significantly promoted international coordination in SDGs monitoring and favored production, access and usability of global and multi-temporal HRLC data products. Accordingly, it has been assessed that currently available open and global HRLC data could potentially be exclusively used for measuring indicators from four of the SDGs (i.e. 6 Clean water and sanitation, 13 Climate action, 14 Life below water, 15 Life on land), and could be used to complement other data types for four other goals (2 Zero hunger, 9 Industry, innovation and infrastructure, 11 Sustainable cities and communities, 12 Responsible consumption and production) [2]. A comprehensive list of SDG indicators to which global HRLC data as well as geospatial information in general have a direct contribution is provided here.
This enhanced global HRLC data availability has led to a wealth of studies, applications and methodology guidelines connected to the use of such data in SDG indicators computing and monitoring [3]. Relevant examples in the literature include the analysis of land-cover efficiency [4], land consumption rate, and forest cover [5], which are key thematic variables e.g. for SDG 15 and 11.
Open and global HRLC data provide different temporal coverages, spatial resolutions, thematic details (or classes), and accuracies. Therefore, the suitability (fit-for-purpose) of global HRLC data for each specific SDG indicator has to be assessed case by case. According to the literature [6] fit-for-purpose global HRLC data should be characterized by a minimum spatial resolution of 1 km and a minimum temporal resolution of 3 to 5 years. The thematic detail is instead indicator-specific. Finally, the accuracy of global HRLC data generally varies across the globe and local validation is strongly suggested depending on study areas and the user-expected accuracy performances.
Some of the most popular and mature, openly available, fit-for-purpose global HRLC datasets are reported and detailed in Table 2.1.1.
Global HRLC Dataset |
Temporal Coverage |
Spatial Resolution |
Thematic detail |
Notes |
|---|---|---|---|---|
GlobeLand30 |
2000, 2010, 2020 |
30 m |
10 main classes: Cultivated land, Forest, Grassland, Shrubland, |
-Provider: National Geomatics Center of China |
Finer Resolution |
2010, 2015, 2017 |
30 m |
10 main classes: Cropland, Forest, Grass, Shrub, Wetland, |
-Provider: University of Tsinghua |
Copernicus |
2015-2019 |
100 m |
10 main classes (cover fractions): Forests, Shrubland, |
-Provider: Copernicus Land Monitoring Service |
ESA-CCI-LC |
1992-present |
300 m |
22 classes aligned with UN-FAO’s Land Cover |
-Provider: ESA Climate Change Initiative |
The Terra and |
2001-present |
500 m |
Multiple classes and classification schemas |
-Provider: United States |
Global Human |
1975-2030 |
Up to 100 m |
Built-up |
-Provider: European |
Freshwater |
2000-2018 |
Up to 300 m |
Permanent and seasonal surface water, |
-Provider: UN Environment |
Footnotes
Other relevant global geospatial datasets for SDGs
Global HRLC data are also used in combination with complementary global geospatial information demonstrating their capabilities of directly supporting the computation of indicators e.g. from SDG 2 [1], 6 [2], and 9 [3].
Actually, land cover is only one of the Global Fundamental Geospatial Data Themes to SDGs identified by the United Nations. These Fundamental Themes include a variety of geospatial datasets which are necessary to spatially represent key phenomena and objects for the realization of economic, social, and environmental targets identified by the SDGs. These datasets include - among others - demography, soil/vegetation parameters, administrative boundaries, infrastructures, transport networks, 2D/3D topographic maps, water bodies, etc.
The suitability (fit-for-purpose) of complementary global geospatial datasets is difficult to assess globally for each specific SDG indicator or country. Spatial and temporal resolutions, thematic details and accuracies should be evaluated case-by-case, using recommendations valid also for the global HRLC data, and their use should be functional to the provision of baseline data in support of national SDGs reporting.
Table 2.2.1 contains information on popular global open geospatial datasets which might be employed in combination with global HRLC data for the computation of land cover-connected SDG indicators.
Global HRLC Dataset |
Temporal Coverage |
Spatial Resolution |
Thematic detail |
Notes |
|---|---|---|---|---|
World Bank Official |
Annual |
/ |
Administrative boundaries (and polygons) including |
-Provider: World Bank |
WorldPop Population Counts |
2000-2020 |
Up to 100 m |
Population |
-Provider: WorldPop - University of Southampton |
Global Human Settlement population |
1975-2030 |
Up to 100 m |
Population |
-Provider: European Commission, |
GHS Urban Centre |
2015, 2019 |
1 km |
Urban centers polygons |
-Provider: European Commission, |
OpenStreetMap |
/ |
/ |
Multiple (including transport networks) |
-Provider: OpenStreetMap |
Footnotes
Tools
QGIS
Among the available free and open-source GIS software, QGIS is the most popular and widely used desktop GIS application. The development of QGIS (formerly Quantum GIS) started in early 2002 and it is still ongoing thanks to the collaborative effort of a large community of developers and users, under the advocacy of the Open Source Geospatial Foundation (OSGeo). QGIS runs on Unix, MacOSX, and Microsoft Windows operating systems. The software is periodically updated through regular releases and bug fixes from volunteer developers coordinated by the QGIS project board. QGIS provides a user-friendly graphical user interface (Fig. 3.1.1) and it supports a variety of data types and formats as well as connection to online GIS resources through OGC service. Many functionalities for map creation and spatial data analysis are integrated into QGIS and can be accessed directly from the software interface. The possibility of accessing the source code allows for adapting the software as needed for specific GIS tasks. Furthermore, QGIS integrates with other free and open-source GIS packages that enable users to extend the capabilities of the software also without any programming skill. QGIS software can be downloaded from the official QGIS website by selecting among the available installation packages the one for your operative system. It is always suggested to work with the Long Term Release (LTR) version of the software which is the current stable and best-maintained release of QGIS.
QGIS main graphical user interface with logo
Important
The reference version for this webbook is QGIS 3.22. You can download and install QGIS on your computer to develop the hands-on exercises proposed in the next section (the QGIS LTR version is suggested).
Key QGIS functionalities for SDG indicators computation
QGIS provides thousands of functionalities and plugins to manipulate and analyze geospatial data. For the computation of SDG indicators using global geospatial datasets, a few key functionalities are needed. A selection of the most relevant ones is reported in Table 3.2.1. The full set of QGIS functionalities can be found in the official QGIS User Guide.
Function |
Description |
Link |
|---|---|---|
Vector and raster data import, |
Basic QGIS operations for data I/O |
|
Vector feature selections |
Selection/extraction of data |
|
Vector attribute table manipulation |
Editing tools for vector layers |
https://docs.qgis.org/3.22/en/docs/user_manual/working_with_vector/attribute_table.html |
Vector geoprocessing |
Geometry processing functionalities |
|
Raster layer clip |
Subsetting functionalities for |
|
Raster calculator |
Computational tool for |
|
Raster layer statistics |
Basic statistics extraction from |
|
Zonal statistics |
Raster statistics spatially constrained |
Case studies
Forest area as a proportion of total land area [15.1.1]
In this exercise you will learn how to calculate the indicator 15.1.1 which is defined as the forest area as a proportion of total land area, and is expressed in percentages. To learn more about this indicator read the metadata provided under this link. The method of computation is given by the following formula:
\(\frac{\text{Forest area (reference year)}}{\text{Land area (reference year)}} * 100\),
thus we need to first define the region and reference year for which the indicator needs to be calculated. After that, we need to estimate the forest area and land area of the region. In this case study we’ll be working on calculating the indicator for Suriname for the reference year 2019.
The first step is to download the vector data of the countries in shapefile format from the WorldBank (Fig. 4.1.1). The file contains polygons of all the world’s countries, and it will come useful in next exercises.
Downloading page of the WorldBank
After downloading the file open QGIS and create a new project, to which add the new shapefile (Fig. 4.1.2).
Adding a vector layer to QGIS interface
Now we need to extract from the world’s map only the polygon of Suriname. To do so we will use the “Select Features by Area or Single Click” and select Suriname as shown in Fig. 4.1.3.
Selecting Suriname by area
After selecting Suriname we will want to create a new shapefile containing only this one polygon. To do so we need to Save the selected features in our working directory and name it (Fig. 4.1.4). Make sure that the “Add saved file to map” box is checked.
Exporting selected features as a new layer
To better visualize the borders of Suriname, change the symbology of the layer. To do so, right click the “Suriname” layer in the layer pane and select “Properties”. After that go to “Symbology” and change it so that the fill color is set to “no brush”, the stroke color is of a bright, contrasting color, and the strike width is quite thick (Fig. 4.1.5).
Changing the symbology of the new Suriname.shp layer
Now that we have the Suriname polygon, which will serve us to derive the land area for the calculations, it’s time to access the forest data. We will use the Global Land Cover data from https://lcviewer.vito.be/2015. After accessing the website follow the instructions shown in Fig. 4.1.6 to retrieve the forest coverage raster file of the Suriname area.
Downloading forest data
After downloading the file, add it to the project by using the “Add Raster Layer” function, which you’ll access analogously to the “Add Vector Layer” function presented before.
Adding the forest raster data layer to the project
After adding the downloaded vector and the raster layer, the view of your project should be similar to the one shown in Fig. 4.1.8.
Expected layer view after adding both the vector and raster layer
To extract the forest areas of Suriname we need to clip the forest raster layer by mask as shown in Fig. 4.1.9.
Clip forest data raster layer by Suriname vector mask layer
The forest raster layer is a Fractional Cover layer, this means that it gives the percentage of a 100 m pixel that is filled with forest. Pixels have values between 0 and 100, in steps of 1%. Since we just need the information whether the forest is present or not, we will reclassify the layer using the Raster Calculator (Fig. 4.1.10). We’ll classify the pixels with a value equal or bigger than 50% as forest (1) and the rest as no forest pixels (0).
Reclassify the clipped raster layer
After reclassifying the raster let’s change the symbology of the new layer for better visualization of the forest area. The expected view after this operation should be as shown in Fig. 4.1.12.
Changing the symbology of forest data reclassified raster layer
Expected layer view after changing the reclassified layer symbology
Since the pixels representing the forest have the value of 1 the sum of all the pixel values will give us the number of forest pixels in the raster file, and knowing the pixel’s size we can calculate the area of the forest in Suriname. To calculate the sum of the pixels values we”ll use the “Zonal Statistics” tool (Fig. 4.1.13).
Zonal statistics tool
To organize the attribute table we can delete some of the unnecessary columns in the Suriname statistics layer (Fig. 4.1.14). It is not a mandatory step, but it provides a more clear view of the table and makes it easier to work with. The “_sum” column is the result of applying the “Zonal statistics” tool in the previous step.
Delete the unwanted columns
The last step is to finally calculate the indicator. We have all the required data in the “Suriname_stats” layer. We will use the field calculator to calculate the indicator. Since we now that the pixel size is 100 m x 100 m we will calculate the forest area as: “_sum” * 100 * 100. The land area is given by the built in function “$area”. The exact formula and result is shown in Fig. 4.1.15.
Calculate the indicator 15.1.1 using the field calculator
The final result suggests that the forest areas contributed approximately 77.42% of the total land area of Suriname in 2019.
Ratio of land consumption rate to population growth rate [11.3.1]
In this section you will learn how to calculate the SDG indicator 11.3.1 which is defined as the ratio of land consumption rate to population growth rate. To be able to derive this indicator we need to first define the population growth and the land consumption rate.
The population growth rate (PGR) is the change of a population in a defined area during a period, expressed as the percentage of the population at the start of that period. Its value depends on the number of births and deaths and on the migration patterns in the given period. In SDG 11.3.1 the PGR is computed for city areas. The land consumption rate (LCR) is defined as the rate at which urbanized land changes through a defined period of time. The LCR is expressed as the percent of the land occupied by the city/urban area at the start of that time. SDG 11.3.defines “land consumption” as the conversion of land from non-urban to urban functions. For the correct calculation of the index it is also needed to agree on what constitutes a city or urban area. For this purpose we will be using the Degree of Urbanization (DEGURBA) method to delineate cities, urban and rural areas for international statistical comparisons endorsed by the United Nations Statistical Commission For more information about the background of the SDG indicator 11.3.1 read the metadata provided under this link The workflow for computing the indicator follows 5 steps:
Decide on the analysis period
Delimitation of the urban area which we will analyze
Computation of the land consumption rate
Computation of the population growth rate
Computation of the ratio of LCR to PGR
In this exercise we will work on a one year interval from 2018 to 2019.
We will use New Delhi as the analysis area. Download from https://ghsl.jrc.ec.europa.eu/ghs_stat_ucdb2015mt_r2019a.php the delimitation of urban areas file as shown in (Fig. 4.2.1).
Downloading the delimitation of urban areas file
Open a new QGIS project and add the downloaded vector (.gpkg) file to it (Fig. 4.2.2).
Adding the vector layer to a new QGIS project
We now want to extract only the delineation of New Delhi. After opening the attribute table of the just added layer, click on the “Select features using an expression” icon (Fig. 4.2.3).
Select features using an expression button in the attribute table
After the Selection by Expression window appears, enter the query: UC_NM_MN = "Delhi [New Delhi]" and click on “Select Features” (Fig. 4.2.4).
Selecting features by expression
After selecting New Delhi, extract it to a new vector layer “NewDelhi.shp” as shown in Fig. 4.2.5.
Saving selected features as new vector layer
Change the symbology of the “NewDelhi.shp” layer to better visualize the delineation of the urban area (Fig. 4.2.6).
Symbology of the New Delhi shapefile vector layer
To compute the land consumption rate in the period 2018-2019 we need the urban areas data from both of the year. We can access the Land Cover raster data from https://lcviewer.vito.be/download. Download the indicated tile from 2018 and 2019 as shown in Fig. 4.2.7.
Land cover raster data download
After downloading both raster layers, add them to your QGIS project (Fig. 4.2.8).
Adding the raster layers into the QGIS project (repeat for the 2019 layer)
Now, clip both layers by mask of the New Delhi boundary, so that we will work only on the area of interest (Fig. 4.2.9).
Clipping the raster layer by the New Delhi boundary mask (repeat for the 2019 layer)
Since we need only the urban areas of New Delhi, we would like to extract from the layers only the pixels representing urban/built up areas. From the Copernicus Global Land Service: Land Cover 100m: version 3 Globe 2015-2019: Product User Manual we know that the map code for the urban/built up class is 50 (Fig. 4.2.10).
Classification of the land cover layer with map codes
We need to reclassify both rasters so that the pixels with the value of 50 have the value of 1 and the rest the value of 0. This step is performed by using the “Raster Calculator” and inputting the expression: (“LC_New_Delhi201x@1” = 50) * 1 + (“LC_New_Delhi201x@1” != 50) * 0, where x is equal to 8 or 9 depending on the layer. The procedure for the 2018 layer is shown in Fig. 4.2.11.
Raster reclassification (repeat for the 2019 layer)
For better visualization, change the symbology of both reclassified layers, so that the render type is set to “Paletted/Unique values” and the color of class 0 is set to transparent, as shown in Fig. 4.2.12.
Symbology properties for the reclassified land cover raster layers (repeat for the 2019 layer)
The desired view after clipping, reclassifying and changing the symbology of the LC layer is as in Fig. 4.2.13.
Land cover layer of the urban areas after the preprocessing steps
To be able to calculate the Land Consumption Rate, we must know the area of the built up zones in both years. To do so, we firstly calculate the sum of the pixels by using the “Zonal Statistics” tool (Fig. 4.2.14). By calculating the sum of all pixels we will actually get the sum of the pixels representing the urban class, as they have value 1 and the rest has value 0.
Calculating the sum of pixel values using “Zonal statistics”
Having the zonal_stats layers with the sum of the urban pixels, we can now calculate the area of the urban zones in New Delhi in both years, knowing that the pixel size is 100x100m. For both zonal_stats layers open the field calculator and add a new field for the built up area in square km. Calculate the value of the field by multiplying the “_sum” field by 0.1 * 0.1 (km). The step for the “zonal_stats_2018” is shown in Fig. 4.2.15, be sure to repeat the step also for “zonal_stats_2019”.
Calculating the urban area in square kilometers (repeat for the 2019 layer)
Finally, it is possible to calculate the Land Consumption Rate (LCR) by applying this formula: LCR = (Vpresent - Vpast)/Vpast, to the calculated values in the previous steps. We will calculate it in a new field in the attribute table of the “zonal_stats_2018” layer as shown in Fig. 4.2.16.
Calculating the Land Consumption Rate
Having calculated the LCR, it is time to calculate the second index needed to calculate SDG 11.3.1 - the Population Growth Rate.
The PGR is calculated using the total population within the urban area for the analysis period using the formula below:
\(PGR = \frac{\ln(\frac{Pop_{(t+n)}}{Pop_t})}{y}\),
where:
\(\ln\) is the natural logarithm value;
\(Pop_t\) is the total population within the urban area in the initial year;
\(Pop_{t+n}\) is the total population within the urban area in the final year;
\(y\) is the number of years between the two measurement periods.
Thus, to calculate the PGR we need the population of New Delhi in 2018 and in 2019. To get this information we first need to download the population layers for India in 2018 and 2019 from https://hub.worldpop.org/geodata/listing?id=29 by clicking “Data & Resources” and then “Download Entire Dataset” (Fig. 4.2.17).
Downloading the population layers for India
Load the population layers to the QGIS project and clip them both to New Delhi boundaries. (Fig. 4.2.18).
Clipping the population layer for the year 2018 to New Delhi’s boundaries. Repeat the process for the 2019 layer
Use the “Zonal Statistics” tool to calculate New Delhi’s population in 2018 and 2019. Select only the “Sum” statistics, which will give us the population, as shown in Fig. 4.2.19.
Calculating New Delhi’s population in 2018. Repeat for the 2019 layer
After the previous step we’ll obtain two layers for both years. To make the computations easier we want to have the population field for both years in one layer. To do so add a new column “population_2019” to the 2018 population zonal statistics layer. Then copy the content of the “_sum” field from the 2019 population zonal statistics layer to the newly created column in the 2018 layer (Fig. 4.2.20).
Adding the population of 2019 to the 2018 zonal statistics layer
To change the name of a field in the attribute table go to the “Properties” of the layer and in the “Field” section activate the editing option (pencil icon). Now you can change the name of the “_sum” column to “population_2018” (Fig. 4.2.21).
Changing the name of a field in a layer’s attribute table
Now that we have the population values for both years in one layer we can calculate the PGR using the “Field Calculator” tool. Create a new field where you’ll calculate the PGR value using the before mentioned formula as shown in Fig. 4.2.22.
Calculating the PGR index
To have both the PGR and the LCR indexes in one layer, add a new field named “PGR” in the attribute table of the “zonal_stats_2018” layer containing the previously calculated LCR (Fig. 4.2.23).
Adding the PGR field into the “zonal_stats_2018” layer
Having the LCR and the PGR in the “zonal_stats_2018” layer we can calculate the SDG indicator 11.3.1 in the same layer by dividing the LCR by the PGR.
To calculate the indicator we will use once again the field calculator as shown in Fig. 4.2.24.
Calculating the SDG indicator 11.3.1
For New Delhi the 11.3.1 indicator was estimated to be around 0.0387 in the 2018-2019 period.
Proportion of the rural population who live within 2 km of an all-season road [9.1.1]
This exercise’s purpose is to show how to compute the SDG indicator 9.1.1 defined as the share of a country’s rural population living within 2 km of an all-season road. This indicator is also commonly known as the Rural Access Index (RAI). To calculate the indicator it is needed to compute the ratio between the rural population within a 2 km buffer of an all-season road and the total rural population of the country. To compute this indicator we need three different geospatial datasets:
The country’s population distribution, that can be retrieved from WorldPop;
From which it is necessary to extract only the rural population, which can be done using a Rural-Urban global definition provided by the Global Human Settlement Layer;
The road network with road types attribute, which can be accessed for all countries from OpenStreetMap;
For more information about the background of the SDG indicator 9.1.1 refer to the metadata provided here link.
In this exercise we will calculate the indicator for Gabon in 2020.
First, we need to download all the necessary data:
From https://hub.worldpop.org/geodata/listing?id=78 download the population raster layer for Gabon in 2020 as shown in Fig. 4.3.1.
From https://ghsl.jrc.ec.europa.eu/ghs_stat_ucdb2015mt_r2019a.php download the Urban Centre Database (Fig. 4.3.2)
Downloading the population raster layer of Gabon 2020
Downloading the delimitation of urban areas file
From https://download.geofabrik.de/africa/gabon.html download the OSM data. After that unpack the .zip file and keep only the files related to the road network (”gis_osm_roads_free”). The procedure is shown in Fig. 4.3.3
Downloading the road network data from OSM.
Open a new QGIS project and add the road network shapefile (Fig. 4.3.4) and the downloaded population grid file (Fig. 4.3.5).
Adding the vector layer of the Gabon road network to the QGIS project
Adding the Gabon population raster layer to the QGIS project
Since we are interested only in all-season, major roads we need to extract only those from the shapefile vector layer. The information about the codes for the roads from the “Format Specification” for OSM data (Fig. 4.3.6) indicates that the major roads have codes starting with 511.
“Format Specification” regarding OSM roads’ codes
Using the selection by expression, we can select the roads with codes that are smaller than 5120 as shown in Fig. 4.3.7. After the selection of major roads we want to export them in a new vector layer that will be called “all_season_roads.shp”.
Selecting major roads in Gabon
The procedure of extracting selected features is presented in Fig. 4.3.8, and the expected view after the procedure is shown in Fig. 4.3.9.
Extracting the selected features
Expected view after extracting only Gabon’s major roads.
We now need to create the 2 km buffer around the major roads, by using the “Buffer” geoprocessing tool (Fig. 4.3.10). Because of reprojecting issues that may be encountered we input the distance in degrees.
Warning
To work with metric buffer distance units, reproject all the layers in the correct UTM zone.
Creating a 2 km buffer around the major roads
Now we have the buffer around the major roads and the population grid of Gabon. We want to extract just the rural population within that buffer, hence we need to exclude the urban areas from the computations. To do so add the urban extent layer previously downloaded (Fig. 4.3.11) and compute the difference between the previously created road buffer and the urban extent layer (Fig. 4.3.12). This will create a new vector layer containing the buffer only in rural areas.
Adding the urban extent vector layer to the QGIS project
Subtracting the urban extent layer from the road buffer layer
Warning
In case of “Invalid Geometry” error: click the wrench icon by the side of the overlay layer and select “Do not Filter (Better Performance)” from the drop down menu in the “Invalid feature filtering” option.
Having this vector layer and the population grid layer we will use them as input layers to calculate the rural population within the 2 km buffer of a major road by using the “Zonal Statistics” tool (Fig. 4.3.13).
Calculating the rural population within the 2 km buffer around major roads.
After this step we have the rural population within a 2 km buffer of a major road (Fig. 4.3.14).
Restult of the “zonal statistics” operation
To calculate the indicator 9.1.1 we also need the total rural population of Gabon, which we will calculate in the next steps. Firstly, we need Gabon’s boundaries to delimitate the urban extent layer ust to our area of interest. To retrieve this data go to https://datacatalog.worldbank.org/search/dataset/0038272 and download the “World Country Polygons - Very High Definition” as shown in Fig. 4.3.15.
Downloading the “World Country Polygons”
Add the downloaded shapefile to the QGIS project (Fig. 4.3.16).
Adding the countries’ boundaries vector layer to the QGIS project
Since we are interested only in the boundaries of Gabon we need to extract them from the vector file. To do so select Gabon’s polygon by area as presented in Fig. 4.3.17.
Selecting Gabon by area
Now extract the selected feature to a new vector layer “Gabon.shp” (Fig. 4.3.18).
Extracting the selected Gabon polygon
For better visualization purposes change the symbology of the new layer in its “Properties”, which can be accessed by right clicking the layer in the layers panel. An example of a clear symbology for a country’s boundaries is presented in Fig. 4.3.19.
Changing the “Gabon.shp” symbology
Now we can clip the urban extent layer to Gabon’s borders by using the “Clip” geoprocessing tool (Fig. 4.3.20).
Clipping the urban extent layer to Gabon’s extent
The urban extent polygons are now limited only to Gabon’s extent. If we open the attribute table now, we can see that for each urban area there’s a distinct row. Since we are interested in the urban areas as a whole we want to have just one record for all the urban areas in Gabon. To do so we need to use the “Dissolve” geoprocessing tool as shown in Fig. 4.3.21.
Dissolving the urban extent layer
The attribute table of the Gabon’s urban extent layer after this operation can be found in Fig. 4.3.22.
The expected result of the “Dissolve” procedure
We can now easily calculate the total urban population with “Zonal Statistics”, in the newly created vector layer the “_sum” field presents the total urban population of Gabon (Fig. 4.3.23).
Calculating the total urban population of Gabon
To calculate the needed rural population of Gabon it is needed to first retrieve the total population, so then we can compute:
\(\text{Rural population} = \text{Total Population} - \text{Urban Population}\);
To calculate from the population grid the total population we will use the “Raster layer statistics” with the population raster layer as input. The output of this process is stored in a temporary .html file which you can find in the right bottom corner of your screen in the “Result Viewer” panel (Fig. 4.3.24).
Calculating the total population of Gabon in 2020 with the “Raster layer statistics”
We now have all the necessary data to calculate the indicator. To make the computations easier and faster we first need to add all the needed values into one layer. Open the attribute table of the “zonal_stats_buffer.shp” layer, which we previously created to calculate the rural population within the 2 km buffer. Start editing and add two new integer fields, one for the total population and the second for the urban population (Fig. 4.3.25).
Adding new fields to the attribute table of “zonal_stats_buffer.shp”
To the “total_pop” field paste the value of the “Sum” from the .html file generated by the “Raster layer statistics” tool (Fig. 4.3.26).
Populating the “total_pop” field
To the “urban_pop” paste the value of the “_sum” field from the attribute table of the “urban_zonal_stats.shp” layer (Fig. 4.3.27).
Populating the “urban_pop” field
Finally, we can calculate the SDG 9.1.1 indicator using the “Field Calculator”. Since the indicator is defined as:
\(\frac{\text{Rural population living in a 2 km buffer from a major road}}{\text{Total rural population}}\)
we can compute it as shown in Fig. 4.3.28.
Calculating of the SDG 9.1.1 indicator
The final result (Fig. 4.3.29) indicated that around 71% of Gabon’s total rural population lives in a 2 km buffer from a major, all-season road.
The attribute table with the final result
Change in the extent of water-related ecosystems over time [6.6.1]
This exercise introduces a possible workflow to calculate the SDG indicator 6.6.1, which tracks how different types of water-related ecosystems change in extent over time. The indicator captures data on different types of freshwater ecosystems and measures their extent change. The indicator 6.6.1 is complex and it focuses on more aspects than just the change in the spatial extent of permanent freshwater, but for the exercise’s simplicity sake just this unit of measure will be presented in this section. This index will be presented by the percentage change (∆) in the area of permanent waters. We will calculate ∆ as:
\(∆ = \frac{(α - β) + (ρ - σ)}{ε + β + σ}\), where:
α - new permanent water (conversion of a no water place into a permanent water place)
β - lost permanent water (conversion of a permanent water place into a no water place)
ρ - seasonal to permanent (conversion of seasonal water into permanent water)
σ - permanent to seasonal (conversion of permanent water into seasonal water)
ε - permanent water surfaces (where permanent water is always observed)
The numerator represents the change trends of the permanent water (if the numerator is positive it means that the area of permanent water places increased with respect to the beginning of the 5 year measurement period). The denominator represents the permanent water surfaces at the beginning of the considered measurement period.
For more information about the background of the SDG indicator 6.6.1 refer to the metadata provided under this link. In this exercise we will focus on the water-related ecosystems in Egypt in a five year period from 2014 to 2018.
To calculate this change index we need first to download the data for:
Coutries border delimitation from https://datacatalog.worldbank.org/search/dataset/0038272 (Fig. 4.4.1)
Downloading the coutries borders layer
Egypt’s water surfaces from “Freshwater Ecosystems Explorer” for the years 2014 and 2018 (Fig. 4.4.2), in which we can find a layer with one band with three possible values:
1 - Not Water;
2 - Seasonal Water;
3 - Permanent Water.
After downloading the file remember to unzip it before working on it.
Downloading the Egypt’s water surfaces layers
Open a new QGIS project and add the countries borders vector layer to it as in Fig. 4.4.3.
Adding the vector layer to the QGIS project
Since we are interested only in Egypt’s water areas, we need to extract its polygon by first selecting it (Fig. 4.4.4) and then extracting the selected features and saving them in a new layer (Fig. 4.4.5).
Selecting Egypt’s polygon
Extracting the selected polygon as a new layer
For better visualization purposes change the symbology on the newly added Egypt layer as portrayed in Fig. 4.4.6.
Changing the symbology of the “Egypt.shp” layer
Now that we have defined the Area Of Interest (AOI) in our project we can add the water extent data (Fig. 4.4.7).
Adding the raster layer of the water extent of Egypt in 2018. Repeat for the 2014 layer
As mentioned before each year’s file has three possible values as specified in Fig. 4.4.8, so we want to represent these classes by changing the symbology of both layers. A proposed symbology for the 2018 layer is presented in Fig. 4.4.9. Make sure to apply the symbology change to both layers with the same palette.
Specifications for the raster water extent layers
Changing the symbology of the water extent layers for better visualization purposes
After adding and changing the symbology of both needed datasets the desired view of the project should be as shown in Fig. 4.4.10.
The desired view after adding and changing the symbology of the data
Now we need to extract from the two raster layers the components α, β, ρ, σ, ε to calculate the change index as defined in Fig. 4.4.11.
The formula of the index to be calculated in this exercise
To do so we will use the “Raster calculator” to create a new raster layer with new values (α, β, ρ, σ, ε) by applying a reclassification scheme presented in Fig. 4.4.12, remembering that the values of the rasters are: 1 - Not Water; 2 - Seasonal Water; 3 - Permanent Water.
The expression “Egypt_classes_2014@1 = 1 AND Egypt_classes_2018@1 = 3” represents the change from a no water place to a water place, which corresponds to α in the equation defined above. After clicking “Run” the process may take a while to complete.
Creating a new raster layer with reclassified values according to the needed information for the index calculation
After the process is done the values of the new raster layer will be: * 1 → α; * 2 → β; * 3 → ρ; * 4 → σ; * 5 → ε.
For better visualization change the symbology of the new raster as shown in Fig. 4.4.13.
Changing the symbology of the new raster
Now that we have a raster showing the changes in permanent water extent we want to know what is their area. To do so we chose the approach of vectorizing the transformation layer for easier computations later on. The process of transforming tre transformation raster layer into a vector layer is heavy and will take a considerable amount of time (YOU CAN DOWNLOAD THE ALREADY VECTORIZED LAYER FROM THE WEBBOOK ZENODO FOLDER). To compute the vector layer we use the “Polygonize (Raster to Vector)” tool from the conversion raster tools (Fig. 4.4.14).
Polygonizing the transformation raster
This vector layer also contains land features that we are not interested in, hence we will extract all the features with values not equal to 0. This will be done by first selecting all such features by expression and then exporting the selected features to a new layer. The workflow is presented in Fig. 4.4.15.
Extracting only the information about water areas from the vector
The layer’s attribute table contains all the polygons as separate. We want to group them by their classification so that the attribute table will contain just five classes representing the component of the final equation. To do that we will use the “Dissolve” geoprocessing tool. It is important to remember what is the field that we are dissolving and set it in the optional settings of the process. The procedure is presented in Fig. 4.4.16.
Dissolving the features of the vector according to the class field (DN)
Now that we have all of our component’s classes grouped we can calculate each class’ area to have the final values for the computation. We will use the “Add Geometry Attributes” vector geometry tool using as input the dissolved transformation vector containing water information. The ellipsoidal parameter should be set in the calculation parameter (Fig. 4.4.17).
Adding geometry attributes to the vector layer
The output vector layer will now contain the area for each type of water transformation (Fig. 4.4.18).
The final view of the attribute table of the vector layer
The values can be plugged in directly to the equation to get the final result of the exercise as the percentage change in spatial extent of permanent water. Remember that in our case: 1 → α; 2 → β; 3 → ρ; 4 → σ; 5 → ε. The area change of permanent water component of the SDG indicator 6.6.1 is equal to approximately 3%, which indicates that in the 2014-2018 period there was an increase (because it’s positive) of permanent water areas of 3%.
Note
Input data for the presented case studies can be downloaded at https://doi.org/10.5281/zenodo.7152401.
Warning
The numerical results of the presented case studies are to be intended as purely demonstrative. Revision and validations of the outlined procedures and results may be required for their reuse outside the domain of this webbook.
Credits
Dr Daniele Oxoli
Daniele Oxoli was born in Como, Italy, in 1990. He received the B.Sc. degree in civil and environmental engineering in 2013, and the M.Sc. degree in geomatics engineering in 2015 from the Politecnico di Milano, Milan, Italy. In 2019, he obtained his Ph.D. degree with honors in Geomatics Engineering from the Politecnico di Milano. He is currently assistant professor at the Geomatics and Earth Observation laboratory (GEOlab) of Politecnico di Milano and he is involved in a number of research projects connected to the use and development of Free and Open Source GIS software and the statistical analysis of spatial data. He is also a Charter Member of the Open Source Geospatial Foundation (OSGeo) and Secretary of the ISPRS WG IV/7 “Intelligent Systems in Sensor Web and IoT”. His research interests include spatial data analysis, spatial statistics, data science and visualization, and Earth Observation.
Dr Sheryl Rose C. Reyes
Sheryl Rose C. Reyes is a Remote Sensing Expert at the United Nations Satellite Centre (UNOSAT) under the United Nations Institute for Training and Research (UNITAR). She received the M.Sc. degree in remote sensing and the B.Sc. degree in geodetic engineering from the University of the Philippines, Diliman. She was also one of the lead authors for the Global Environment Outlook 6 for Youth (GEO-6 for Youth) report, which is the first fully interactive e-publication of the United Nations Environment Programme (UNEP), written by youth for youth to inform, engage, educate, and lead youth towards environmental action. In addition to her professional experience, she was also the former President of the International Society for Photogrammetry and Remote Sensing Student Consortium (ISPRS SC) from 2016 - 2022, the representation of the youth and students in ISPRS. Her research interests include remote sensing, land cover change modeling, GIS and ecosystem services.
Dr Shu Peng
National Geomatics Center of China, Beijing
Prof. Maria Antonia Brovelli
Degree with honors in Physics, PhD in Geodesy and Cartography. She is Professor of “GIS” and “The Copernicus Green Revolution for sustainable development” at Politecnico di Milano (PoliMI) and a member of the School of Doctoral Studies in Data Science at “Roma La Sapienza” University. From 2006 to 2011 she lectured GIS at the ETH of Zurich and from 1997 to 2011 she was the Head of the Geomatics Laboratory of PoliMI (Campus Como). From 2011 to 2016 she was the Vice-Rector of PoliMI for the Como Campus. Currently, she is the coordinator of the Copernicus Academy Network for the PoliMI and the Head of the GEOLab, the Interdepartmental Lab where 7 Departments of POLIMI are contributing. She is Vice President of the ISPRS Technical Commission on Spatial Information Science, former member of ESA ACEO (Advisory Committee of Earth Observation); co-chair of the United Nations Open GIS Initiative, chair of the UN-GGIM (Global Geospatial Information Management) Academic Network, mentor of the PoliMI Chapter of YouthMappers (PoliMappers), one of the three curators of the geospatial series of the AI for Good, organized by ITU in partnership with 40 UN Sister Agencies. Her research activity is in the field of geomatics. Her interests have been various, starting from geodesy, radar-altimetry and moving later to GIS, webGIS, geospatial web platform, VGI, Citizen Science, Big Geo Data, geoAI. She is participating and leading research on these topics within the frameworks of both national and international projects and scientific networks. One of her main interests is in Open-Source GIS, where she is playing a worldwide leading role.
Prof. Serena Coetzee
University of Pretoria, South Africa
Dr Ivana Ivánová
Dr Ivana Ivánová is a Senior Lecturer in Spatial Sciences and FrontierSI Research Fellow in Spatial Information Infrastructures. Her research interests and expertise are in spatial data quality spatial resources and spatial information infrastructures. Dr Ivánová holds an engineering and doctoral degree, both from the Slovak University of Technology in Bratislava, in geodesy and cartography with specialization in Geoinformatics. She researched and lectured at several universities – Slovak University of Technology in Bratislava (2000-2007), University of Twente in the Netherlands (2007-2013) and São Paulo State University in Brazil (2014-2017). Dr Ivánová has extensive experience in standardization. From 2004 - 2007, she represented the Slovak national standardization organization in the CEN/TC 287 Geographic information’s Outreach Group. She has been overseeing and reviewing the adaptation of ISO 19100 series of norms into a national legal framework for geographic information. She currently leads development of ISO 19157-1 and ISO 19157-3 standards on geographic information quality. At ISO/TC211 Dr Ivánová is co-convening the Group on Ontology Maintenance and at Open Geospatial Consortium (OGC) she co-chairs the Data Quality Domain Working Group. Ivana represents Curtin University in the IT-004 Geographic information/Geomatics working group at Standards Australia. She regularly contributes to the work of several geospatial (incl. OSGeo, Geo4all, FOSS4G) and research (incl. RDA, ESIP and ARDC) communities.
Dr Darshana Rawal
Has done PhD in Geography, M.Sc. in Geomatics and Master of computer Applications. She has 16 years of teaching and research experience. Her area of specialization is ‘Information Science and Geospatial Technology’. The specific application areas are Disaster, Health, Sustainable Environment, Urban Ecosystems that include City Planning, Infrastructure, Environment; IT-based web development, which are the key subjects of my expertise. She teaches various subjects named GIS, Web-GIS, Crowd Sourcing and Participatory-GIS, and Digital Image Processing. She is a member of various international Societies and organizations, those are: Member of the GeoforAll 2023; Co-Chair of WG V/3 Open Source Promotion and Web-based Resource Sharing; Member of the International Advisory Editorial Board of the Geo-Progress Journal; Member of the Scientific Council for the section “Methodological and Technical Issues of Geographic Information and Spatial Analysis; of GeoProgress Journal; Member, UN Open GIS Initiative of OSGeo, Member, FOSS4G an Academic/Scientific Track (AT) , First Prize for ‘Most Active Congress Contributor’ in XXIII ISPRS Congress, Prague, Czech Republic, 2016 and Excellence Award for a Project on’’ Use of MongoDB for Social Media Database Management” during International Geoinformatics Summer School at Wuhan University, China, 2015.
Prof. Giuseppina Vacca
University of Cagliari, Italy
Prof. Sisi Zlatanova
University of New South Wales, Built Environment, Red Centre Building, Australia
Julia Anna Leonardi
Julia was born in Warsaw, Poland, in 1998. She received her B.Sc. degree in Geography with specialization in geoinformatics from the University of Warsaw. Her B.Sc. thesis was on assessment of different Machine Learning algorithms on hyperspectral images for high mountain vegetation classification. She is now persuing a M.Sc. in Geoinformatics Engineering at Politecnico di Milano.
Contacts
Daniele Oxoli, PhD
Politecnico di Milano - GEOlab
Dept. of Civil and Environmental Engineering
Building 3a, P.zza Leonardo da Vinci 32, 20133 Milano (IT)
daniele.oxoli@polimi.it
Important
This project was funded by the International Society of Photogrammetry and Remote Sensing (ISPRS) within the Educational and Capacity Building Initiative 2022.
Note
All materials in this web-book are from public literature and the Internet. Unintended use herein of the copyrighted material is purely for educational purposes
License Creative Commons Attribution 4.0 License (CC BY 4.0)