1. The GIS Process¶

Before we start, it would be useful to give a brief overview of a process that can be used to solve any GIS problem. The way to go about it is:

1. State the Problem
2. Get the Data
3. Analyze the Problem
4. Present the Results

2. The Problem¶

Let’s start off the process by deciding on a problem to solve. For example, you are an estate agent and you are looking for a residential property in Swellendam for clients who have the following criteria:

1. It needs to be in Swellendam
2. It must be within reasonable driving distance of a school (say 1km)
3. It must be more than 100m squared in size
4. Closer than 50m to a main road
5. Closer than 500m to a restaurant

3. The Data¶

To answer these questions, we’re going to need the following data:

1. The residential properties (buildings) in the area
2. The roads in and around the town
3. The location of schools and restaurants
4. The size of buildings

All of this data is available through OSM and you should find that the dataset you have been using in lab 2 can also be used here.

If you want to download data from another area you can find more information in QGIS official manual.

4. Follow Along: Start a Project and get the Data¶

We first need to load the data to work with. This time we will load data from a GeoPackage file, which is a container that allows you to store GIS data (layers) in a single file. Unlike the ESRI Shapefile format (e.g. all the .shp datasets you used earlier), a single GeoPackage file can contain various data (both vector and raster data) in different coordinate reference systems, as well as tables without spatial information; all these features allow you to share data easily and avoid file duplication.

1. Download this GeoPackage file.

2. Create a new map.

3. Click on the Open Data Source Manager button.

4. On the left click on the GeoPackage tab.

5. Click on the New button and browse to the lab5tutorial1.gpkg file you downloaded before.

6. Select the file and press Open. The file path is now added to the Geopackage connections list, and appears in the drop-down menu.

   You are now ready to add any layer from this GeoPackage to QGIS.

7. Click on the Connect button. In the central part of the window you should now see the list of all the layers contained in the GeoPackage file.

8. Select any of the layers and click on the Add button.

   You will see that the layer has been added to the Layers panel with features displayed on the map canvas.

9. Using the same approach you can add all other layers from the file (you can also add them all simultaneously by selecting them with Shift key).

10. If you want you can add a background map. Open the Browser and load the OSM background map from the XYZ Tiles menu.

    

11. Zoom to the layer extent to see Swellendam, South Africa

Before proceeding we should filter the roads layer in order to have only some specific road types to work with.

Some of the roads in OSM dataset are listed as unclassified, tracks, path and footway. We want to exclude these from our dataset and focus on the other road types, more suitable for this exercise.

Moreover, OSM data might not be updated everywhere and we will also exclude NULL values.

1. Right click on the roads layer and choose Filter….

2. In the dialog that pops up we can filter these features with the following expression:

   "highway" NOT IN ('footway','path','unclassified','track') OR "highway" != NULL
                 

   The concatenation of the two operators NOT and IN means to exclude all the unwanted features that have these attributes in the highway field.

   != NULL combined with the OR operator is excluding roads with no values in the highway field.

   You will note the icon next to the roads layer that helps you remember that this layer has a filter activated and not all the features are available in the project.

The map with all the data should look like the following one:

5. Try Yourself Convert Layers’ CRS¶

Because we are going to be measuring distances within our layers, we need to change the layers’ CRS. To do this, we need to select each layer in turn, save the layer to a new one with our new projection, then import that new layer into our map.

You have many different options, e.g. you can export each layer as a new Shapefile, you can append the layers to an existing GeoPackage file or you can create another GeoPackage file and fill it with the new reprojected layers. We will show the last option so the training_data.gpkg will remain clean. But feel free to choose the best workflow for yourself.

1. Right click the roads layer in the Layers panel;

2. Click Export –> Save Features As…;

3. In the Save Vector Layer As dialog choose GeoPackage as Format;

4. Click on … of File name parameter and name the new GeoPackage as vector_analysis;

5. Change the Layer name as roads_34S;

6. Change the CRS parameter to WGS 84 / UTM zone 34S;

7. Finally click on OK:

   

   This will create the new GeoPackage database and fill it with the roads_34S layer.

8. Repeat this process for each layer, creating a new layer in the vector_analysis.gpkg GeoPackage file with _34S appended to the original name and removing each of the old layers from the project.

   Note

   When you choose to save a layer to an existing GeoPackage, QGIS will append that layer in the GeoPackage.

9. Once you have completed the process for each layer, right click on any layer and click Zoom to layer extent to focus the map to the area of interest.

Now that we have converted OSM’s data to a UTM projection, we can begin our calculations.

6. Follow Along: Analyzing the Problem: Distances From Schools and Roads¶

QGIS allows you to calculate distances from any vector object.

1. Make sure that only the roads_34S and buildings_34S layers are visible, to simplify the map while you’re working

2. Click on the Processing ‣ Toolbox to open the analytical core of QGIS. Basically: all algorithms (for vector and raster) analysis are available within this toolbox.

3. We start by calculating the area around the roads_34S by using the Buffer algorithm. You can find it expanding the Vector Geometry group.

   

   Or you can type buffer in the search menu in the upper part of the toolbox:

   

4. Double click on it to open the algorithm dialog

5. Set it up like this

   

6. The default Distance is in meters because our input dataset is in a Projected Coordinate System that uses meter as its basic measurement unit. You can use the combo box to choose other projected units like kilometers, yards, etc.

   Note

   If you are trying to make a buffer on a layer with a Geographical Coordinate System, Processing will warn you and suggest to reproject the layer to a metric Coordinate System.

7. By default Processing creates temporary layers and adds them to the Layers panel. You can also append the result to the GeoPackage database by:

   1. clicking on the … button and choose Save to GeoPackage…
   2. naming the new layer roads_buffer_50m
   3. and saving it in the vector_analysis.gpkg file
   

8. Click on Run and then close the Buffer dialog.

Now your map will look something like this:

If your new layer is at the top of the Layers list, it will probably obscure much of your map, but this gives you all the areas in your region which are within 50m of a road.

However, you’ll notice that there are distinct areas within your buffer, which correspond to all the individual roads. To get rid of this problem:

1. Uncheck the roads_buffer_50m layer and re-create the buffer using the settings shown here:

   

   Note that we’re now checking the Dissolve result box

2. Save the output as roads_buffer_50m_dissolved

3. Click Run and close the Buffer dialog again

Once you’ve added the layer to the Layers panel, it will look like this:

Now there are no unnecessary subdivisions.

7. Try Yourself Distance from schools¶

Use the same approach as above and create a buffer for your schools.

It needs to be 1 km in radius. Save the new layer in the vector_analysis.gpkg file as schools_buffer_1km_dissolved.

8. Follow Along: Overlapping Areas¶

Now we have areas where the road is 50 meters away and there’s a school within 1 km (direct line, not by road). But obviously, we only want the areas where both of these criteria are satisfied. To do that, we’ll need to use the Intersect tool. You can find it in Vector Overlay group within Processing ‣ Toolbox.

1. Set it up like this:

   
   • The input layers are the two buffers
   • The saving location is, once again, the vector_analysis.gpkg GeoPackage
   • And the output layer name is road_school_buffers_intersect

2. Click Run.

   In the image below, the blue areas show us where both distance criteria are satisfied at once!

   

3. You may remove the two buffer layers and only keep the one that shows where they overlap, since that’s what we really wanted to know in the first place:

   

9. Follow Along: Extract the Buildings¶

Now you’ve got the area that the buildings must overlap. Next, you want to extract the buildings in that area.

1. Look for the menu entry Vector Selection ‣ Extract by location within Processing ‣ Toolbox

2. Set up the algorithm dialog like in the following picture

   

3. Click Run

4. You will see that the algorithm reports an error "Feature (1) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.". That means that the layer buildings_34S has features with problematic geometry. In order to proceed, we need to fix the problematic geometry.

   

5. First, let's analyse what features have wrong geometry. Temporarily close Extract by location window. Look for the menu entry Vector geometry ‣ Check validity, select the buildings layer, and run the algorithm.

   

6. As a result, the algorithm will report that it has found three invalid features and three new layers will appear in the list of layers: Valid output, Invalid output, and Error output.

   

7. Invalid output is of more importance for us, so remove the other two layers. Invalid output layer contains all features that have invalid geometry. Select some appropriate distinct color for this layer and put it on top to make it visible. There are only two features in this layer, so you can zoom to either of them.

   

8. Now, if you will open Attributes table of the Invalid output layer and scroll to the right, you will find a new column in the end with name errors. This is a column added by the Check validity algorithm and it is intended to clarify what is exactly wrong with the geometry. For both of our features it says "Nested shells". This means that one polygon of a multi-polygon feature is completely within another polygon (please note that the buildings layer has multi-polygon geometry type, which means that one feature can consist of several separate polygons).

   

9. If you enable editing mode and select the feature, you'll probably find that the feature indeed consists of two different polygons: the inner part is completely contained within the outer part.

   

10. There are several possibilities of how to fix the invalid geometry. One approach would be to manually remove the inner part or make a hole in the outer part with the shape of the inner part. However, this can be cumbersome if there are many features with invalid geometry. In this scenario, you may use Vector geometry ‣ Fix geometries algorithm that will automatically fix all features with invalid geometry (however, without guarantees that the result will satisfy you). Find the algorithm in the Processing toolbox and set it like in the picture below.

    

11. Now, you have a new buildings layer with fixed geometries. If you check, you will find that the Fix geometries algorithm simply cutted holes and removed inner polygons in the features with invalid geometry. While generally this might be not appropriate, for the purpose of our task it is good enough. You can now remove buildings_34S from layer list.

    

12. Now you can run Extract by location algorithm again, but this time using the buildings layer with fixed geometry. If you see "No spatial index exists for input layer, performance will be severely degraded", you may go into menu Vector ‣ Data Management Tools ‣ Create Spatial Index..., select the fixed buildings layer in the window and run the algorithm to build the index to improve the speed of QGIS.

    

13. After running the Extract by location algorithm again, you’ll probably find that not much seems to have changed. If so, move the well_located_houses layer to the top of the layers list, then zoom in.

    

    The red buildings are those which match our criteria, while the buildings in green are those which do not.

14. Now you have two separated layers and can remove buildings_34S_fixed from layer list.

10. Try Yourself Further Filter our Buildings¶

We now have a layer which shows us all the buildings within 1km of a school and within 50m of a road. We now need to reduce that selection to only show buildings which are within 500m of a restaurant.

Using the processes described above, create a new layer called houses_restaurants_500m which further filters your well_located_houses layer to show only those which are within 500m of a restaurant.

11. Follow Along: Select Buildings of the Right Size¶

To see which buildings are of the correct size (more than 100 square meters), we first need to calculate their size.

1. Select the houses_restaurants_500m layer and open the Field Calculator by clicking on the button in the main toolbar or within the attribute table

2. Set it up like this

   

   We are creating the new field AREA that will contain the area of each building square meters.

3. Click OK. The AREA field has been added at the end of the attribute table.

4. Click the edit mode button again to finish editing, and save your edits when prompted.

5. Build a query as earlier in this lesson

   

6. Click OK.

Your map should now only show you those buildings which match our starting criteria and which are more than 100m squared in size.

12. Try Yourself¶

Save your solution as a new layer, using the approach you learned above for doing so. The file should be saved within the same GeoPackage database, with the name solution.

1. Follow Along: The Tools and the Data¶

You can find all the network analysis algorithms in the Processing ‣ Network Analysis menu. You can see that there are many tools available:

Open the project from this zip-archive, it contains two layers:

1. network_points
2. network_lines

As you can see the network_lines layer has already a style that helps to understand the road network.

The shortest path tools provide ways to calculate either the shortest or the fastest path between two points of a network, given:

• start point and end point selected on the map
• start point selected on the map and end points taken from a point layer
• start points taken from a point layer and end point selected on the map

Let’s start.

2. Calculate the shortest path (point to point)¶

The Network analysis ‣ Shortest path (point to point) allows you to calculate the shortest distance between two manually selected points on the map.

In this example we will calculate the shortest (not fastest) path between two points.

In the following image we choose these two points as starting and ending point for the analysis:

1. Open the Shortest path (point to point) algorithm

2. Select network_lines for Vector layer representing network

3. Let Shortest in the Path type to calculate parameter

4. Click on the … button next to the Start point (x, y) and choose the location tagged with Starting Point in the picture. The menu is filled with the coordinates of the clicked point.

5. Do the same thing but choosing the location tagged with Ending point for End point (x, y)

6. Click on the Run button:

   
   
   

7. A new line layer is created representing the shortest path between the chosen points. Uncheck the network_lines layer to see the result better:

   
   
   

8. Let’s open the attribute table of the output layer. It contains three fields, representing the coordinates of the starting and ending points and the cost.

   We chose Shortest as Path type to calculate, so the cost represent the distance, in layer units, between the two locations.

   In our case, the shortest distance between the chosen points is around 1000 meters:

   
   
   

Now that you know how to use the tool, feel free to change them and test other locations.

3. Try Yourself Fastest path¶

With the same data of the previous exercise, try to calculate the fastest path between the two points.

How much time do you need to go from the start to the end point?

4. Follow Along: Advanced options¶

Let’s explore some more options of the Network Analysis tools. In the previous exercise we calculated the fastest route between two points. As you can imagine, the time depends on the travel speed.

We will use the same layers and same starting and ending points of the previous exercises.

1. Open the Shortest path (point to point) algorithm

2. Fill the Input layer, Start point (x, y) and End point (x, y) as we did before

3. Choose Fastest as the Path type to calculate

4. Open the Advanced parameter menu

5. Change the Default speed (km/h) from the default 50 value to 4

   
   
   

6. Click on Run

7. Once the algorithm is finished, close the dialog and open the attribute table of the output layer.

   The cost field contains the value according to the speed parameter you have chosen. We can convert the cost field from hours with fractions to the more readable minutes values.

8. Open the field calculator by clicking on the icon and add the new field minutes by multiplying the cost field by 60:

   
   
   

That’s it! Now you know how many minutes it will take to get from one point to the other one.

5. Shortest map with speed limit¶

The Network analysis toolbox has other interesting options. Looking at the following map:

we would like to know the fastest route considering the speed limits of each road (the labels represent the speed limits in km/h). The shortest path without considering speed limits would of course be the purple path. But in that road the speed limit is 20 km/h, while in the green road you can go at 100 km/h!

As we did in the first exercise, we will use the Network analysis ‣ Shortest path (point to point) and we will manually choose the start and end points.

1. Open the Network analysis ‣ Shortest path (point to point) algorithm

2. Select network_lines for the Vector layer representing network parameter

3. Choose Fastest as the Path type to calculate

4. Click on the … button next to the Start point (x, y) and choose the location tagged with Start Point in the picture. The menu is filled with the coordinates of the clicked point.

5. Do the same thing but choosing the location tagged with End point for End point (x, y)

6. Open the Advanced parameters menu

7. Choose the speed field as the Speed Field parameter. With this option the algorithm will take into account the speed values for each road.

   
   
   

8. Click on the Run button

9. Turn off the network_lines layer to better see the result

   
   
   

As you can see the fastest route does not correspond to the shortest one.

6. Service area (from layer)¶

The Network Analysis ‣ Service area (from layer) algorithm can answer the question: given a point layer, what are all the reachable areas given a distance or a time value?

Given a distance of 250 meters we want to know how far we can go on the network from each point of the network_points layer.

1. Uncheck all the layers except network_points

2. Open the Network Analysis ‣ Service area (from layer) algorithm

3. Choose network_lines for Vector layer representing network

4. Choose network_points for Vector layer with start points

5. Choose Shortest in Path type to calculate

6. Enter 250 in the Travel cost parameter

7. Click on Run and then close the dialog

   
   
   

   The output layer represents the maximum path you can reach from the point features given a distance of 250 meters:

   
   
   

Overview of the task¶

Given the locations of all known significant earthquakes, find out the nearest populated place for each location where the earthquake happened.

Get the data¶

We will use NOAA’s National Geophysical Data Center’s Significant Earthquake Database as our layer representing all major earthquakes. Download the tab-delimited earthquake data.

Natural Earth has a nice Populated Places dataset. Download the simple (less columns) dataset

For convenience, you may directly download a copy of both the datasets from the links below:

signif.txt

ne_10m_populated_places_simple.zip

Data Sources: [NGDC] [NATURALEARTH]

Procedure¶

1. Locate the downloaded ne_10m_populated_places_simple.zip file in the Browser panel and expand it. Drag the ne_10m_populated_places_simple.shp file to the canvas.

2. You will see a new layer ne_10m_populated_places_simple loaded in the Layers panel. This layer contains the points representing populated places. Now we will load the earthquakes layer. This layer comes as a Tab Serepated Values (TSV) text file. To load this file, click the Open Data Source Manager button on the Data Source Toolbar. You can also use Ctrl + L keyboard shortcut.

3. Click the ... button next to File name and browse to the downloaded signif.txt file. Once loaded, the File Format and Geometry Definition fields should be auto-populated with correct values. Click Add followed by Close.

4. Zoom around and explore both the datasets. Each yellow point represents the location of a significant earthquake and each red point represents the location of a populated place. Our goal is to find out the nearest point from the populated places layer for each of the points in the earthquake layer.

5. Before we do the analysis, we need to clean up our inputs. The signif layer contains many records without a valid geometry. These records were imported with a NULL geometry. So let’s remove these records first. Go to Processing ‣ Toolbox.

6. Search for and locate the Vector geometry ‣ Remove null geometries tool. Double-click to launch it.

7. Select signif as the Input layer and click Run. Once the processing finishes, click Close.

8. You will see a new layer caled Non null geometries loaded into the Layers panel. We will use this layer instead of the original signif layer in further analysis. Un-check the signif layer in the Layers panel to hide it. Now it is time to perform the nearest neighbor analysis. Search and locate the Vector analysis ‣ Distance to nearest hub (line to hub) tool. Double-click to launch it.

9. In the Distance to Nearest Hub (Line to Hub) dialog, select Non null geometries as the Source points layer. Select ne_10m_populated_places_simple as the Distination hubs layer. Select name as the Hub layer name attribute. The tool will also compute straight-line distance between the populated place and the nearest earthquake. Set Kilometers as the Measurement unit. Click Run. Once the processing finishes, click Close.

10. Back in the main QGIS window, you will see a new line layer called Hub distance loaded in the Layers panel. This layer has line features connecting each earthquake point to the nearest populated place. Right-click the Hub distance layer and select Open Attribute Table.

11. Scroll right to the last columns and you will see 2 new attributes called HubName and HubDist added to the original earthquake features. This is the name the distance to the nearest neighbor from the populated places layer.

12. Our analysis is complete. We can now explore another tool that can also do a similar analysis. Distance Matrix is a powerful tool that allows you to not only compute distance to the nearest point, but to all the points from another layer. We can use this method as an alternative to the Distance to nearest hub tool. Un-check the Hub distance layer to hide it. Search and locate the Vector analysis ‣ Distance matrix tool.

13. In the Distance matrix dialog, set Non null geometries as the Input point laeyer and I_D as the Input unique ID field. Set ne_10m_populated_places_simple as the Target point layer and name as the Target unique ID field. Select Linear (N*k x 3) distance matrix as the Output matrix type. The key here is to set the Use only the nearest (k) target points parameter to 1 - which will give you only the nearest neighbor in the output. Click Run to start the matrix calculation. Once the processing finishes, click Close.

14. Once the processing finishes, a new layer called Distance matrix will be loaded. Note that the output of this tool is a layer containin MultiPoint geometries. Each feature contains 2 points - source and target. Open the Attribute Table for the layer. You will see that each feature has attributes mapping the earthquake to its nearest populated place. Note that the distance here is in the layer’s CRS units (degrees).

15. At this point, you can save your results in the format of your choice by right-clicking the layer and selecting Export ‣ Save Features As. If you want to visualize the results better, we can easily create a hub-spoke rendering from the feature’s geometry. Right-click the Distance matrix layer and select Properties.

16. In the Properties dialog, switch to the Symbology tab. Click on the Simple marker sub-renderer and select Geometry generator as the Symbol layer type. Set LineString / MultiLineString as the Geometry type. Click the Expression button.

17. Here we can enter an expression to create a line geometry from the 2 points within each multi-point source geometry. Enter the following expression.

make_line(point_n( $geometry, 1), point_n( $geometry, 2))
    

18. Back in the Symbology tab, set the style of the line as per your liking and click OK.

19. You will see the Distance matrix layer now rendered with lines instead of points. Note that we did not have to create a new layer for this visualization. The layer still contains MultiPoint geometries, but it is dynamically rendered as lines based on the expression.