In parts 1—6 of this series we worked on implementing the famous boids model in Python and ended up with a pair of scripts that output a CSV file of boid locations which we can then import into QGIS. What if I told you there was a better way? QGIS plugins can be developed in Python, which is mostly why we are using one of my least-preferred languages here, and in this section we will set up a basic implementation of our earlier model to be run from QGIS and output our results as a memory layer directly into our map. We’ll follow Ujaval Gandhi’s excellent tutorial on setting up a Python plugin for QGIS using the plugin builder (GeoApt LLC).
Getting started
Follow the tutorial linked above to get started with building a plugin, but stop just after copying our plugin to the profile directory. You should now have a plugin which loads in QGIS but only gives you a blank dialog when run. Copy the boid.py (or whatever you named your module from the previous chapters) file into the directory alongside the plugin.
Add boid.py to the end of the “python_files:” line in pb_tool.cfg the [files] section of pb_tool.cfg:
. . . [files] # Python files that should be deployed with the plugin python_files: __init__.py ``your_plugin``.py ``your_plugin``_dialog.py boid.py . . .
Recompile the directory using pb_tool compile. This should be the only time you need to do this step.
Import the Boid object into your plugin by adding an import statement into the top of “your_plugin“.py. Also import numpy:
import numpy as np from .boid import Boid
The structure of our plugin builder plugin
The plugin builder plugin will take a Qt GUI file and display it. This file can be used to set up our analysis. The way we set up this particular plugin, the main code will run when the OK button is clicked. This behaviour is implemented at the very end of the run method in “your_plugin“.py, in the block of code starting with if result:. You will see that this currently just contains pass.
Jamming our model in
The naive approach to running our model is to take everything we put in our main.py file into a new method in this file. We can then replace the pass at the end of the run method with a call to our method. I put mine right before the run method and called it runBoids:
. . .
def runBoids(self):
# Study area bounds:
. . .
f.open("positions.csv")
. . .
# Let's loop!
for tS in timeSteps:
for boid in school:
. . .
. . .
def run(self):
. . .
if result:
self.runBoids()
If we do this, however, we will have put a fancy GUI and new dependencies on our earlier script. We probably want to add our results to a layer directly into our project.
A mapping approach
We can load our results as a layer directly using PyQGS, QGS’s built-in Python API. Before we do this we need some to add some imports to the top of our file. Update your imports to look like this:
from qgis.PyQt.QtCore import (
QSettings,
QTranslator,
QCoreApplication,
QVariant
)
from qgis.PyQt.QtGui import (
QIcon,
QColor
)
from qgis.PyQt.QtWidgets import QAction
from qgis.core import (
Qgis, # Send messages to the user when the model run completes
QgsApplication, # Run parts of QGIS
QgsPointXY, # Use the point geometry
QgsProject, # Open the current project, add layers
QgsCoordinateReferenceSystem, # Get and set our CRS
QgsField, # Add fields to our points
QgsFields, # IBID
QgsFeature, # Add features to a layer
QgsGeometry, # Add geometry to features
QgsVectorDataProvider, # Create and modify a vector layer
QgsVectorLayer, # IBID
QgsVectorFileWriter # If we want to save our output
)
from qgis.core.additions.edit import edit
# Imports for our model
import numpy as np
from .boid import Boid
Creating a memory vector layer
The fastest way to implement our vector layer is as a memory layer. Don’t worry, it will still be slower than molasses in the northern hemisphere in January. This is Python after all. This risk of losing our results is present with a memory layer if QGIS crashes or we suffer a power outage, but what’s GIS without a little danger?
Creating a memory layer is pretty easy. We’ll simply create a vector layer object with geometry “Point”, name “Temporary_Boids”, and type “memory”
vLayer = QgsVectorLayer("Point", "Temporary_Boids", "memory")
In order to modify the vector layer we just created, we’ll need its data provider. We can access that directly from the layer object.
prov = vLayer.dataProvider()
QGIS will handle the geometry directly so we don’t need an X or Y field as we did in the CSV. There are still important fields that we might want for displaying our layer, however. We’ll add “boidID”, “Direction”, and “Time”. Fields are added as attributes the QgsField object takes two arguments: name and type. Our boidID will be a string, Direction will be a double, and Time will be an integer. We’ll then update the fields in our layer:
prov.addAttributes([QgsField("boidID", QVariant.String),
QgsField("Direction", QVariant.Double),
QgsField("Time", QVariant.Int)])
vLayer.updateFields()
Setting our layer CRS
By default our memory layer has no CRS. In my experience QGIS takes this to mean it is in WGS84 by default, and we’d have to set it at the end. For now let’s just use the project CRS (remember to set it first—from my examples we have been using WGS84 / UTM Zone 10N—EPSG:32610). We can get the CRS from the project using QgsProject.instance()’s crs() method, and set the layer’s CRS using its setCrs() method.
project = QgsProject.instance() crs = project.crs() . . . vLayer.setCrs(crs)
Adding features
Adding features to our vector map is a four step process. We’ll want to create a feature object, set the object’s geometry, and then add it’s attributes. Finally, we’ll add the feature to our vector layer’s data provider:
fet = QgsFeature() fet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(boid.position[0], boid.position[1]))) fet.settAttributes([str(boid.boidID), direction, tS]) prov.addFeatures([fet])
Now, we probably didn’t add a point within the extents of our vector layer, so it won’t be terribly useful. We’ll want to update the vector layer’s extents:
vLayer.updateExtents()
Adding the layer to our project
If you try these functions from QGIS, you might notice that stuff happens (ie, memory usage increases, QGIS locks up temporarily), but nothing shows up in our project. We still need to add our memory layer to our project to be able to see it. We’ll get our project instance, then check that our layer is valid, and then add it as a map layer to our project:
project = QgsProject.instance()
if vLayer.isValid():
project.addMapLayer(vLayer)
Implementing our control script in this method
Now we’ve covered our geometry options, we need to actually run our boids. This will function much the same way as our original main.py script, modified to use our QGIS geometry.
def runBoids(self):
# A simple implementation of boids, responsible for
# calling our boid object and updating it's position,
# then adding the boid positions to our map.
# Get our project instance
project = QgsProject.instance()
# Get our project CRS
crs = project.crs()
# Create our layer
vLayer = QgsVectorLayer("Point", "temporary_boids", "memory")
prov = vLayer.dataProvider()
# Add our fields
prov.addAttributes([QgsField("boidID", QVariant.String),
QgsField("Direction", QVariant.Double),
QgsField("Time", QVariant.Int)])
vLayer.updateFields()
# Set our CRS
vLayer.setCrs(crs)
# Establish our study area boundaries
xMin = 585706
yMin = 4073066
xMax = 594870
yMax = 4078604
# Study area size
xRange = xMax - xMin
yRange = yMax - yMin
# How many boids?
nB = 100
school = [Boid(np.random.random() * 10 + xMin + xRange/2,
np.random.random() * 10 + yMin + yRange/2,
xMin, yMin, xMax, yMax) for _ in range(nB)]
# How many time steps?
timeSteps = 1000
# Let's loop!
for tS in range(timeSteps):
for boid in school:
boid.behave(school)
boid.update()
direction = np.arctan2(boid.velocity[0], boid.velocity[1]) * 57.29578
# Create a feature and add it
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(boid.position[0], boid.position[1])))
fet.setAttributes([str(boid.boidID), direction, tS])
prov.addFeatures([fet])
vLayer.updateExtents()
# Add our map layer to the project
if vLayer.isValid():
project.addMapLayer(vLayer)
Taking control of our parameters
So we now have a working plugin that runs Boids directly in QGIS. We still need to change our model parameters in the code to change how the model runs. What if I want to run my boids somewhere other than the Soquel Canyon Marine Protected Area? We have a whole GUI feature we can take advantage of.
Changing our parameters
The first step is to change our runBoids method to take all of the important control information as parameters. We can then delete the calls to those parameters in the code itself. Don’t worry about sampleFreq just yet, we will cover that shortly:
def runBoids(self, nB, xMin, yMin, xMax, yMax, timeSteps, sampleFreq):
# A simple implementation of boids, responsible for
# calling our boid object and updating it's position,
# then adding the boid positions to our map.
# Get our project instance
project = QgsProject.instance()
# Get our project CRS
crs = project.crs()
# Create our layer
vLayer = QgsVectorLayer("Point", "temporary_boids", "memory")
prov = vLayer.dataProvider()
# Add our fields
prov.addAttributes([QgsField("boidID", QVariant.String),
QgsField("Direction", QVariant.Double),
QgsField("Time", QVariant.Int)])
vLayer.updateFields()
# Set our CRS
vLayer.setCrs(crs)
# Establish our study area boundaries
# xMin = 585706
# yMin = 4073066
# xMax = 594870
# yMax = 4078604
# Study area size
xRange = xMax - xMin
yRange = yMax - yMin
# How many boids?
# nB = 100
school = [Boid(np.random.random() * 10 + xMin + xRange/2,
np.random.random() * 10 + yMin + yRange/2,
xMin, yMin, xMax, yMax) for _ in range(nB)]
# How many time steps?
# timeSteps = 1000
# Let's loop!
for tS in range(timeSteps):
for boid in school:
boid.behave(school)
boid.update()
direction = np.arctan2(boid.velocity[0], boid.velocity[1]) * 57.29578
# Create a feature and add it
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(boid.position[0], boid.position[1])))
fet.setAttributes([str(boid.boidID), direction, tS])
prov.addFeatures([fet])
vLayer.updateExtents()
# Add our map layer to the project
if vLayer.isValid():
project.addMapLayer(vLayer)
Now we need to add all of these parameters to our call to runBoids at the end of the run method.
def run(self):
. . .
if result:
self.runBoids(nB, nX, nY, mX, mY, tS, sF)
This won’t run yet, however. We need to ask our user for these parameters. We’ll add these to our GUI.
Adding prompts to our GUI
We’re going to ask our user for the relevant method parameters. To do this, fire up QT Creator and open “your_plugin“_dialog_base.ui.
We’re going to start by adding seven line edit boxes. Rename these line edit boxes to “numBoids”, “minX”, “minY”, “maxX”, “maxY”, “tStep”, and “sFreq” respectively. Again, we’ll address the sFreq box shortly.
Next we’ll arrange these boxes and add labels. I put numBoids at the top, followed by tStep and sFreq. I put minX and minY adjacent to eachother above maxX and maxY.
Now we’ll add some text labels. I added a label next to each of the numBoids, tStep, and sFreq boxes and changed their text to “Number of Boids”, “Number of time steps”, and “Sampling frequency”, respectively. Additionally I added labels next to minX and minY and maxX and maxY labeled “Southwest corner” and “Northeast corner” respectively. I put labels above the x and y fields and changed their text to “X” and “Y”, respectively. Finally I added a label above the min and max boxes and changed the text to “study area bounds”.
Finally we can save the file.
Adding the text input to our code
Now we will handle the text from the text boxes we just added.
if result:
nB = int(self.dlg.numBoids.text())
nX = int(self.dlg.minX.text())
nY = int(self.dlg.minY.text())
mX = int(self.dlg.maxX.text())
mY = int(self.dlg.maxY.text())
tS = int(self.dlg.tStep.text())
sF = int(self.dlg.sFreq.text())
self.runBoids(nB, nX, nY, mX, mY, tS, sF)
Adding some functionality to our code
One thing we haven’t discussed yet is that our model simulates animal behaviours in discrete time. Animals react quickly to their environment, so quickly in fact that our 1 second canonical time makes our fish among the slowest-reacting animals that could exist. Ecologists usually don’t want to know the exact position of every animal every second, especially if we are limited in hard disk space, but if we increase the time interval our animals become very slow to react to their environment, potentially influencing their ecology. Since our boids track their own position independently of the data being saved, we can implement some code to only keep a fraction of the boids’ positions. This compromise can speed up our processing time significantly and save us some precious disk (and memory) space. We’ll use the modulus of our time step over our desired sampling frequency to get the points we actually save.
If we treat our time step as one second, which I am canonically doing here, setting sampleFreq to 60 will save the individuals’ locations every minute. 3600 will be every hour. 86400 will be every day. The model will still calculate new locations every second but we won’t save near as much data.
def runBoids(self, nB, xMin, yMin, xMax, yMax, timeSteps, sampleFreq):
. . .
for tS in range(timeSteps):
for boid in school:
boid.behave(school)
boid.update()
# Check if we want to store this location
if tS % sampleFreq == 0:
direction = np.arctan2(boid.velocity[0], boid.velocity[1]) * 57.29578
# add a feature
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(boid.position[0],boid.position[1])))
fet.setAttributes([str(boid.boidID), direction, tS])
prov.addFeatures([fet])
vLayer.updateExtents()
. . .
Telling QGIS users we are done
This code takes a while to run and it’s nice to let our user know that we finished successfully. To indicate a successful finish we can use the QGIS message bar.
self.iface.messageBar().pushMessage("Success", "Boids generated successfully.", level=Qgis.Success, duration=10)
Conclusion
If you’ve made it this far, congratulations on setting up an agent-based model in QGIS! I won’t recap all of the code here since including the plugin builder generated code makes it a bit much, but check out the project on my github if you want something quick to download.
Is this model useful? Not as written, unless you want to study the behaviours of a fish that only cares about the locations of conspecifics, doesn’t eat, has no habitat, and never changes depth. It could be a good starting point for other models though. YMMV. There is a lot more that can be done here: we could update our agent so that the parameters of our individual boids can be changed from the dialog, we could add some methods to our agent object so that it interfaces directly with the map layers and could possibly navigate or evaluate habitat from real data. Boids could vary individually, have random walk behaviours, interact with currents and chemical gradients, interact with agents from other species, exhibit diurnal behaviours, or even die! The modeling world is your oyster.