Overview of the task¶

We will take several global vector layers and clip them to the extent of Africa in a single batch command.

Get the data¶

Natural Earth has several global vector layers. Download the following layers

• Admin 0 - Countries
• Railroads
• Ports
• Urban Areas

For convenience, you may directly download a geopackage containing the above layers from below:

ne_global.gpkg

Data Source: [NATURALEARTH]

Procedure¶

1. In the QGIS Browser Panel, locate the directory where you saved your downloaded data. Expand the zip or the gpkg entry and select the ne_10m_admin_0_countries layer. Drag the layer to the canvas.

2. You will see the layer loaded in the Layers panel. As our task is to clip the global layers to the boundary of Africa, we need to first prepare a layer containg features only from that continent. Let’s look at the attribute table to see what column can be used to query features belonging to a particular continent. Right-click the ne_10m_admin_0_countries layer and select Open Attribute Table.

3. In the Attribute Table window, as you scroll horizontally, you will see that the data contains an attribute called CONTINENT. We can use this attribute to filter features.

4. Close the attribute table and return to the main QGIS window. Right-click the ne_10m_admin_0_countries layer and select Filter.

5. In the Query Builder window, select the CONTINENT field and click Sample. This will populate the Values panel with the a subset of values of that attribute from features. This step is useful to get an understanding of what type of values are present in your dataset. We can see that our dataset contains contains a value called Africa among others.

6. Now enter the expression in the Provider specific filter expression textbox. You can click the CONTINENT label, followed by = button and Africa label. Or you can type the following expression in the textbox. Click OK after entering the expression.

"CONTINENT" = 'Africa'
              

7. You will see that the map canvas now contains only the countries from Africa. Note the filter icon next to the layer name indicating that that layer a filter applied to it. If you wanted to see and use all the features from the layer, you can click the filter icon and clear the expression. For now, we will keep the filter so we can clip other layers to Africa.

8. We are now ready to run the batch process to clip the layers. Locate the natural earth global layers ne_10m_railroads, ne_10m_ports and ne_10m_urban_areas in the QGIS Browser panel. Hold the Ctrl key and click each layer to select them. Once selected, drag them to the canvas.

9. Once the layers are loaded, you will notice that they are global layers and have features spanning all the countries. Now, it’s time to start our batch clip process to clip these layers to Africa. Open Processing ‣ Toolbox.

10. Browse all available algorithms and find the Clip tool from Vector overlay ‣ Clip. You may also use the Search box to easily find the algorithm as well.

11. Right-click the Clip algorithm and select Execute as Batch Process.

12. In the Batch Processing dialog, the first tab is Parameters where we define our inputs. Each row in the table represents 1 processing task. Click Add row button to add a new row. As our task involves 3 layers, add 3 rows.

13. Click the ... next to the first row in the Input layers column. Select Select from Open Layers.

14. In the Multiple selection dialog, check the 3 layers that we want to clip and click OK.

15. You will notice that the Input layer columns will be auto-populated with all layers you had selected. Next, we need to select the layer containing the boundary to clip our input layers. Click the ... button for the first row under Overlay layer column and select Select form Open Layers.

16. In the Multiple selection dialog, check ne_10m_admin_0_countries and click OK.

17. Since the clip layer is the same for all our inputs, a handy shortcut is to double-click the column header Overlay layer and the same layer will be auto-filled for all the rows.

18. Next, we need to define our outputs. Click the ... buton next to the first row in the Clipped column.

19. Browse the the directory where you want your output layers. Type the filename as clipped_ and click Save.

20. You will see a new Autofill settings dialog pop up. Select Fill with parameter values as the Autofill mode. Select Input layer as the Parameter to use. This setting will add the input file name to the output along with the specified clipped_ filename. This is important to ensure all the output files have unique names and they do not overwrite each other.

21. Now we are ready to start the batch procesing. Make sure to check Load layers on completion and click Run.

22. The clip algorithm will run for each of the inputs and create output files are we have specified. Once the batch process finishes, click Close to return to QGIS.

23. Back in the main QGIS window, you will see the layers added to QGIS canvas. As you will notice, all the global layers are properly clipped to the continent boundary that we had specified.

24. We have accomplished the task of clipping multiple layers in a batch. But QGIS3 has another handy feature that will help you save and deliver the result in a more efficient way. If you wanted to deliver the clipped layers to someone, you would zip the individual files outside of QGIS. A better option to to package the output layers in a single Geopackage. In the QGIS Browser, locate the clipped output layers. You may have to click the Refresh button to see the newly added files.

25. Hold Ctrl key and select the layers. Drag them to the canvas to load them in QGIS.

25. In Processing Toolbox, locate the Database ‣ Package layers tool.

27. In the Package layers dialog, click the ... button next to Input layers.

28. In the Multiple selection dialog, check the clipped_ne_10m_ports, clipped_ne_10m_railroads and clipped_ne_10m_urban_areas layers. Click OK.

29. Once the input layers are selected, click the ... next to Destination Geopackage and choose Save To File.

30. Enter the output file name as clipped_layers.

31. Click Run to start the packaging process.

32. Once the process finishes, you will see a new geopackage file in your QGIS Browser containing all the clipped output layers. This is a single file on your computer that contains all the output layers.

Overview of the task¶

We will take a point layer of maritime piracy incidents and create a processing model to produce a density map by aggregating them over a global hexagonal grid.

Get the data¶

National Geospatial-Intelligence Agency’s Maritime Safety Information portal provides a shapefile of all incidencts of maritine piracy in the form on Anti-shipping Activity Messages. Download the Arc Shape file version of the database.

Natural Earth has several global vector layers. Download the 10m Physical Vectors - Land containing Land polygons.

For convenience, you may directly download a copy of the above layers from below:

ASAM_shp.zip

ne_10m_land.zip

Data Source: [NGA_MSI] [NATURALEARTH]

Procedure¶

1. In the QGIS Browser Panel, locate the directory where you saved your downloaded data. Expand the ne_10m_land.zip and select the ne_10m_land.shp layer. Drag the layer to the canvas. Next, locate the ASAM_shp.zip file. Expand it and select the asam_data_download/ASAM_events.shp layer and drag it on to the canvas.

2. The ASAM_events.shp layer does not have projection information associated with it, so you will be prompted to select a CRS in the Coordinate Reference System Selector. Here, the points are in the Latitude and Longitude coordinates, so select the WGS 84 CRS and click OK.

3. Once the layer is loaded, you can see the individual points representing incidents of piracy locations. Let’s start building our Processing model to process these layers. Go to Processing ‣ Graphical Modeler....

4. In the Processing Modeler dialog, locate the Model Properties panel. Enter piracy hexbin as the Name of the model and projects as the Groups. Click the Save button.

5. Save the model as piracy_hexbin.

6. Now we can start building a graphical model of our processing pipeline. The Processing modeler dialog contains a left-hand panel and a main canvas. On he left-hand panel, locate the Inputs panel listing various types of input data types. Scroll down and select the + Vector Layer input. Drag it to the canvas.

7. Enter Input Points as the Parameter name and Point as the Geometry type. This input represents the piracy incidents point layer.

8. Next, drag another + Vector Layer input to the canvas. Enter Base Layer as the Parameter name and Polygon as the Geometry type. This input represents the natural earth global land layer.

9. As we are generating a global hexagonal grid, we can ask the user to supply us the grid size as an input instead of hard-coding it as part of our model. This way, the user can quickly experiment with different grid sizes without changing the model at all. select a + Number input and drag it to the canvas. Enter Grid Size as the Parameter name and click OK.

10. Now that we have our user inputs defined, we are ready to add processing steps. All of the processing algorithms are available to you under the Algorithms tab. The first step in our pipeline will be to reproject the base layer to the Project CRS. Search for Reproject layer algorithm and drag it to the canvas.

11. In the Reproject layer dialog, select Base Layer as the Input layer. Check the Use project CRS as the Target CRS. Click OK.

12. In the Processing Modeler canvas, you will notice a connection appear between the + Base Layer input and the Reproject layer algorithm. This connection indicates the flow of our processing pipeline. Next step is to create a hexagonal grid. Search for the Create grid algorithm and drag it to the canvas.

13. In the Generate grid dialog, choose Hexagon (polygon) as the Grid type. Select Extent of 'Reprojected' from algorithm 'Reproject Layer' as the Grid extent. Click the 123 button under the Horizonal spacing label and choose Model input.

14. Select Grid Size input for Using model input. Repeat the same process for Vertical Spacing. Click OK.

15. At this point, we have a global hexagonal grid. The grid spans the full extent of the base layer, including land areas and places where there are no points. Let’s filter out those grid polygons where there are no input points. Search for Extract by location algorithm and drag it to the canvas.

16. For Extract features from, select 'Grid' from algorithm 'Generate Grid', Where the features (geometric predicate) as Intersect and By compating to the features from as Input points. Click OK.

17. Now we have only those grid polygons that contain some input points. To aggregate these points, we will use Count points in polygon algorithm. Search and drag it to the canvas.

18. Select 'Extracted (location)' from algorithm 'Extract by location' as the value for Polygons. The Points layer would be Input Points. At the bottom, name the Count output layer as Aggregated. Click OK.

19. The model is now complete. Click the Save button.

20. Switch to the main QGIS window. You can find your newly created model in the Processing Toolbox under Models ‣ projects ‣ piracy_hexbin. Now it is time to run and test the model. As our goal is to aggregate the input points over hexagonal grids, it is important that the grids are generated using a equal-area projection. This will ensure that regardless of the location of the grid, it will cover exactly the same area. Our model doesn’t explicitely ask for a CRS, but uses whatever CRS is set as the Project CRS. Let’s choose a global equal area projection as the Project CRS. Go to Project ‣ Properties.

21. In the Project Properties dialog, switch to the CRS tab. We will use a global Mollweide projection for this exercise which is a equal area projection. Search for Mollweide in the Filter box and select World_Mollweide EPSG:54009 as the CRS. Click OK.

22. You will see the layers getting reprojected on-the-fly to the selected CRS. Locate the piracy_hexbin model in the Processing Toolbox and double-click it.

23. Our Base Layer is the ne_10m_land and the Input Points layer is ASAM_events. The Grid Size needs to be specified in the units of the selected CRS. The World_Mollweide CRS unit is meters, so we specify 100000 m (100 Kms) as the Grid Size. Click Run to start the processing pipeline. Once the process finishes, click Close.

24. You will see a new layer Aggregated loaded as the result of the model. As you explore, you will notice the layer contains an attribute called NUMPOINTS containing the number of piracy incidents points contained within that grid feature. Let’s style this layer to display this information better. Right-click the Aggregated layer and select Properties.

25. Switch to the Symbology tab. Select Graduated symbology and NUMPOINTS as the Column. Click Change.. next to Symbol label.

26. Select Simple fill symbol and check the Transparent Stroke button under Stroke color. This is to make the hexagon edges transparent.

27. Click the dropdown next to Color ramp and select the Viridis ramp. Click the dropdown again and select Invert Color Ramp to reverse the order of color.

28. The Graduated symbology will divide the values in the selected column into distinct classes and assign a different color to each of the classes. Select Natural Breaks (Jenks) as the Mode and click Classify and click OK.

29. Back in the main QGIS window, turn off the ASAM_events layer. You will see a nice visualization of piracy hotspots across the globe.

Now that you have encoded the full data pipeline in the model, it is easy to reproduce your results. A model also allows you to experiment quickly without manually repeating each intermediate step every time. If your inputs change over time, say an updated database of piracy is released after a few months, you can run your model on that input to generate a similar visualization without having to remember each step.

Overview of the task¶

In this tutorial, we will identify the suitable areas for development, that are

• Close to roads, and
• Away from water bodies, and
• Not in a protected region.

Get the data¶

We will use vector data layers from OpenStreetMap (OSM). OSM is a global database of freely available base map data. Geofabrik provides daily updated shapefiles of OpenStreetMap datasets.

We will be using the OSM data layers for the state of Assam in India. Geofabrik India shapefiles were downloaded, clipped to the state boundary and packaged in a single GeoPackage file. You can download a copy of the geopackage from the link below:

assam.gpkg

Data Source: [GEOFABRIK]

Procedure¶

1. Browse to the downloaded assam.gpkg file in QGIS Browser. Expand it and drag each of the 5 individual data layers to the map canvas. You will see boundary, roads, protected_regions, water_polygons and water_polylines layers loaded in the Layers panel.

2. First step in the overlay analysis, is to convert each data layer to raster. An important consideration is that all rasters must be of the same extent. We will use the boundary layer as the extent for all the rasters. Go to Processing ‣ Toolbox. Search for and locate the GDAL ‣ Vector conversion ‣ Rasterize (vector to raster)` algorithm. Double-click to launch it.

3. In the Rasterize (vector to raster) dialog, select roads as the Input layer. We want to create an output raster where pixel values are 1 where there is a road and 0 where there are no roads. Enter 1 as the A fixed value to burn. The input layers are in a projected CRS with meters are the unit. Select Geoferenced units as the Output raster size units. We will set the resolution of the output raster to be 15 meters. Select 15 as both Width/Horizontal resolution and Height/Vertical resolution. Next, click the ... button next to Output extent and select boundary for Use extent for.

4. Scroll down further and click the arrow button in the Assign a specific nodata value to output bands. That field should now be set to Not set. This is important because when raster calculator (which we will use later) encounters a pixel with nodata value in any layer, it sets the output to nodata as well, resulting is wrong output. Set the Rasterized output raster as raster_roads.tif and click Run.

5. Once the processing finishes, you will see a new layer raster_roads loaded in the Layers panel. The default styling will show pixels with roads as white and others as black. We want to convert other 4 vector layers to rasters as well. Rather than running the rasterize algorithm one-by-one, we can use the built-in batch-processing functionality to convert them all at once. See batch-processing tutorial to learn more about batch processing. Right-click the Rasterize (vector to raster) algorithm and select Execute as Batch Process.

6. In the Batch Processing dialog, click the ... button in the first row of the Input layer column. Select boundary, protected_regions, water_polygons and water_polylines layers and click OK.

7. Fill in the parameters with the same values we used in the roads layer. After filling the first-row of the parameter, use the Autofill ‣ Fill Down button to add the same value for all layers.

8. In the last column Rasterized, click the ... button in the first row. Select Fill with parameter values as the Autofill mode and Input layer as the Parameter to use. Click OK.

9. Browse to a directory on your computer and name the layer as raster_. The batch processing interface will autocomplete the name with the layer name and fill in all rows. Make sure the Load layers on completion box is checked and click OK.

10. Once the processing finishes, you will have 4 new raster layers loaded in the Layers panel. You will notice that we have 2 water related layers - both representing water. We can merge them to have a single layer representing water areas in the region. Search for and locate Raster analysis ‣ Raster calculator algorithm in the Processing Toolbox. Double-click to launch it.

11. Enter the following expression in the Expression box. You can click on the appropriate layer in the Layers box to auto insert the layer names. This expression means that we want to sum the pixel values in the first band of both the water rasters. Click the ... button next to Reference layer(s) and select ``raster_water_polygons` as the reference layer. Name the output raster_water_merged.tif and click Run.

"raster_water_polygons@1" + "raster_water_polylines@1"
    

12. The resulting merged raster will have pixels with value 1 for all areas with water. But you will notice that there are some regions where there was both a water polygon and a water polyline. Those areas will have pixels with value 2 - which is not correct. We can fix it with a simple expression. Open Raster analysis ‣ Raster calculator algorithm again.

13. Enter the following expression which will assign the value 1 that match the expression and 0 where it doesn’t. Click the ... button next to Reference layer(s) and select ``raster_water_merged` layer. Name the output raster_water.tif and click Run.

"raster_water_merged@1" > 0
    

14. The resulting layer raster_water now has pixels with only 0 and 1 values.

15. Now that we have layers representing road and water pixels, we can generate proximity rasters. These are also known as Euclidean distances - where each pixel in the output raster represents the distance to the nearest pixel in the input raster. This resulting raster can be then used to determine suitable areas which are within certain distance from the input. Search for and locate the GDAL ‣ Raster analysis ‣ Proximity (raster distance) algorithm. Double-click to launch it.

16. In the Proximity (Raster Distance) dialog, select raster_roads as the Input layer. Choose Georeferenced coordinates as the Distance units. As the input layers are in a projected CRS with meters as the units, enter 5000 (5 kilometers) as the Maximum distance to be generated. Make sure the Nodata value to use for the destination proximity raster value is Not set. Name the output file as roads_proximity.tif and click Run.

17. Once the processing is over, a new layer roads_proximity will be added to the Layers panel. To visualize it better, let’s change the default styling. Click the Open the Layer Styling panel button in the Layers panel. Change the Max value to 5000 under Color gradient.

18. Repeat the Proximity (Raster Distance) algorithm for the raster_water layer with same parameters and name the output water_proximity.tif.

19. Once the processing finishes, you can apply the similar styling as before to visualize the results better. If you click around the resulting raster, you will see that it is a continuum of values from 0 to 5000. To use this raster in overlay analysis ,we must first re-classify it to create discrete values. Open Raster analysis ‣ Raster calculator algorithm again.

20. We want to give higher score to pixels that are near to roads. So let’s use the following scheme.

• 0-1000m –> 100

• 1000-5000m –> 50

• >5000m –> 10

  Enter the following expression that applies the above criteria on the input. Click the ... button next to Reference layer(s) and select ``roads_proximity` layer. Name the output roads_reclass.tif and click Run.

  100*("roads_proximity@1"<=1000) + 50*("roads_proximity@1">1000)*("roads_proximity@1"<=5000) + 10*("roads_proximity@1">5000)
      

  

21. Once the re-classification process finishes, a new layer roads_proximity will be added to the Layers panel. This layer has only 3 different values, 10, 50 and 100 indicating relative suitability of the pixels with regards to distance from roads. Open Raster analysis ‣ Raster calculator algorithm again.

22. Repeat the re-classification process for the water_proximity layer. Here the scheme will be reverse, where pixels that are further away from water shall have higher score.

• 0-1000m –> 10

• 1000 -5000m —> 50

• >5000m –> 100

  Enter the following expression that applies the above criteria on the input. Click the ... button next to Reference layer(s) and select ``water_proximity` layer. Name the output water_reclass.tif and click Run.

  100*("water_proximity@1">5000) + 50*("water_proximity@1">1000)*("water_proximity@1"<=5000) + 10*("water_proximity@1"<1000)
      

  

23. Now we are ready to do the final overlay analysis. Recall that our criteria for determining suitability is as follows - close to roads, away from water and not in a protected region. Open Raster analysis ‣ Raster calculator. Enter the following expression that applies these criteria. Note that we are multiplying the result with raster_boundary@1 at the end to discard pixel values outside of the state boundary. Click the ... button next to Reference layer(s) and select raster_boundary layer. Name the output overlay.tif and click Run.

("roads_reclass@1" + "water_reclass@1")*("raster_protected_regions@1"  !=  1 )*"raster_boundary@1"
    

24. Once the processing finishes, the resulting raster overlay will be added to the Layers panel. The pixel values in this raster range from 0 to 200 - where 0 is the least suitable and 200 is the most suitable area for development. Click the Open the Layer Styling panel button in the Layers panel.

25. Select singleband_pseudocolor renderer and the Spectral color ramp. Click Classify to apply the color ramp to the raster.

26. Click on the default label values next to each color and enter appropriate labels. The labels will also appear as the legend under the overlay layer.

27. Raster layers are rectangular grids. We want to hide pixels outside the state boundary. An easy way to achieve this is applying an Inverted Polygons rendered to the vector boundary layer. Scroll down in the Layers panel and locate the boundary layer. Select Inverted Polygons as the renderer and leave other options to default.

28. For the effect of the renderer to show, it needs to be at the top of the Table of Contents. Right-click the boundary layer and select Move to Top.

29. Check the layer and the map canvas would update to show the overlay raster clipped to the boundary layer. This is the final output that shows areas within the state that are suitable for development.