Creating IDW and Spline Interpolation Maps Using QGIS

Performing IDW and Spline Interpolation Methods Using QGIS

Learning Outcomes of this Tutorial

Performing the following tasks within 'Quantum' GIS (QGIS): Importing point data files (shapefile, CSV); importing satellite imagery basemaps; labeling point data using attribute information; creating interpolated raster surfaces using IDW and spline interpolation methods (point shapefile being used for the model input), symbolizing interpolation raster surfaces, clipping raster surfaces with digitized polygons, and exporting the raster surfaces as a map.

Purpose

The purpose of this tutorial is to show users how to perform Inverse Distance Weight (IDW) and Spline interpolations, using several tools and plugins within QGIS, with the goal of creating interpolation maps. For the knowledge level of this tutorial, it is assumed that the user has some experience with GIS systems, whether it be from commercial products such as ESRI ArcMap or MapInfo or from open source packages such as QGIS, GRASS GIS, etc. However, this tutorial incorporates screen-shots for each step of the procedures and should be easy to follow along to.

Introduction

The purpose of this project was to explore toolsets and applications within different ‘open source’ geographic information systems (GIS) programs. The term ‘open source’ refers to the fact that all of these programs are free of charge to use, compared to more commonly known/taught (and user friendly) commercial packages such as ArcGis, which have costly license subscriptions. For a good listing of open course GIS programs available for download, visit the Open Source GIS page. This project was completed as a partial requirement (final course project) of the GEOM-4008 GIS course at Carleton University. For a complete list of other tutorials created for this course by students in the past, visit the CUOSG page.

Data

The data set being used for the purpose of this tutorial is a lake water geochemical dataset, consisting of 25 sampled locations on Frame Lake (Yellowknife, North West Territories), provided by the Patterson Research Group at Carleton University. This dataset has measured concentrations of major chemical elements (in PPM) and trace elements (PPM and PPB). The purpose of collecting this data is to map the concentration distributions of harmful elements, such as As, Pb, Hg and Cr within the lake, in order to determine which parts of the lake need the most remediation work done before trying to reintroduce fish populations into the lake.

However, this is unpublished scientific research data and therefore cannot be shared yet. However, any point data-set can be used for the purpose of this tutorial, so long as the points have a geographic location and an attribute that can be used for the interpolation (elevation, chemical concentration, etc). Some good sources of data for this tutorial include: Geobase, Geogratis, GIS Data Depot and Free GIS Data .

About QGIS

As described on their website, QGIS is a "A Free and Open Source Geographic Information System", and the latest version (v2.6) of the software is available for download on their Official Website. Like other GIS platforms it can be used to edit, create, visually represent, analyze and export a variety of geospatial data and it is licensed under the GNU General Public License. QGIS is an official project of the Open Source Geospatial Foundation (OSGeo). It runs on Linux, Unix, Mac OSX, Windows and Android and supports numerous vector, raster, and database formats and functionalities.

Acquiring QGIS (v.2.2) for this Tutorial

The slightly older version of QGIS that will be used for the purpose of this tutorial is v.2.2 which can be downloaded from this link link. Once you have downloaded the software from the QGIS website, follow the step-by-step installation instructions on your computer. With QGIS v.2.2 now installed on your computer and a suitable point dataset to be interpolated saved on your hard drive, you are now ready to start the tutorial.

Tutorial

Loading shapefile and CSV file data

• Launch QGIS (version 2.2) once you have installed it on your computer, and downloaded a suitable dataset to work with. There are two simple ways to import the point data into QGIS, depending on the type of data that you downloaded (shapefile, CSV file, excel spreadsheet). If your data is in the form of an excel spreadsheet (where the points have associated latitudes, longitudes and attribute data), save it as a comma-separated value (CSV) file before trying to import it into QGIS.

• To import a shapefile, click on the ‘Layer’ tab in the upper left hand corner of the screen, and select the ‘Add Vector Data’ tab.

• A window will then prompt you to navigate to your saved shapefile dataset on your hard drive. Afterwards click ‘open’ and you should see your point data on the main display.

• To import a CSV file, click on the ‘Layer’ tab in the upper left hand corner of the screen, and select the ‘Add Delimited Text Layer’ tab.

• A window will then prompt you to navigate to your saved CSV file on your hard drive. Assign a layer name for the output point-data, for this tutorial ‘FrameLake_Dataset’ was assigned. Select the ‘longitude’ field to represent the x-values and the ‘latitude’ field to represent the y-values. You may be also prompted to provide a map projection for the data (sometimes this is done automatically). If so, select a generic geographic projection (create hyperlink) such as WGS84 (hyperlink). Then click ‘ok’ to load the xy point data.

• The next step when loading a CSV file is exporting the loaded point data (from the previous step) to a shapefile. To do this in QGIS, right-click on the loaded point data layer and click ‘Save As’ tab.

• A window will then prompt you to pick the format of the file being saved, choose ‘ESRI Shapefile’. Browse for a location to save the shapefile on your hard drive, and for this tutorial, FrameLake_Data was assigned for the shapefile name. You will also need to specify a UTM projected coordinate system (hyperlink) for the output shapefile, in order to run interpolations on it in later steps. For North American datasets, one can use the NAD83 UTM projections, just be sure to select the correct UTM zone for the map area of your dataset. For the Yellowknife dataset used in this tutorial, NAD83-UTM zone 11 was chosen for the projection (browsed for NAD83(NSRS2007) / UTM zone 11N). And lastly, check the ‘Add saved file to map’ box. Then click ‘ok’.

Adding basemaps and labeling points

• Once you have successfully loaded your shapefile data, there are a few steps to perform in order to help you visualize your dataset, which includes loading a base map (used for creating a clipping layer in future steps) and labeling your sample points with attribute data (concentrations of harmful elements in this tutorial).

• To load a satellite image basemap in QGIS, you need to install a plugin to do this task. Click on the ‘Plugins’ tab and select ‘Manage and Install Plugins’.

• A window will appear that will show what plugins have been already installed, and what plugins are available for download. (Note. Plugins are the open source equivalent of tools within the Arc Toolbox in ESRI ArcMap). In the ‘Search’ bar, type ‘OpenLayers’ and then select the ‘OpenLayers Plugin’ from the list of available plugins. Then click ‘Install Plugin’.

• After the plugin has installed, close the plugin installation window. Then click on the Plugins tab again, hover over the OpenLayers Plugin, and select ‘Add Bing Aerial Layer’.

• The basemap will then appear on your map area. Drag the ‘Bing Aerial’ layer to the bottom of the map layer in the ‘Layers’ window. Then right-click on the ‘FrameLake_Data’ shapefile and select ‘Zoom to layer extent’. You will now notice that the lake water samples are now distributed around the arial extent of Frame Lake on the Bing Aerial layer.

• The next step will be to label the data points with the arsenic concentration measurement at each station. To do this, right-click on the FrameLake_Data shapefile and select ‘Properties’.

• The properties window will then appear. Click on the ‘Labels’ tab, check the ‘Label this layer with’ checkbox and then select ‘AsPPM’ from the dropdown box of attribute data. In this window you can specify the font, style, size, color, etc. of the labels. I changed the colour to white (to contrast the lake) and then clicked ‘apply’ and ‘ok’.

• Now the As concentrations are displayed for each point on the map and the points can be viewed in relation to a basemap. Having this information available will be helpful for the next portion of the tutorial which deals with performing the interpolations.

Creating the IDW interpolated raster surfaces

• The first interpolation method that will be investigated in this tutorial is the IDW (hyperlink) interpolation. Click on the ‘Processing’ tab and select ‘Toolbox’.

• The ‘Processing Toolbox’ panel will then appear on the right side of the map view. You will notice that QGIS draws tools from 8 different GIS sources (GDAL/OGR, GRASS, SAGA, etc.) including its own QGIS tool selection. The easiest way to find an appropriate IDW interpolation tool is to type ‘inverse distance’ into the upper search column. You will notice two IDW tools to choose from, one by GRASS GIS and the other by SAGA GIS. Having looked at both, I would recommend using the SAGA IDW tool, as it is more user friendly. Double click on the SAGA ‘Inverse distance weighted’ tool.

• The SAGA IDW inputs window will then appear, which allows you to input the data layer being interpolated and customize the interpolation based on your needs. For the ‘Points’ input use the ‘FrameLake_Data’ shapefile layer. For the ‘Attribute’, select AsPPM from the dropdown selection. Use the default values for ‘Target Grid’, ‘Distance Weighting’, ‘Inverse Distance Power’ and ‘Exponential and Gaussian Weighting Bandwidth’, as they are appropriate for the data distribution in this study. For more information on these parameters, consult the following help page. For the ‘Search Range’, use the ‘search radius (local)’ option. For the ‘Search Radius’, I assigned a radius of 500 m, which ensured that every point in the dataset had at least 1 neighboring point to compare to (feel free to assign larger and smaller values to compare the effect on the output raster). ‘All directions’ was selected for the ‘Search Mode’ because the data is fairly irregularly dispersed on the map. The default value of 10 was used for the ‘Maximum number of points’ input due to the low number of points. A ‘cell size’ of 1m was selected to produce a high resolution output raster layer. And lastly, assign a location to save the output raster ‘Grid’ (I called mine IDW_500m), and check the ‘open output file after running algorithm’ checkbox. Click ‘Run’.

• The IDW interpolation tool will then output a black and white interpolated raster surface onto the map area. Drag the grid layer bellow the point shapefile layer and rename it from ‘Grid’ to ‘IDW_500m’ by right-clicking on the raster grid and selecting ‘rename’. The next step will be to reclassify the IDW_500m interpolation into 6 fields, each with a given range of values determined by the user. Unfortunately, this basic task cannot be done using the QGIS raster properties, so it will have to be done using a separate tool.

• In the processing toolbox window, search for the SAGA tool called ‘Reclassify grid values’, and then double click on it from the selection of tools.

• In the ‘reclassifying grid values’ window, select the IDW_500m grid raster. For the ‘Method’ option, select ‘simple table’ as this method allows you to input a range of values for a given field. For the ‘Lookup Table’ click on the button to the right of “Fixed Table 3x3’.

• The As concentration data for the 25 points were classified into 6 fields using a Jenks Natural Breaks Optimization (hyperlink) method in another program prior to starting this tutorial (out of the scope of this tutorial). If you are using data which is different from the data used in this tutorial, you will need to decide on how best to classify your data. In the Fixed Table window, you have the option of adding or removing as many rows as you like. For this tutorial we are going to assign 6 rows with the following values (see image) for the minimum, maximum and new values. Then click ok.

• Keep the remainder of the inputs in the ‘Reclassify grid values’ window at their default settings, and assign a filename of ‘IDW_500m_Reclass’ for the output ‘reclassified grid’. Also make sure to check the box which opens the output file after running the algorithm. Then click run. The new reclassified grid will then appear in the map window as a greyscale raster surface.

• The next step is to assign a colour ramp to this reclassified raster layer. To do this, right-click on the ‘reclassified grid’ and select ‘properties’. Select the ‘Style’ tab, and then pick ‘Singleband pseudocolor’ in the ‘render type’ dropdown box. Then in the ‘Generate new color map’ area, choose the ‘spectral’ colour ramp, check the ‘invert’ checkbox, select the ‘Equal Interval’ mode and 6 classes. Assign a min value of 140 and max value of 1340 and then click the ‘classify’ button. This will create a 6 class colour ramp (left hand box in the image) which needs to be reclassified based on the values of the ‘reclassified grid’.

• Double click on the numbers in the ‘value’ column, and replace them with the assigned values of the reclassified grid, as shown in the following image. Also update the ‘label’ field with the appropriate concentration ranges (as this is not done automatically).

• After modifying the colour ramp values to reflect the correct range of values from the reclassified grid, click ‘apply’ and ‘ok’ and you will now have coloured map of the IDW_500m interpolation surface for As concentrations in Frame Lake with a matching legend. The last step of this mapping process is to create a polygon of the lake which will be used as a mask to clip the interpolation surface.