What to expect from this Tutorial

Going through this tutorial using ' Quantum' GIS ([1]): You will be gain the following skills.

Importing CSV file
creating interpolated raster surfaces using IDW
symbolizing the interpolated raster surfaces
clipping raster surfaces with base map polygons
exporting the raster surfaces as a finished map.

Purpose

The purpose of this tutorial is to expand your ability to run important analysis tools such as Interpolation and Image analysis (IDW) [[2]], using several proprietary tools [3], with the end goal of creating a finished IDW interpolation map. The difficulty level for this tutorial is an entry-level QGIS user, But it is assumed that the user has some experience with other GIS systems Such as ArcPro or ArcMap. I chose to create this tutorial to try and bridge the gap between this software and Esri’s products. Making it easier to learn how to use the software and carry out the functions that the user is already familiar with.

This tutorial will be incorporating some screen-shots for each step to show the path I took to create all the layers that are in the finished product

Introduction

The purpose of this tutorial is to learn how to use the proprietary tools in an open-source platform. I decided to do this in QGIS as the software is often looked at as the main competitor to Esri as the design and application are similar to that of ArcPro and ArcMap. This tutorial looks to try and bridge the gap between this software and Esri’s products, to make it easier to learn how to use the software and carry out the functions that the user is already familiar with. Although the software looks similar differences between QGIS and Esri’s products which are often taken as the industry standard. This all ties into making this tutorial practical and accessible for all users picking up this software for the first time, or just wanting to expand their knowledge. This was chosen to help show off the strength of open-source software in a practical way encouraging those to use it more often.

open-source

open-source ([4]) is a term that refers to having accessible, free and redistributable programs available to the public. This allows for individuals to work on the code together in order to create a model that encourages open collaboration and develops code together.

GIS

In our case, a license subscription to Esri can be extremely expensive to have just one account leading to a lot of overhead costs for full teams. This is why having open-source and freely available software is vital for users to know to bring down costs where possible. To look at GIS programs available for free download, you can visit the [5] page to look at all the projects currently being worked on by the community.

Downloading QGIS (3.22)

for this tutorial, you will need to have QGIS installed [[6]] on your desktop or laptop. at the time of creating this tutorial, QGIS 3.22 is the most recent stable version and the version all of this tutorial is being done in. Once you have selected the version you wish to download based on your operating system, follow the guide that is installed in the package. With QGIS 3.22 now installed on your computer and a proper path for accessing the application on your hard drive, you are now able to start the tutorial!

This tutorial was created as a result of the final project given out in the GEOM-4008 GIS course at Carleton University. If you wish to look at other tutorials created by other students, this can be found at [[7]]

Data

The data I selected for this tutorial goes along with the theme of keeping everything open source and using publicly available data from the Canadian government. To ensure the longevity of this data we will be using monthly climate summaries from Stat CAN as it is very extensive and easily accessible data for all provinces and territories in Canada. This in-depth look at provincial weather will be the best and easily malleable data set to use for our project. For this case, we will be looking at Alberta's average mean temperature. This will give us a good well-distributed data set that we will be able to use to create an IDW model and help us guess the values of areas that we do not know.

Downloding the data

Go to [[8]] and scroll to the bottom of the page
being by clicking the year 2 years before your current date; (as of today it's Decmeber 2021 so I would use December 2019 etc)
Select your area of interest (for this assessment I used Alberta as it has the most climate stations in the country)
ensure that you are downloading the CVS data and click the download button when ready.
Be sure to also download the Legend(txt) file as it will allow you to see what all the data is meant to represent
We now need our boundary file, to get this go to [[9]] and download the ArcShape file and click the provinces under "Cartographic Boundary File"
this will bring you to a second screen that says "Archived – gpr_000b11a_e.zip (ZIP version, 40,968.0 kb)" click this and download this for the project

Editing the data

With this data having so many variables we need to edit it to make it easier to work with:

Open your CVS file using Excel and make sure that you select "separated by a comma" to format the file correctly.
beginning with the variable to include in the set, use your legend(txt) file to locate the mean temperature and delete all other variables not dealing with location or identification
right click on the "Clim_ID" column to select the entire row
under the home tab of Exel click the sort/filter drop-down menu and select largest to smallest.
You should see a bunch of stations that are all formatted to the left of the cell and will contain letters unlike the rest of the data set. Select all these rows and delete them from your file (do not worry, they are not needed and can create some issues when loading it into QGIS)
Once you have done this save your data and you are ready to bring it into QGIS
1. note: If you are having issues with bringing the data into QGIS make sure that all of your climate ID's are correctly formatted as this can stop the data from being able to be read by QGSI

Tutorial

Loading CSV File and Provincial Boundary File

Launch QGIS (version 3.22 or newer) and ensure that you have downloaded the correct data sets. and create a new project file

at the top of your toolbar, click the plugin tab and click the "manage and install plugins" item.

In the "all" tab, search for the "QuickMapServices" plugin and install it to QGIS. once this is done restart QGIS

once you are back into a new project layer, click the web tab at the top and go to the QuickMapServices dropdown and select oms standard. (this should add in a base map to your project)

we will now add in our 2 layers starting with the weather stations. by clicking on the layer tab at the top of the page go to the "add layer" tab and scroll to the "Delimited Text layer"

in the new window, click the 3 dots at the top right corner of the screen (you will be finding the CVS file that we edited earlier) once you have found it make sure that the "Comma" and "tab" boxes are checked off in the "custom delimiters" file format

in the Geometry definition option set your x field to "Long" and your y field to "Lat" so that QGIS knows how to read your coordinates

click on the layer tab at the top of the page go to the "add layer" tab and scroll to the "add Vector layer" tab

in the new window, click the 3 dots at the top right corner of the screen (you will be finding the boundary file we downloaded) in the package you downloaded click the file with the .shp end to it and add it to your project
1. once it is loaded in right-click on it and select "Open attribute table"
2. once open click the pencil at the top of the window. then select all layers besides Alberta and click the " delete selected feature" tap on the toolbar you clicked the pencil on
3. once this is complete correctly you should only have the Alberta boundary file on your project layer (and the weather stations)

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.

Creating a IDW model in QGIS

Contents