# Weighted Multi-Criteria Overlay Analysis using GDAL Tools

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

# Fixing Rasters with Missing Data using QGIS, GDAL and Python

When working with raster data, you may sometimes need to deal with data gaps. These could be the result of sensor malfunction, processing errors or data corruption. Below is an example of data gap (i.e. no data values) in aerial imagery.

# Reclassifying Rasters using gdal_calc

gdal_calc is one of the lesser used tools among the GDAL utilities. There aren’t many examples of using it in the wild and some advanced features are not well documented. I am finding myself using it a lot lately and have discovered some really powerful use cases.

# Convert between Orthometric and Ellipsoidal Elevations using GDAL

When working with elevation data, sometimes you may discover that 2 datasets from different providers have very different elevation values for the same location. A common reason for this being each dataset being referenced to a different surface.

# Creating Geospatial PDFs with GDAL Tools

GeoPDF is a unique data format that brings the portability of PDF to geospatial data. A GeoPDF document can present raster and vector data and preserve the georeference information. This can be a useful format for non-GIS folks to consume GIS data without needing GIS-software. While GeoPDF is a proprietary format, we have a close alternative in the open Geospatial PDF format. GDAL has added support for creating Geospatial PDF documents from version 1.10 onwards. In this post, I will show how to create a GeoPDF document containing multiple vector layers.

## Get the Tools

### Windows

OsGeo4W is the best way to install GDAL on Windows. The default installation gives your GDAL tools with PDF format support. You can use the GDAL tools via the OsGeo4W Shell included in the install.

### Mac

KyngChaos providers a convenient GDAL installer for Mac. You also need to install the additional GeoPDF plugin to enable support for PDF format.

Once installed, add the path to GDAL library to your `.bash_profile` file to be able to use the commands easily from the terminal. Launch a Terminal and type in the following commands.

```echo 'export PATH=/Library/Frameworks/GDAL.framework/Programs:\$PATH' >> ~/.bash_profile
source ~/.bash_profile
```

### Linux

Installation instructions will vary with the distribution. On Ubuntu, you can install the `gdal-bin` package.

`sudo apt-get install gdal-bin`

## Verify GDAL Install

If you already have GDAL installed, or just installed it, run the following command in a terminal to verify that your GDAL installation is working and has support for GeoPDF format.

```gdalinfo --formats | grep -i pdf
```

If you see Geospatial PDF printed in the output – you are all set. If you do not get any output or get an error, your install is not correctly configured.

## Get the Data

For this example, I chose to use OpenStreetMap Metro Extracts from MapZen. Download the shapefiles (OSM2PGSQL SHP format) for the city of your choice. I am using the extract for Bangalore city in this example. Unzip the downloaded file to a folder on your computer.

## Procedure

The process for creating a GeoPDF file from a bunch of shapefiles is the matter of running a single `gdal_translate`command. But we need to prepare the data and figure out the correct command-line options. So follow along to understand how you can arrive at the final command – or simply scroll to the end to see the final command-line.

latuviitta.org has a comprehensive overview of all the options available for GeoPDF creation via GDAL. The follow steps are adapted and simplified version of that guide.

• First step is to create a `.vrt` file that can hold all the vector layers we want in the PDF. If you just need a single layer in the PDF, you can skip creating the `.vrt` file and directly reference the layer in place of the VRT. Note the <SrcSQL> tag in the VRT file. This is for filtering out all features where the ‘name’ field is empty. You can leave that out or modify to suit your dataset. Name this file `osm.vrt` and save it on the same folder with your data.
``````<OGRVRTDataSource>
<SrcDataSource>bengaluru_india_osm_line.shp</SrcDataSource>
<SrcLayer>bengaluru_india_osm_line</SrcLayer>
<SrcSQL dialect="sqlite">SELECT name, highway, geometry from bengaluru_india_osm_line where name is not NULL</SrcSQL>
<GeometryType>wkbLineString</GeometryType>
<LayerSRS>WGS84</LayerSRS>
</OGRVRTLayer>
<OGRVRTLayer name="pois">
<SrcDataSource>bengaluru_india_osm_point.shp</SrcDataSource>
<SrcLayer>bengaluru_india_osm_point</SrcLayer>
<SrcSQL dialect="sqlite">SELECT name, geometry from bengaluru_india_osm_point where name is not NULL</SrcSQL>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>WGS84</LayerSRS>
</OGRVRTLayer>
</OGRVRTDataSource>``````
• GeoPDF is a raster format that can overlay vectors on top. So we need a raster layer as the base. If you have some satellite imagery or scanned raster for the area, you can use it as the base layer, or we can create an empty raster for the extent of the vector layer. `ogrtindex` command creates a bounding box polygon from the given input layers. `gdal_rasterize` command then fills this polygon with the given value and creates a raster. The `-tr` option specifies the pixel resolution of the raster in degrees. You can tweak that to get the output size you need. `cd` to the directory where you have extracted the vector layers and run the following commands.
```cd Users\Ujaval\Downloads\bengaluru_india.osm2pgsql-shapefiles

ogrtindex -accept _different_schemas extent.shp osm.vrt

gdal_rasterize -burn 255 -ot Byte -tr 0.0001 0.0001 extent.shp bangalore.tif```
• Now we can convert the empty `bangalore.tif` raster to a PDF – overlaying the vector layers from the `osm.vrt` file.
`gdal_translate -of PDF -a_srs EPSG:4326 bangalore.tif bangalore.pdf -co OGR_DATASOURCE=osm.vrt -co OGR_DISPLAY_FIELD="name"`
• Once the conversion finishes, you can open the resulting `bangalore.pdf` file in any PDF viewer. Opening it in Adobe Acrobat viewer, you can see the map data layers. You can browse the features in the layer panel, search for any attribute value and zoom/pan the map.
• Another popular use of GeoPDF files is to use it as offline base maps using programs such as Avenza PDF Maps. Loading the `bangalore.pdf` file on Avenza Maps on your mobile phone, you can use the GPS to view your current location or trace a GPS route on top. Search also works across layers in the PDF.

You can download the sample bangalore.pdf Geospatial PDF format file for exploring the format yourself.

# Spatial Joins on the Command Line

GDAL and OGR libraries come with handy command-line tools. These tools are quite powerful and can save you a lot of effort if you know how to use them. Here I will show how to use the ogrinfo and ogr2ogr tools to perform spatial joins. A single command can do complex operations on your spatial data and save you a lot of clicking-around and data-munging in a GIS.

## Get the Tools

The best way to get the command-line tools on Windows is via the OSGeo4W Installer. If you are on Linux or Mac, see these instructions to get the package for your platform.

## Get the Data

Review the data and problem statement from the Performing Spatial Joins tutorial. Download the Borough Boundaries and Nursing Homes shapefiles.

## Procedure

OGR command line tools accept only 1 input. But we have 2 inputs for the spatial join. An easy way to fix this, is to use a VRT file. A VRT file allows us to specify multiple inputs and pass them to the command-line tool as layers of a single input.

Unzip the input shapefiles in a single folder on your drive. Create a file named input.vrt in the same folder with the following content.

```<OGRVRTDataSource>
<OGRVRTLayer name="boroughs">
<SrcDataSource>nybb.shp</SrcDataSource>
<SrcLayer>nybb</SrcLayer>
</OGRVRTLayer>
<OGRVRTLayer name="nursinghomes">
<SrcDataSource>OEM_NursingHomes_001.shp</SrcDataSource>
<SrcLayer>OEM_NursingHomes_001</SrcLayer>
</OGRVRTLayer>
</OGRVRTDataSource>```

Open the OSGeo4W shell and cd to the directory containing the shapefiles and the vrt file. Run the ogrinfo command to check if the VRT file is correct.

`ogrinfo input.vrt `

OGR tools can run SQL queries on the input layers. We will use the ST_INTERSECTS function to find all nursing homes that intersect the boundary of a borough and use the SUM function to find the total nursing home capacity of a borough. Run the following command.

```ogrinfo -sql "SELECT b.BoroName, sum(n.Capacity) as total_capacity from
boroughs b, nursinghomes n WHERE ST_INTERSECTS(b.geometry, n.geometry) group
by b.BoroName" -dialect SQLITE input.vrt```

You can see that in a single command we got the results by doing a spatial join that takes a lot of clicking around in a GIS environment. We can do a reverse spatial join as well. We can join the name of the Borough to each feature of the Nursing Homes layer. Using the ogr2ogr tool we can write out a shapefile from the resulting join. Note that we are adding a geometry column in the SELECT statement which results in a spatial output. Run the following command:

```ogr2ogr -sql "SELECT n.Name, n.Capacity, n.geometry, b.BoroName from
boroughs b, nursinghomes n WHERE ST_INTERSECTS(b.geometry, n.geometry)"
-dialect SQLITE output.shp input.vrt```

Open the output.shp in a GIS to verify that the new shapefile as attributes joined from the intersecting borough. You can use ogrinfo command to check that as well.

`ogrinfo -al output.shp `