Creating IDW and Spline Interpolation Maps Using QGIS

From CUOSGwiki
Revision as of 18:55, 14 December 2014 by TimothyChadwick (talk | contribs)
Jump to navigationJump to search


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, symbolizing the interpolated raster surfaces, clipping raster surfaces with digitized polygons, and exporting the raster surfaces as a finished 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 finished 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 system (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 usually more 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 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.

Data

The dataset 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 is comprised of measured concentrations of major chemical elements (in PPM) and trace elements (PPM and PPB). The purpose of collecting this data was 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, precipitation measurements, etc). Some good sources of data that can be used 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.

QGISBanner.png







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 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.
Image1.png
















  • 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.

Image2.png


  • 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.
Image3.png















  • 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 geographic coordinate system such as WGS-84 or NAD-83. Then click ‘ok’ to load the xy point data.
Image4.png









  • 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.
Image5.png









  • 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 UTM projected coordinate system 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 (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’.
Image6.png









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’.


Image7.png







  • 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’.
Image8.png









  • 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’.
Image9.png









  • 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.
Image10.png









  • 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’.
Image11.png












  • 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’.
Image12.png








  • 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.
Image13.png









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’.
Image14.png





  • 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.
Image15.png






  • 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’.
Image16.png














  • 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.
Image17.png









  • 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.
Image18.png









  • 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’.
Image19.png












  • 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.
Image20.png






  • 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.
Image21.png









  • 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’.
Image22.png









  • 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).
Image23.png











  • 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.
Image24.png









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’.
Image25.png






  • 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.
Image26.png













  • 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.
Image27.png












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








  • 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.
Image29.png










  • 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.
Image30.png
















  • 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.
Image31.png









  • 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.
Image32.png










  • 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.
Image33.png










  • 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.
Image34.png











Creating the Spline Interpolated Raster Surfaces

  • The next interpolated raster surface that will be created in this tutorial is a thin plate spline (hyperlink) 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.

Image35.png


  • 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 get achieve smoother value 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. Assign an output filename of Spline_1200m and run the tool.
Image36.png













  • 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 (use local hyperlink). For the output filename, use Spline_1200m_Reclass. Then click run and you should get the following output raster image :
Image37.png











  • 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 (insert local hyperlink). 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.
Image38.png










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










Exporting the Final Maps

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

Image41.png


  • 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 graphic editing tools within this print composer, so feel free to explore and play with them before exporting the final map.


Image40.png











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









  • There are a number of formats available to export the final map as, which including: image 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 (hyperlink) for final touch-ups if needed. To export the map, click on the ‘composer’ tab and select the ‘Export as SVG’ tab.

Image44.png


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

Conclusion

In conclusion, this tutorial showed 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 an IDW and Thin Plate Spline interpolation surface. 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.

Resources

Retrieved from "https://dges.carleton.ca/CUOSGwiki/index.php?title=Creating_IDW_and_Spline_Interpolation_Maps_Using_QGIS&oldid=5338"