Multi-criteria Overlay Analysis is the process of the allocation of land to suit a specific objective on the basis of a variety of attributes that the selected areas should possess. Although this is a common GIS operation, it is best performed in the raster space.

This post outlines the typical workflow to take source vector data, transform them to appropriate rasters, re-classify them and perform mathematical operations to do a weighted suitability analysis. The post uses the command-line utilities provided by the open-source GDAL library. If you want to do such analysis in QGIS, please check my tutorial at Multi Criteria Overlay Analysis (QGIS3)

We will work with crime and infrastructure data for the city of London and find suitable areas to build new parking facilities that can help reduce bicycle thefts. Our analysis will apply the following 3 criteria. The proposed parking must be

1. In a bicycle theft hotspot
2. Close to a bicycle route
3. Far from existing parking facilities

Want to follow along?

• Install GDAL. I prefer installing using conda.
• The source data layers are provided in the multicriteria.gpkg. Download the file and save it to your computer. Open a terminal and switch to the folder containing the geopackage.

Step 1: Rasterize vector layers

The step in the overlay analysis is to convert each vector data layer to raster. This is accomplished using the gdal_rasterize command. Let’s check the content of the data layers.

ogrinfo multicriteria.gpkg

An important consideration is that all rasters must be of the same extent. We can inspect use the boundary layer and find out the extent in the CRS of the layer. We will specify this extent using the -te option with gdal_rasterize command.

ogrinfo -so multicriteria.gpkg boundary

Let’s start by rasterizing the cycling_routes layer. We want to create an output raster where the pixel values are 1 where there is a route and 0 where there are no routes. This value is specified using the -burn option. We can set the resolution of the output raster to be 10 meters using the -tr option. Once the command finishes, each pixel’s value will be set to 1 or 0 depending on whether a route line intersects it or not.

gdal_rasterize -ot Int16 -burn 1 -tr 10 10 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l cycling_routes routes.tif

We can rasterize the cycle_parking layer with the same options.

gdal_rasterize -ot Int16 -burn 1 -tr 10 10 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l cycle_parking existing_parking.tif

Now only the bicycle_thefts layer need to be rasterized. This is a dense point layer with locations of individual thefts – with overlapping points. We can create a raster layer by adding the total number of points within each pixel and calculate the density of thefts. This type of approach in mapping point density is called Binning or Aggregation. We can use the -add option to increment the value to the pixel for each input feature.

gdal_rasterize -burn 1 -add -tr 100 100 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l bicycle_thefts thefts.tif -ot Int16

Step 2: Generate proximity (Euclidean distance) rasters

Now that we have all the input layers converted, we can generate proximity rasters. These are also known as Euclidean distance rasters – where each pixel in the output raster represents the distance to the nearest pixel in the input raster. The distance can be specified in pixels or geographic units using the -distunits option. This resulting raster can be then used to determine suitable areas which are within certain distance from the input.

gdal_proximity.py routes.tif routes_proximity.tif \
 -ot Int16 -distunits GEO

gdal_proximity.py existing_parking.tif \ 
 existing_parking_proximity.tif -ot Int16 -distunits GEO

Step 3: Re-classify raster values

To use the proximity rasters in overlay analysis, we must first re-classify it to create discrete values. Here we use the range 0-100 and assign an appropriate value based on the input range.

For point density raster representing Thefts, we assign higher value to more thefts. For Routes raster, we assign higher value to pixels that are close to the routes. Lastly, for the Existing Parking raster, we assigned higher value to pixels that are far away. The final assignment look like below.

The following commands uses gdal_calc.py to take the continuous values from the input rasters and create re-classified rasters with discrete values.

gdal_calc.py -A thefts.tif \
  --outfile thefts_reclass.tif \
  --calc="100*(A>20) + 50*(A>10)*(A<=20) + 10*(A<10)"

gdal_calc.py -A routes_proximity.tif \
  --outfile routes_reclass.tif \
  --calc="100*(A<=50) + 50*(A>50)*(A<=100) + 10*(A>100)"

gdal_calc.py -A existing_parking_proximity.tif \
  --outfile existing_parking_reclass.tif \
  --calc="100*(A>100) + 50*(A>50)*(A<=100) + 10*(A<50)"

Step 4: Overlay analysis

Now we are ready to do the final overlay analysis. Before we proceed, it is important to make sure that all the rasters have the same pixel resolution. We downsample the existing_parking_reclass.tif and routes_reclass.tif layers so their resolution matches thefts_reclass.tif layer.

gdalwarp -r average -tr 100 100 \ 
 existing_parking_reclass.tif parking_reclass_resampled.tif
gdalwarp -r average -tr 100 100 \ 
 routes_reclass.tif routes_reclass_resampled.tif

We are now ready to run the overlay analysis and generate the suitability map.

Simple OVERLAY

If we choose to give equal weight (importance) to each criteria – we can we can simply add all 3 rasters and normalize the pixel values so they are in the 0-100 range. We set the NoData value in the output using the --NoDataValue option.

gdal_calc.py \
  -A thefts_reclass.tif \
  -B routes_reclass_resampled.tif \
  -C parking_reclass_resampled.tif --outfile suitability.tif \
  --calc="(A + B + C)/3" --NoDataValue=0

Weighted Overlay

In a weighted-overlay analysis different criteria are assigned weights on their relative importance. Let’s assume that the proximity to cycling routes is not as important as the other two factors in the final decision. We can give more important to other criteria as below.

• Theft Density: Weight=2
• Routes: Weight=1
• Existing Parking: Weight=2

We can use the the similar expression, but we must take care to normalize the output by dividing the total with the sum of all weights.

gdal_calc.py \
  -A thefts_reclass.tif \
  -B routes_reclass_resampled.tif \
  -C parking_reclass_resampled.tif --outfile suitability.tif \
  --calc="(2*A + B + 2*C)/5" --NoDataValue=0

The resulting raster represents the suitability of new bicycle parking facilities according to the chosen criteria and their relative weights.

You can also run this analysis as a batch process by putting the commands in a shell script and running it. Below is a complete script to run the weighted multi-criteria overlay analysis.

gdal_rasterize -burn 1 -add -tr 100 100 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l bicycle_thefts thefts.tif -ot Int16
gdal_rasterize -ot Int16 -burn 1 -tr 10 10 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l cycling_routes routes.tif
gdal_rasterize -ot Int16 -burn 1 -tr 10 10 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l cycle_parking existing_parking.tif
gdal_proximity.py routes.tif routes_proximity.tif -ot Int16 -distunits GEO
gdal_proximity.py existing_parking.tif existing_parking_proximity.tif -ot Int16 -distunits GEO
gdal_calc.py -A thefts.tif --outfile thefts_reclass.tif --calc="100*(A>20) + 50*(A>10)*(A<=20) + 10*(A<10)"
gdal_calc.py -A routes_proximity.tif --outfile routes_reclass.tif --calc="100*(A<=50) + 50*(A>50)*(A<=100) + 10*(A>100)"
gdal_calc.py -A existing_parking_proximity.tif --outfile existing_parking_reclass.tif --calc="100*(A>100) + 50*(A>50)*(A<=100) + 10*(A<50)"
gdalwarp -r average -tr 100 100 existing_parking_reclass.tif parking_reclass_resampled.tif
gdalwarp -r average -tr 100 100 routes_reclass.tif routes_reclass_resampled.tif
gdal_calc.py -A thefts_reclass.tif -B routes_reclass_resampled.tif -C parking_reclass_resampled.tif --outfile suitability.tif --calc="(2*A + B + 2*C)/5" --NoDataValue=0

Interested in learning GDAL Tools? Check out my free course Mastering GDAL Tools.

Data Credits

The datasets used in this post were sourced from the following open-data portals.

• Bicycle Parking: Cycling Infrastructure data: Public Tfl data. Contains OS data © Crown copyright and database rights 2016 and Geomni UK Map data © and database rights [2019]
• Cycling Thefts: London Individual crime and anti-social behavior (ASB) incidents, including street-level location information and subsequent police and court outcomes associated with the crime. Published by Single Online Home National Digital Team under Open Government Licence v3.0. Downloaded from https://data.police.uk/
• Westminster Boundary: Statistical GIS Boundary Files for London. Contains National Statistics data © Crown copyright and database right [2015], Contains Ordnance Survey data © Crown copyright and database right [2015]