Creating IDW and Spline Interpolation Maps Using QGIS

Learning Outcomes of this Tutorial

By the end of this tutorial, you will be able to perform 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 Inverse Distance Weight) and spline interpolation methods
• symbolizing the interpolated raster surfaces
• clipping raster surfaces with digitized polygons
• exporting the raster surfaces as a finished map

Purpose

The purpose of this tutorial is to introduce users to IDW and Spline interpolations, using several tools and plugins within QGIS, with the goal of creating interpolation maps.

IDW uses the measured values surrounding the prediction location to predict a value for any unsampled location, based on the assumption that things that are close to one another are more alike than those that are farther apart. For more information on IDW visit How IDW works.

Spline Interpolates a raster surface from points using a two-dimensional minimum curvature spline technique. The resulting smooth surface passes exactly through the input points.

IDW and Spline interpolations are similar in that they are methods that create surfaces based on similarity of the data or the degree of smoothness. However, spline differs from IDW in that it passes through each sample point. IDW does not pass through any points.

Introduction

QGIS is a free and open source geographic information system (GIS) that is developed and supported by thousands of users and organizations around the world. It's core features include a suite of vector and raster tools that are compatible with many of the most popular geospatial file formats

This tutorial was completed as a partial requirement of the GEOM-4008 GIS course at Carleton University as the final course project. For a complete list of other tutorials created for this course by students in the past, visit the CUOSG page.

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 step-by-step instructions and screenshots for each step of the procedures which should be easy to follow along at any skill level.

Open Source

The term ‘open source’ refers to accessible, free and redistributable programs available to the public. Compared to more common commercial packages such as ArcGIS, which have costly license subscriptions, open source programs like QGIS are free of charge to use. Additionally, open source programs and data allows individuals to work together to code. This type of model encourages open collaboration. If you would like, feel free to look at available open source technology at OSGeo.

About QGIS

As described on their website, QGIS is a "A Free and Open Source Geographic Information System", and the latest version of the software (v3.26 Buenos Aires as of September 2022) is available for download on their Official Website. Similar to 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.

Installing QGIS (v.2.2)

There are several versions of QGIS available for download. For the purposes of this tutorial, we will be using the long-term release version (QGIS 3.22).

• Follow this link QGIS Download to the QGIS site.
• Ensure that you are on the “Installation Download” tab.
  • There are three versions of the QGIS installer available in this tab. The first is a network installer which is not needed for our application ,the other two are the latest release (Version 3.26) and the long-term release (Version 3.22). We will be using the long-term release because it provides the most stability and is more widely used due to it receiving minor updates less frequently.
• The slightly older version of QGIS that will be used for the purpose of this tutorial is v.3.22 which can be downloaded from this link.
• Once you have downloaded the software from the QGIS website, follow the step-by-step installation instructions on your computer. After the QGIS v.3.22 LTR is downloaded and ready to execute, double click the QGIS shortcut icon on your Desktop or search for it in the Start Menu. The QGIS window opens, and we are ready to begin.

Example 1 IDWs of Chemical Concentrations

Data

The dataset selected for the purpose of this tutorial is publicly available from the Government of British Columbia (GoBC) from the Regional Geochemical Survey. The joint federal-provincial Regional Geochemical Surveys (RGS) have been carried out in British Columbia since 1976 as part of the National Geochemical Reconnaissance (NGR) program to aid exploration and development of mineral resources. The British Columbia Geological Survey (BCGS) maintains provincial geochemical databases capturing information from multi-media surveys. The latest release (as of September 2022) augments the database with new RGS data compiled from BCGS and Geoscience BC publications between 2016 and 2019. For more information visit the Regional Geochemical Survey Page.

This tutorial uses data provided by the GoBC. However, any point dataset can be used for the purpose of this tutorial, as long as the points have a geographic location and an attribute that can be used for the interpolation (elevation, chemical concentration, precipitation measurements, etc).

Downloading the Data

1. Go to the Regional Geochemical Survey Data Page, scroll down to the Data and Resources section and click download to the right of the 2020 database Excel file labeled RGS2020 data (see figure 2).
2. It is recommended that you also download the metadata file associated with the data if you would like to know an explanation of the various headings.
3. We will also be using a boundary shapefile to create the extent of the survey data’s raster from Administrative Boundaries in Canada, Download the Prepackaged Shapefiles then add the ‘geo_political_region_2’ files to your project folder.

Editing the Data

This data has very large amounts of variables. As such we need to edit it to make it easier to work with:

1. Open the Excel file (or other spreadsheet software- if opening a CSV file make sure that you select "separated by a comma" to format the file correctly).
2. Beginning with the variable to include in the set, in this case we are using arsenic concentration in parts per billion (ppb) in filtered, acidified water by ICP-MS.(As_ppb). Use the metadata file to locate the chemical concentration of interest and delete all other variables not dealing with location or identification (see figure 3).

Fig x. Delete sheets and variables that are not of interest

1. Once you have done this save your data as a CSV file (File, export, change file type, CSV, Save As) and you are ready to bring them into QGIS.

Example 2 Meteorological Data

Loading Shapefiles 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 the appropriate geographic coordinate system for your dataset, such as WGS-84 or NAD-83. 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 the ‘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 assign the filename 'FrameLake_Data' for the output shapefile. You will also need to specify a projected coordinate system for the output shapefile, in order to run interpolations on it in later steps. For North American datasets, one would commonly use the NAD83 UTM projections, making 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 (you can browse for NAD83(NSRS2007) / UTM zone 11N using the CRS 'Browse' button). And lastly, check the ‘Add saved file to map’ box. Then click ‘ok’.

Adding Basemaps and Labeling Data 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 (going to be used to create 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 then appear that shows what plugins have already been 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 sample stations are now distributed around the areal 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’ check-box and then select ‘AsPPM’ from the drop-down box of attribute data. In this window you can specify the font, style, size, color, etc. of the labels. I changed the font colour to white in this tutorial to contrast the lake colour, 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 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 drop-down 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’ check-box. 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 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 terms of data ranges. 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’ drop-down box. Then in the ‘Generate new color map’ area, choose the ‘spectral’ colour ramp, check the ‘invert’ check-box, 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.

Creating Clipping Polygon and Clipping Interpolation Surface

• This portion of the tutorial will provide instructions on how to create a clipping mask for the interpolation raster surface. The idea is to digitize a polygon of the lake boundary and then clip the interpolated surfaces with it, so that the interpolation does not extend onto land (which may have very different As concentrations compared to water). In order to do so, a new shapefile needs to be created. Under the ‘Layer’ tab, click on ‘New’ and then select ‘New Shapefile Layer’.

• In the ‘New Vector Layer’ window, select the option to create a ‘polygon’ in the ‘type’ field. Then specify the correct projection for the shapefile (for the purpose of this tutorial use NAD83 / UTM zone 11N). Seeing that this shapefile is just going to be used to create a clipping mask polygon, no ‘new attributes’ need to be assigned. Click 'ok'. A window will then appear which allows you to specify a filename and location to save the new shapefile. Call is ‘Clip_Polygon’, and save it.

• You will now notice your new polygon shapefile in the list of layers. The next step will be digitizing a polygon of Frame Lake. Right-click on the ‘Clip_Polygon’ layer and select the ‘Toggle editing’ button.

• Then click the ‘Add Feature’ button to start creating the polygon.

• Start digitizing the polygon by clicking on the edge of the lake and following the shoreline all around the lake. Use as many digitizing points as possible to ensure an accurate polygon.

• Once you have finished digitizing your polygon, right click on the last point. A window will appear that will allow you to specify an id for the polygon, assign it a value of ‘1’. Then click 'ok'. Now that we have our Lake polygon, we can save the edits made to the shapefile and close the editor. Click on the ‘current edits’ button and select the ‘save for all layers’ button. Then right-click on the Clip_Polygon shapefile layer and click on the ‘Toggle editing’ button to stop editing the polygon.

• Now we are ready to clip the ‘Reclassified grid’ using the ‘Clipping_Polygon’. To do so, search for the SAGA tool ‘Clip grid with polygon’ in the processing toolbox, and open it. In the ‘clip grid with polygon’ tool window, assign the ‘Reclassified Grid’ as the input raster layer and the ‘Clip_Polygon’ as the input polygon. Specify an output filename ‘IDW_500m_Clipped’ and click run.

• You will notice that the reclassified grid is now clipped, but the default original greyscale colouring has returned. There is a quick way to re-symbolize the clipped raster with the appropriate colour ramp.

• Right-click on the ‘reclassified grid’, which still has its intact colour ramp and open the properties menu. Click on the ‘export color map file’ button, and save the text file to hard drive. Then re-open the properties menu for the ‘IDW_500m_Clipped’ raster, select the ‘style’ tab, and then pick ‘Singleband Pseudocolor for the ‘render type’. Then click on the ‘load colour map from file button’ (to the left of the ‘export color map to file’ button) and open the saved text file from the earlier step. You will now notice that your legend is back to the original format.

• And now the map is complete and ready to be checked for accuracy. When you inspect the As concentration point values and the interpolated surface value under each point (looking at the legend), you will see that all the points fall in the correct interpolated value range area, indicating a good interpolation job. Now that the IDW interpolation map is complete, we will now repeat the process using the 'thin plate spline' interpolation tool.

Creating the Spline Interpolated Raster Surfaces

• The next interpolated raster surface that will be created in this tutorial is a thin plate spline interpolation surface, using the Frame_Lake point dataset as the input. Use the same map window as the IDW interpolation, just turn off the ‘IDW_500m_clipped’ raster layer so that just the data points are visible again. In the processing toolbox search for ‘Thin Plate Spline (local)’ which will result in one SAGA spline tool with that title. Open the SAGA ‘Thin plate spline (local)’ tool.

• Use the ‘FrameLake_Data’ for the input data points. Again, use AsPPM for the input attribute. Use the default value of 0.000100 for the ‘Regularization’ input. For the search radius, more sample locations needs to be used for this thin plate spline interpolation compared to the previously done IDW interpolation, in order to obtain smoother interpolated boundaries. Therefore, a search radius of 1200m was assigned for this tutorial. The remainder of the inputs were left with their default value, except cell-size which was changed to 1.0. For more information about how splines work, visit the following website. Assign an output filename of Spline_1200m and run the tool.

• Rename the output grid layer from ‘Grid’ to ‘Spline_1200m’. The next step will be to reclassify this output grid with the same range of As concentration value used in the IDW interpolation tutorial. Open the ‘Reclassify grid values’ tool and follow the instructions from the IDW tutorial to reclassify the raster grid Reclassifying tutorial. For the output filename, use Spline_1200m_Reclass. Then click run and you should get the following output raster image :

• Next we will clip the ‘Spline_1200m_Reclass’ raster with the lake outline shapefile, using the same technique that was used for the IDW interpolation surface raster Clipping tutorial. Use the ‘Spline_1200m_Reclass’ as the input, and specify an output filename of ‘Spline_1200m_Clip’. Now that the spline interpolation surface has been clipped, the last step is to assign a colour ramp to it. Use the same procedure and range values that were used for the IDW interpolation tutorial Symbolizing legend tutorial.

• And that is all for the Spline portion of the tutorial. The resulting final spline interpolation surface should look like this:

• And now the map is complete and ready to be checked for accuracy. When you inspect the As concentration point values and the interpolated surface value under each point (looking at the legend), you will see that all the points fall in the correct interpolated value range area, indicating a consistent interpolation job. The final portion of this tutorial will deal with exporting finished maps from QGIS.

Exporting the Final Maps

• Under the ‘Project’ tab, select “New Print Composer’.

• Assign a ‘print composer tittle’ that will be used for the new print composer (For example: ‘Spline Map’). Under the ‘Layout’ tab, you can add a number of map elements to the print composer including the actual map, a legend, scale bars, etc., depending on the nature and purpose of the final map. For this tutorial, experiment by adding the spline interpolation map, map title (under ‘add label’), a legend and scale bar. A few tips: Make sure you zoom to the layer extent of the spline interpolation map first before adding it as a map to the print composer, text editing for the ‘label’ is done in a text editing box on the right hand side of the print composer window. There is a wide variety of graphics editing tools within this print composer, so feel free to explore and play with them before exporting the final map.

• An example final output map is shown below. The last step of this tutorial will be exporting the map.

• There are a number of formats available to export the final map as, which including: JPEG file, PDF or as an SVG file. For the purpose of this tutorial, we will export the map as an SVG file, so that it can be loaded into an open source graphics editor such as Inkscape for final touch-ups if needed. To export the map, click on the ‘composer’ tab and select the ‘Export as SVG’ tab.

• Assign a file name to the exported SVG file, and click save. And that concludes this tutorial.

Conclusion

In conclusion, this tutorial aimed to show the user how to import point data into QGIS, run interpolations on the point data, clip and symbolize the data and export the interpolated data as a finished map. I hope this tutorial provided adequate instructions on how to download the software and data, and to perform the various tasks needed to create IDW and Thin Plate Spline interpolation maps. Unfortunately, I did not have the time to create a tutorial on kriging using an open source GIS program, but this is certainly still something I would like to do in the future. Please feel free to comment on anything in this tutorial that is poorly explained or does not make any sense, and feel free to contact me by email if you have further questions.

Resources

• ArcGIS Resources (August 11, 2012). ArcGIS Help 10.1:Understanding interpolation analysis. http://resources.arcgis.com/en/help/main/10.1/index.html#//009z0000006w000000

• ArcGIS Resources (August 11, 2012). ArcGIS Help 10.1:How IDW works. Retrieved on December 14, 2014 from http://resources.arcgis.com/en/help/main/10.1/index.html#//009z00000075000000

• ArcGIS Resources (August 11, 2012). ArcGIS Help 10.1:How spline works. Retrieved on December 14, 2014 from http://resources.arcgis.com/en/help/main/10.1/index.html#//009z00000078000000

• GeoBase (December 08, 2014). Geobase database homepage. Retrieved on December 14, 2014 from http://www.geobase.ca/geobase/en/

• GeoGratis (2014). Natural Resources Canada: GeoGratis database homepage. Retrieved on December 14, 2014 from http://www.geogratis.gc.ca/geogratis/Home?lang=en

• GIS Data Depot (2014). Geo-Community: The GIS Data Depot. Retrieved on December 14, 2014 from http://data.geocomm.com/

• Free GIS Datasets (December 11, 2014). Free GIS Datasets homepage. Retrieved on December 14, 2014 from http://freegisdata.rtwilson.com/

• QGIS. (2014). QGIS: A free and Open Source Geographic Information System. Retrieved December 14, 2014 from http://www.qgis.org/en/site/index.html

• QGIS. (2014) Download QGIS for your Platform. Retrieved on December 14, 2014 from http://www.qgis.org/en/site/forusers/download.html

• QGIS. (2014). Features of QGIS. Retrieved on December 14, 2014 from http://www.qgis.org/en/site/about/features.html