• Peter Kimberley, Gradata Systems Pty Ltd, Canberra, Australia Vietnam War Missing in Action (MIA) Support
• HIdenori Watanave, The University of Tokyo, Tokyo, Japan, Digital Archives of War and Disasters
• Ethan Berg, Agora World, Philadelphia, PA, USA, GeoForAll: Simplifying 3D Geospatial Metaverse Creation

So there we go! Working together collaboratively with Lutra Consulting on another great addition to QGIS 3D Functionality thanks to Cesium Ecosystem Grants. Stay tuned on our social channels to find out when it will be released in QGIS.

Annotation and Graphic Layer Improvements

Text Annotations along Curves

For those working with curved annotations, we’ve got you covered! SLYR now supports the conversion of text annotations along curves in QGIS. With this enhancement, you’ll get accurate conversion of any curved text and text-along-line annotations from MXD and APRX documents. This has been a long-requested feature which we can now introduce thanks to enhancements coming in QGIS 3.32.

ArcGIS Pro Graphics Layer Support

SLYR now supports the conversion of ArcGIS Pro graphics layers, converting all graphic elements to their QGIS “Annotation Layer” equivalents. If you’ve spent hours carefully designing cartographic markup on your maps, you can be sure that SLYR will allow you to re-use this work within QGIS!

Enhanced Page Layout Support

Improved QGIS to ArcGIS Pro Conversions

New tools

The latest SLYR release introduces some brand new tools for automating your conversion workflows:

All these exciting additions to SLYR are available today to SLYR license holders. If you’re after one-click, accurate conversion of projects from ESRI to QGIS, contact us to discuss your licensing needs.

As described on our SLYR page, we also provide some of the conversion tools for free use via the SLYR “Community Edition”. We’re proud to announce that we’ve just hit the next milestone in the Community Edition funding, and will now be releasing all of SLYR’s support for raster LYR files to the community edition! This complements the existing support for vector LYR files and ESRI style files available in the community edition. For more details on the differences between the licensed and community editions, see our product comparison.

Did you know that MovingPandas also supports local image coordinates? Indeed, it does.

In today’s post, we will explore how we can use this feature to analyze bicycle tracks extracted from video footage published by Michael Szell @mszll:

• Dataset: https://zenodo.org/record/7288616
• Data description: https://arxiv.org/abs/2211.01301

The bicycle trajectory coordinates are stored in two separate lists: xs_640x360 and ys640x360:

This format is kind of similar to the Kaggle Taxi dataset, we worked with in the previous post. However, to use the solution we implemented there, we need to combine the x and y coordinates into nice (x,y) tuples:

df['coordinates'] = df.apply(
    lambda row: list(zip(row['xs_640x360'], row['ys_640x360'])), axis=1)
df.drop(columns=['xs_640x360', 'ys_640x360'], inplace=True)

Afterwards, we can create the points and compute the proper timestamps from the frame numbers:

def compute_datetime(row):
    # some educated guessing going on here: the paper states that the video covers 2021-06-09 07:00-08:00
    d = datetime(2021,6,9,7,0,0) + (row['frame_in'] + row['running_number']) * timedelta(seconds=2)
    return d
def create_point(xy):
    try: 
        return Point(xy)
    except TypeError:  # when there are nan values in the input data
        return None
new_df = df.head().explode('coordinates')
new_df['geometry'] = new_df['coordinates'].apply(create_point)
new_df['running_number'] = new_df.groupby('id').cumcount()
new_df['datetime'] = new_df.apply(compute_datetime, axis=1)
new_df.drop(columns=['coordinates', 'frame_in', 'running_number'], inplace=True)
new_df

Once the points and timestamps are ready, we can create the MovingPandas TrajectoryCollection. Note how we explicitly state that there is no CRS for this dataset (crs=None):

trajs = mpd.TrajectoryCollection(
    gpd.GeoDataFrame(new_df), 
    traj_id_col='id',  t='datetime', crs=None)

Plotting trajectories with image coordinates

Similarly, to plot these trajectories, we should tell hvplot that it should not fetch any background map tiles (’tiles’:None) and that the coordinates are not geographic (‘geo’:False):

If you want to explore the full source code, you can find my Github fork with the Jupyter notebook at: https://github.com/anitagraser/desirelines/blob/main/mpd.ipynb

The repository also contains a camera image of the intersection, which we can use as a background for our trajectory plots:

bg_img = hv.RGB.load_image('img/intersection2.png', bounds=(0,0,640,360)) 

One important caveat is that speed will be calculated in pixels per second. So when we plot the bicycle speed, the segments closer to the camera will appear faster than the segments in the background:

To fix this issue, we would have to correct for the distortions of the camera lens and perspective. I’m sure that there is specialized software for this task but, for the purpose of this post, I’m going to grab the opportunity to finally test out the VectorBender plugin.

Georeferencing the trajectories using QGIS VectorBender plugin

Let’s load the five test trajectories and the camera image to QGIS. To make sure that they align properly, both are set to the same CRS and I’ve created the following basic world file for the camera image:

1
0
0
-1
0
360

Then we can use the VectorBender tools to georeference the trajectories by linking locations from the camera image to locations on aerial images. You can see the whole process in action here:

After around 15 minutes linking control points, VectorBender comes up with the following georeferenced trajectory result:

Not bad for a quick-and-dirty hack. Some points on the borders of the image could not be georeferenced since I wasn’t always able to identify suitable control points at the camera image borders. So it won’t be perfect but should improve speed estimates.

This post is part of a series. Read more about movement data in GIS.

Sorry, this entry is only available in Dutch.

Get your QGIS-ass to ‘s-Hertogenbosch!

Today’s post is a geeky deep dive into how to leverage DVC (not just) data version control to track QGIS geoprocessing workflows.

“Why is this great?” you may ask.

DVC tracks data, parameters, and code. If anything changes, we simply rerun the process and DVC will figure out which stages need to be recomputed and which can be skipped by re-using cached results.

This can lead to huge time savings compared to re-running the whole model

You can find the source code used in this post on my repo https://github.com/anitagraser/QGIS-resources/tree/dvc

I’m using DVC with the DVC plugin for VSCode but DVC can be used completely from the command line, if you prefer this appraoch.

Basically, what follows is a proof of concept: converting a QGIS Processing model to a DVC workflow. In the following screenshot, you can see the main stages

1. The QGIS model in the upper left corner
2. The Python script exported from the QGIS model builder in the lower left corner
3. The DVC stages in my dvc.yaml file in the upper right corner (And please ignore the hello world stage. It’s a left over from my first experiment)
4. The DVC DAG visualizing the sequence of stages. Looks similar to the QGIS model, doesn’t it ;-)

Besides the stage definitions in dvc.yaml, there’s a parameters file:

random-points:
  n: 10
buffer-points:
  size: 0.5

And, of course, the two stages, each as it’s own Python script.

First, random-points.py which reads the random-points.n parameter to create the desired number of points within the polygon defined in qgis3/data/test.geojson:

import dvc.api

from qgis.core import QgsVectorLayer
from processing.core.Processing import Processing
import processing

Processing.initialize()

params = dvc.api.params_show()
pts_n = params['random-points']['n']

input_vector = QgsVectorLayer("qgis3/data/test.geojson")
output_filename = "qgis3/output/random-points.geojson"

alg_params = {
    'INCLUDE_POLYGON_ATTRIBUTES': True,
    'INPUT': input_vector,
    'MAX_TRIES_PER_POINT': 10,
    'MIN_DISTANCE': 0,
    'MIN_DISTANCE_GLOBAL': 0,
    'POINTS_NUMBER': pts_n,
    'SEED': None,
    'OUTPUT': output_filename
}
processing.run('native:randompointsinpolygons', alg_params)

And second, buffer-points.py which reads the buffer-points.size parameter to buffer the previously generated points:

import dvc.api
import geopandas as gpd
import matplotlib.pyplot as plt

from qgis.core import QgsVectorLayer
from processing.core.Processing import Processing
import processing

Processing.initialize()

params = dvc.api.params_show()
buffer_size = params['buffer-points']['size']

input_vector = QgsVectorLayer("qgis3/output/random-points.geojson")
output_filename = "qgis3/output/buffered-points.geojson"

alg_params = {
    'DISSOLVE': False,
    'DISTANCE': buffer_size,
    'END_CAP_STYLE': 0,  # Round
    'INPUT': input_vector,
    'JOIN_STYLE': 0,  # Round
    'MITER_LIMIT': 2,
    'SEGMENTS': 5,
    'OUTPUT': output_filename
}
processing.run('native:buffer', alg_params)

gdf = gpd.read_file(output_filename)
gdf.plot()

plt.savefig('qgis3/output/buffered-points.png')

With these things in place, we can use dvc to run the workflow, either from within VSCode or from the command line. Here, you can see the workflow (and how dvc skips stages and fetches results from cache) in action:

If you try it out yourself, let me know what you think.

The QGIS conda packages have been around for a while. One of their use cases, for example, is to allow Linux users to easily install multiple versions of QGIS.

Similarly, we’ve seen posts on using PyQGIS in Jupyter notebooks. However, I find the setup with *.bat files rather tricky.

This post presents a way to set up a conda environment with QGIS that is ready to be used in Jupyter notebooks.

The first steps are to create a new environment and install QGIS. I use mamba for the installation step because it is faster than conda but you can use conda as well:

(base) PS C:\Users\anita> conda create -n qgis python=3.9
(base) PS C:\Users\anita> conda activate qgis
(qgis) PS C:\Users\anita> mamba install -c conda-forge qgis=3.28.2 
(qgis) PS C:\Users\anita> qgis

If we now try to import the qgis module in Python, we get an error:

(qgis) PS C:\Users\anita> python
Python 3.9.15 | packaged by conda-forge | (main, Nov 22 2022, 08:41:22) [MSC v.1929 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'qgis'

To fix this error, we need to get the paths from the Python console inside QGIS:

import sys
sys.path
['H:/miniconda3/envs/qgis/Library/./python', 'C:/Users/anita/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', ... ]

This list of paths can be configured as the defaults for our qgis environment using conda develop:

(qgis) PS C:\Users\anita> conda activate base
(base) PS C:\Users\anita> mamba install conda-build -c conda-forge
(base) PS C:\Users\anita> conda develop -n qgis [list of paths from qgis python console] 

With this setup, the import should now work without errors:

(base) PS C:\Users\anita> conda activate qgis
(qgis) PS C:\Users\anita> python
Python 3.9.15 | packaged by conda-forge | (main, Nov 22 2022, 08:41:22) [MSC v.1929 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis

The example Jupyter notebook covers running a QGIS Processing algorithm and visualizing the results in the notebook using GeoPandas:

Head over to Github to find the full instructions: https://github.com/anitagraser/QGIS-resources/blob/master/qgis3/notebooks/hello-world.ipynb

Thanks to the sponsoring of the Swiss QGIS User Group, starting from QGIS 3.26 is it possible to override field names in the layer export dialog. Previous to that, QGIS would always export with the technical names from the database, whereas now it’s possible to override with the alias defined in QGIS or any custom name. One use for this in Switzerland — a highly polyglot country — is an export with translated names.

This is done via an additional column “Export name”. For convenience we also added a tri-state checkbox to toggle export names to their alias defined in the layer configuration or back to the field name. If a name is changed by hand the checkbox shows a mixed state.

QField is a community-driven open-source project. It is free to share, use and modify and it will stay like that. The very essence of a community is to help and support each other. And that’s where YOU come into play. To make it work we need your support!

For those who don’t know much about the concept of open source projects, a bit of background. Investing in open-source projects is a technical and ethical decision for OPENGIS.ch. Open source is a technological advantage, as we receive input from many developers worldwide who are motivated to work out the best possible software. It prevents our customers from vendor lock-in and allows complete ownership and control of the developed software. And finally, not only financially independent businesses and people should benefit from professional software but also those who might not have the financial means to pay for features, and licences. 

You are not a developer, but you still like to use QField and support it? Good news. You don’t have to be a developer to use, contribute or recommend the app. There are plenty of things that need to be done to help QField to remain the powerful software it is right now and become even better. Here are a few suggestions on how you can give something back.

1. Review the app ★★★★★ in google’s play store or apple’s app store. 
2. Let the world know about it! It doesn’t matter if you’re on Twitter, LinkedIn, Instagram or any other social media platform. Show and tell about where QField helped you. We appreciate every post and we promise to like, share and comment.
3. Write about your experience and please let us know. Be it in your blog or as a new success story. Insights into field projects are extremely valuable. It helps us to make the app even more efficient for your work, and it helps others to understand the range of applications for QField.
4. Register for a paid QFieldCloud account. QFieldCloud allows to synchronize and merge the data collected in QField. QFieldCloud is hosted by the makers of QField and by getting an account you help QField too.

Do you want to do something that is more hands-on and directly linked to the app? No problem. 

4. Help with the documentation. You can document features, or improve the documentation in English. Read the how-to guide to get started.
5. And if you are multilingual you might consider translating the documentation or the app in your language.
6. Become a beta tester and be the first to report a bug! When something doesn’t work properly it might be a bug. The quicker we know about it, the faster it can be resolved.
7. You can ask and answer questions on gis.stackexchange and help others on the user discussions platform.
8. If you are a developer and you want to get involved in QField development, please refer to the individual documentation for QField, QFieldCloud and QFieldSync.

And now finally for those of you who have the financial means, you can either sponsor a feature or subscribe to one of the monthly sponsorships. By doing so you help get freshly baked QField versions straight to everyone’s devices.

Nothing in it for you? In that case, just drop by to say thank you or have a hot or cold beverage with us next time you meet OPENGIS.ch at a conference and you might make our day!
Want to know more about the idea of community-driven open-source projects and the QGIS project in particular? Check out Nyall Dawson’s blog post about how to effectively get things done in open source!

This blog post is about QGIS relations and how they are edited in the attribute form with widgets in general, as well as some plugins that override the relations editor widget to improve usability and solve specific use cases. The start is quite basic. If you are already a relation hero, then jump directly to the plugins.

QGIS Relations in General

Let’s have a look at a simple example data model. We have four entities: Building, Apartment, Address and Owner. In UML it looks like this:

• A building can have none or multiple apartments, but an apartment must to be related to a building. This black box on the left describes the relation strength as a composition. An apartment cannot exist without a building. When a building is demolished, all apartments of it are demolished as well.
• An apartment needs to be owned by at least one owner. An owner can own none or more apartments. This is a many-to-many relation and this means, it will be normalized by adding a linking (join) table in between.
• A building can have an address (only one – no multiple entrances in this example). An address can refer to one building.
  Why not making one single table on a one-to-one relation? To ensure their existence independently: When a building is demolished, the address should persist until the new building is constructed.

Creating Relations in QGIS

In QGIS we have now five layers. The four entities and the linking table called “Apartment_Owner”.

Open Project > Properties… > Relations

With Discover Relations the possible relations are detected from the existing layers according to their foreign keys in the database. In this example no CASCADE is defined in the database what means that the relations strength is always “Association”.

Where would “Composition” make sense?

Of course in the relation “Apartment” to “Building”, to ensure that when a feature of “Building” is deleted, the children (“Apartment”) are deleted as well, because they cannot exist without a building. Also a duplication of a feature of “Building” would duplicate the children (“Apartment”) as well.

But as well on the linking (join) table “Apartment_Owner” and its relation to “Apartment” and “Owner” a composition would make sense. Because when a feature of “Apartment” or “Owner” is deleted, the entry in the linking table should be deleted as well. Because this connection does not exist anymore and otherwise this would lead to orphan entries in the linking table.

Walk through the widgets

To demonstrate the relation widgets Relation Editor, Relation Reference and Value Relation we make a walk through the digitizing process.

Relation Editor

First we create a “Building” and call it “Garden Tower”. Then we add some “Apartments”.

The “Apartments” are created in the widget called Relation Editor. This shows us a list (similar to the QGIS Attribute Table) of all children (“Apartment”) referencing to this “Building”. We have here activated the possibilities to add, delete and duplicate child-features.

In the widget settings (Right-click on the layer > Properties… > Attribute Form) we see that there are other possibilities to link and unlink child-features as well as zoom to the current child-feature (what only would make sense when they have a geometry).

As well we can set here the cardinality. This will become interesting when we go to the “Owner” to “Apartment” relation. But let’s first have a look at the opposite of what we just did.

Relation Reference

When we open now a feature of “Apartment”, we see that we have a drop down to select the “Building” to reference to.

On the right of this drop down we can see some buttons. Those are for the following functionalities (from left to right):

• Open the form of the current parent feature (in our case the “Building” feature called “Garden Tower”)
• Add a new feature on the parent layer (in our case “Building”)
• Highlight the parent layer (in our case “Building”) on the map
• Select the parent feature (in our case “Building”) on the map to reference it

In the settings (Right-click on the layer > Properties… > Attribute Form) we see that we choose the configured relation to connect the child (“Apartment”) to the parent (“Building”). This won’t be needed with the widget Value Relation.

Value Relation

The Value Relation does not require a relation at all. We simply choose the “parent” layer (“Building”) its primary key as the key (“t_id”) and a descriptive field as the value (“Description”).

The result shows us a drop down as well to select the parent.

It is much easier to configure, but you can see the limitations. There are no such functionalities to control the parent feature like add, identify on map etc. As well you need to be careful to fulfill the foreign key constraint (you have to choose the correct field to link with). All this is given, when you build a Relation Reference on an existing relation.

Many-to-Many Relations

Now we link some “Owner” to our “Apartment”. We could create new ones like we did it for the “Apartment” in “Building” or we can link existing ones. For linking we choose the yellow link-button on the top of the Relation Editor.

This dialog looks similar to the Relation Editor widget. You have just to select the “Owner” you want to link to the “Apartment” by checking the yellow box. It’s a very powerful tool, but people are often confused about the load of functionality here and the selection that can be difficult to get used to (yellow boxes vs. blue index selection). For this case we extended the Relation Editor widget with a plugin.

Anyway after that we linked our features of the layer “Owner”.

Have you seen the linking table in between? Well, me neither. It’s completely invisible for the end user. This because of the cardinality setting I mentioned already. When we choose the linked table “Owner” instead of “Many to one relation”, then we can create and link the other parent (“Owner”) directly.

One-to-One Relation

A one-to-one relation like we have here between “Building” and “Address” is created in the database more or less like a normal one-to-many relation. This means one of the tables (in our case “Address”) has a foreign key pointing to the parent table (“Building”). There are tricks to fulfill the one-to-one maximum cardinality (like e.g. by setting a UNIQUE constraint on this foreign key column) but still in the QGIS user interface it looks like a one-to-many relation. It’s displayed in a normal Relation Editor widget.

Solutions could be so called “Joins”. Go to the settings (Right-click on the layer > Properties… > Joins)

Here you can join a layer of your choice and add the fields of this other layer (in our case “Address”) to your current feature form (of “Building”). So it appears to the user that it’s the same table containing fields of “Building” and “Address”.

Negative point about those joins are, that they are fault prone. You have to be careful with default values (e.g. on primary keys) of the joined layer. You cannot expect a fully reliable feature form like you have it in the Relation Editor. Here as well, we extended the Relation Editor widget with a plugin.

Plugins for Relation Editor Widgets

Since QGIS 3.18 the base class of the Relation Editor Widgets became abstract, what opened the possibility to use it in PyQGIS and derive it to super nice widgets handling specific use cases and improving the usability.

Linking Relation Editor Widget

As mentioned before, the QGIS stock dialog to link children is full of features but it can be overwhelming and difficult to use. Mostly because of the two selection possibilities in the list. A blue selection is for the currently displayed feature, and a yellow checkbox selection is for the features to be actually linked.

In collaboration with the Model Baker Group we wanted to improve the situation. But as we where unsure how the end solution should look like, so we decided to experiment in a plugin. The result is a link manger dialog, in which features can be linked and unlinked by moving them left and right. The effective link is created or destroyed when the dialog is accepted.

Find more information on the repository https://github.com/opengisch/qgis-linking-relation-editor

Ordered Relation Editor Widget

Sometimes the order of the children play a role on the project, and you want to have them displayed following that. For that there is the Ordered Relation Editor Widget. You can configure a field in the children to be used to order them. In the given example the field Floor was used to order Apartments. Reordering the fields by Drag&Drop would change the value of the configured field. Display name and optionally a path to an icon to be shown on the list can be configured by expression in the Attribute Form tab in the layer properties (Right-click on the layer > Properties… > Attribute Form).

Document Management System Widgets

Often in QGIS projects there is the need to deal with external documents. This could be for example pictures, documentations or reports about some features. To support that we added two new tables in the project:

• Documents each document is represented by a row in this table. The table has following fields:
  • id
  • path is the filename of the document.
• DocumentsFeatures this is a linking (join) table and permits to link a document with one or more features in more layers. The table has following fields:
  • id
  • document_id id of the document.
  • feature_id id of the feature.
  • feature_layer layer of the feature.

Thanks to a QGIS feature named Polymorphic Relations we can link a document with features of multiple layers. The polymorphic relation can evaluate an expression to decide in which table will be the feature to link. Here a screenshot of the relation configuration:

After this configuration in the layers “Apartment” and “Building” it will be possible to link children from the “Documents” table. The document management plugin provides two widgets to simplify the handling of the relation. In the feature side widget the documents are displayed as a grid or list. If possible a preview of the contend is shown and you can add new documents via Drag&Drop from the system file manager. Double-click on a document will open it in the default system viewer.

The second widget is meant to be used in the Feature Form of the “Documents” table, and it permits to handy see, for each document, with which feature from which layer it is linked.

Find more information on the repository https://github.com/opengisch/qgis-document-management-system-plugin

That’s it

Well then. We hope that all the beginners reading this article received some light on QGIS Relations and all the advanced user some inspiration on the immense possibilities you have with QGIS ?

The QGIS 3.28 release is an extremely exciting release for all users who work in mixed software workplaces, or who need to work alongside users of ESRI software. In this post we’ll be giving an overview of all the new tools and features introduced in 3.28 which together result in a dramatic improvement in the workflows and capabilities in working with ESRI based formats and services. Read on for the full details…!

Before we begin, we’d like to credit the following organisations for helping fund these developments in QGIS 3.28:

• Naturstyrelsen, Denmark
• Provincie Gelderland, Netherlands
• Uppsala Universitet, Department of Archaeology and Ancient History
• Gemeente Amsterdam
• Provincie Zuid-Holland, Netherlands

FileGeodatabase (GDB) related improvements

The headline item here is that QGIS 3.28 introduces support for editing, managing and creating ESRI FileGeodatabases out of the box! While older QGIS releases offered some limited support for editing FileGeodatabase layers, this required the manual installation of a closed source ESRI SDK driver… which unfortunately resulted in other regressions in working with FileGeodatabases (such as poor layer loading speed and random crashes). Now, thanks to an incredible reverse engineering effort by the GDAL team, the open-source driver for FileGeodatabases offers full support for editing these datasets! This means all QGIS users have out-of-the-box access to a fully functional, high-performance read AND write GDB driver, no further action or trade-offs required.

Operations supported by the GDAL open source driver include:

• Editing existing features, with full support for editing attributes and curved, 3D and measure-value geometries
• Creating new features
• Deleting features
• Creating, adding and modifying attributes in an existing layer
• Full support for reading and updating spatial indexes
• Creating new indexes on attributes
• “Repacking” layers, to reduce their size and improve performance
• Creating new layers in an existing FileGeodatabase
• Removing layers from FileGeodatabases
• Creating completely new, empty FileGeodatabases
• Creating and managing field domains

On the QGIS side, the improvements to the GDAL driver meant that we could easily expose feature editing support for FileGeodatabase layers for all QGIS users. While this is a huge step forward, especially for users in mixed software workplaces, we weren’t happy to rest there when we  had the opportunity to further improve GDB support within QGIS!

So in QGIS 3.28 we also introduced the following new functionality when working with FileGeodatabases:

FileGeodatabase management tools

QGIS 3.28 introduces a whole range of GUI based tools for managing FileGeodatabases. To create a brand new FileGeodatabase, you can now right click on a directory from the QGIS Browser panel and select New – ESRI FileGeodatabase:

After creating your new database, a right click on its entry will show a bunch of available options for managing the database. These include options for creating new tables, running arbitrary SQL commands, and database-level operations such as compacting the database:

You’re also able to directly import existing data into a FileGeodatabase by simply dragging and dropping layers onto the database!

Expanding out the GDB item will show a list of layers present in the database, and present options for managing the fields in those layers. Alongside field creation, you can also remove and rename existing fields.

Field domain handling

QGIS 3.28 also introduces a range of GUI tools for working with field domains inside FileGeodatabases. (GeoPackage users also share in the love here — these same tools are all available for working with field domains inside this standard format too!) Just right click on an existing FileGeodatabase (or GeoPackage) and select the “New Field Domain” option. Depending on the database format, you’ll be presented with a list of matching field domain types:

Once again, you’ll be guided through a user-friendly dialog allowing you to create your desired field domain!

After field domains have been created, they can be assigned to fields in the database by right-clicking on the field name and selecting “Set Field Domain”:

Field domains can also be viewed and managed by expanding out the “Field domains” option for each database.

Relationship discovery

Another exciting addition in QGIS 3.28 (and the underlying GDAL 3.6 release) is support for discovering database relationships in FileGeodatabases! (Once again, GeoPackage users also benefit from this, as we’ve implemented full support for GeoPackage relationships via the “Related Tables Extension“).

Expanding out a database containing any relationships will show a list of all discovered relationships:

(You can view the full description and details for any of these relationships by opening the QGIS Browser “Properties” panel).

Whenever QGIS 3.28 discovers relationships in the database, these related tables will automatically be added to your project whenever any of the layers which participate in the relationship are opened. This means that users get the full experience as designed for these databases without any manual configuration, and the relationships will “just work”!

Dataset Grouping

Lastly, we’ve improved the way layers from FileGeodatabases are shown in QGIS, so that layers are now grouped according to their original dataset groupings from the database structure:

Edit ArcGIS Online / Feature Service layers

While QGIS has had read-only support for viewing and working with the data in ArcGIS Online (AGOL) vector layers and ArcGIS Server “feature service” layers for many years, we’ve added support for editing these layers in QGIS 3.28. This allows you to take advantage of all of QGIS’ easy to use, powerful editing tools and directly edit the content in these layers from within your QGIS projects! You can freely create new features, delete features, and modify the shape and attributes of existing features (assuming that your user account on the ArcGIS service has these edit permissions granted, of course). This is an exciting addition for anyone who has to work often with content in ArcGIS services, and would prefer to directly manipulate these layers from within QGIS instead of the limited editing tools available on the AGOL/Portal platforms themselves.

This new functionality will be available immediately to users upon upgrading to QGIS 3.28 — any users who have been granted edit capabilities for the layers will see that the QGIS edit tools are all enabled and ready for use without any further configuration on the QGIS client side.

Filtering Feature Service layers

We’ve also had the opportunity to introduce filter/query support for Feature Service layers in QGIS 3.28. This is a huge performance improvement for users who need to work with a subset of a features from a large Feature Service layer. Unfortunately, due to the nature of the Feature Service protocol, these layers can often be slow to load and navigate on a client side. By setting a SQL filter to limit the features retrieved from the service the performance can be dramatically increased, as only matching features will ever be requested from the backend server. You can use any SQL query which conforms to the subset of SQL understood by ArcGIS servers (see the Feature Service documentation for examples of supported SQL queries).

What’s next?

While QGIS 3.28 is an extremely exciting release for any users who need to work alongside ESRI software, we aren’t content to rest here! The exciting news is that in QGIS 3.30 we’ll be introducing a GUI driven approach allowing users to create new relationships in their FileGeodatabase (and GeoPackage!) databases.

At North Road we’re always continuing to improve the cross-vendor experience for both ESRI and open-source users through our continued work on the QGIS desktop application and our SLYR conversion suite. If you’d like to chat to us about how we can help your workplace transition from a fully ESRI stack to a mixed or fully open-source stack, just contact us to discuss your needs.

In the previous post, we — creatively ;-) — used MobilityDB to visualize stationary IOT sensor measurements.

This post covers the more obvious use case of visualizing trajectories. Thus bringing together the MobilityDB trajectories created in Detecting close encounters using MobilityDB 1.0 and visualization using Temporal Controller.

Like in the previous post, the valueAtTimestamp function does the heavy lifting. This time, we also apply it to the geometry time series column called trip:

SELECT mmsi,
    valueAtTimestamp(trip, '2017-05-07 08:55:40') geom,
    valueAtTimestamp(SOG, '2017-05-07 08:55:40') SOG
FROM "public"."ships"

Using this SQL query, we again set up a — not yet Temporal Controller-controlled — QueryLayer.

To configure Temporal Controller to update the timestamp in our SQL query, we again need to run the Python script from the previous post.

With this done, we are all set up to animate and explore the movement patterns in our dataset:

This post is part of a series. Read more about movement data in GIS.

Our SLYR tool is the complete solution for full compatibility between ArcMap, ArcGIS Pro and QGIS. It offers a powerful suite of conversion tools for opening ESRI projects, styles and other documents directly within QGIS, and for exporting QGIS documents for use in ESRI software.

A lot has changed since our last SLYR product update post, and we’ve tons of very exciting improvements and news to share with you all! In this update we’ll explore some of the new tools we’ve added to SLYR, and discuss how these tools have drastically improved the capacity for users to migrate projects from the ESRI world to the open-source world (and vice versa).

ArcGIS Pro support

The headline item here is that SLYR now offers a powerful set of tools for working with the newer ArcGIS Pro document formats. Previously, SLYR offered support for the older ArcMap document types only (such as MXD, MXT, LYR, and PMF formats). Current SLYR versions now include tools for:

Directly opening ArcGIS Pro .lyrx files within QGIS

LYRX files can be dragged and dropped directly onto a QGIS window to add the layer to the current project. All the layer’s original styling and other properties will be automatically converted across, so the resultant layer will be an extremely close match to the original ArcGIS Pro layer! SLYR supports vector layers, raster layers, TIN layers, point cloud layers and vector tile layers. We take great pride in just how close the conversion results are to how these layers appear in ArcGIS Pro… in most cases you’ll find the results are nearly pixel perfect!

In addition to drag-and-drop import support, SLYR also adds support for showing .lyrx files directly in the integrated file browser, and also adds tools to the QGIS Processing Toolbox so that users can execute bulk conversion operations, or include document conversion in their models or custom scripts.

ArcGIS Pro map (mapx) and project (aprx) conversion

Alongside the LYRX support, we’ve also added support for the ArcGIS Pro .mapx and .aprx formats. Just like our existing .mxd conversion, you can now easily convert entire ArcGIS Pro maps for direct use within QGIS! SLYR supports both the older ArcGIS Pro 2.x project format and the newer 3.x formats.

Export from QGIS to ArcGIS Pro!

Yes, you read that correctly… SLYR now allows you to export QGIS documents into ArcGIS Pro formats! This is an extremely exciting development… for the first time ever QGIS users now have the capacity to export their work into formats which can be supplied directly to ESRI users. Current SLYR versions support conversion of map layers to .lyrx format, and exporting entire projects to the .mapx format. (We’ll be introducing support for direct QGIS to .aprx exports later this year.)

We’re so happy to finally provide an option for QGIS users to work alongside ArcGIS Pro users. This has long been a pain point for many organisations, and has even caused organisations to be ineligible to tender for jobs which they are otherwise fully qualified to do (when tenders require provision of data and maps in ArcGIS compatible formats).

ArcGIS Pro .stylx support

Alongside the other ArcGIS Pro documents, SLYR now has comprehensive support for reading and writing ArcGIS Pro .stylx databases. We’ve dedicated a ton of resources in ensuring that the conversion results (both from ArcGIS Pro to QGIS and from QGIS to ArcGIS Pro) are top-notch, and we even handle advanced ArcGIS Pro symbology options like symbol effects!

Take a look below how even very advanced ArcGIS Pro style libraries convert beautifully to QGIS symbol libraries:

ArcMap Improvements

While we’ve been focusing heavily on the newer ArcGIS Pro formats, we’ve also improved our support for the older ArcMap documents. In particular, SLYR now offers more options for converting ArcMap annotation layers and annotation classes to QGIS supported formats. Users can now convert Annotation layers and classes directly over to QGIS annotation layer or alternatively annotation classes can be converted over to the OGC standard GeoPackage format. When exporting annotation classes to GeoPackage the output database is automatically setup with default styling rules, so that the result can be opened directly in QGIS and will be immediately visualised to match the original annotation class.

Coming soon…

While all the above improvements are already available for all SLYR license holders, we’ve got many further improvements heading your way soon! For example, before the end of 2022 we’ll be releasing another large SLYR update which will introduce support for exporting QGIS projects directly to ArcGIS Pro .aprx documents. We’ve also got many enhancements planned which will further improve the quality of the converted documents. Keep an eye on this blog and our social media channels for more details as they are available…

You can read more about our SLYR tool at the product page, or contact us today to discuss licensing options for your organisation.

Thanks to the recent popularity of the “30 Day Map Challenge“, the month of November has become synonymous with beautiful maps and cartography. During this November we’ll be sharing a bunch of tips and tricks which utilise some advanced QGIS functionality to help create beautiful maps.

One technique which can dramatically improve the appearance of maps is to swap out rectangular inset maps for more organic shapes, such as circles or ovals.

Back in 2020, we had the opportunity to add support for directly creating circular insets in QGIS Print Layouts (thanks to sponsorship from the City of Canning, Australia!). While this functionality makes it easy to create non-rectangular inset maps the steps, many QGIS users may not be aware that this is possible, so we wanted to highlight this functionality for our first 30 Day Map Challenge post.

Let’s kick things off with an example map. We’ve shown below an extract from the 2032 Brisbane Olympic Bid that some of the North Road team helped create (on behalf of SMEC for EKS). This map is designed to highlight potential venues around South East Queensland and the travel options between these regions:

Circles featured heavily in previous Olympic bid maps (such as Budapest) where we took our inspiration from. This may, or may not, play a part in using the language of the target map audience – think Olympic rings!

Budapest Olympics 2024 Masterplan

Step by Step Guide to Creating a Circle Inset

Firstly, prepare a print layout with both a main map and an inset map. Make sure that your inset map is large enough to cover your circular shape:

From the Print Layout toolbar, click on the Add Shape button and then select Add Ellipse:

Draw the ellipse over the middle of your inset map (hint: holding down Shift while drawing the ellipse will force it to a circular shape!). If you didn’t manage to create an exact circle then you can manually specify the width and height in the shape item’s properties. For this one, we went with a 50mm x 50mm circle:

Next, select the Inset Map item and in its Item Properties click on the Clipping Settings button:

In the Clipping Settings, scroll down to the second section and tick the Clip to Item box and select your Ellipse item from the list. (If you have labels shown in your inset map you may also want to check the “force labels inside clipping shape” option to force these labels inside the circle. If you don’t check this option then labels will be allowed to overflow outside of the circle shape.)

Your inset map will now be bound to the ellipse!

Here’s a bit more magic you could add to this map – in the Main Map’s properties, click on Overviews and set create one for the Inset map – it will nicely show the visible circular area and not the rectangle!

Bonus Points: Circular Title Text!

For advanced users, we’ve another fun tip…and when we say fun, we mean ‘let’s play with radians’! Here we’re going to create some title text and a wedged background which curves around the outside of our circular inset. This takes some fiddly playing around, but the end result can be visually striking! Here we’re going to push the QGIS print layout “HTML” item to create some advanced graphics, so some HTML and CSS coding experience is advantageous. (An alternative approach would be to use a vector illustration application like Inkscape, and add your title and circular background as an SVG item in the print layout).

We’ll start by creating some curved circular text:

First, add a “HTML frame” to your print layout:

HTML frames allow placement of dynamic content in your layouts, which can use HTML, CSS and JavaScript to create graphical components.

In the HTML item’s “source” box, add the following code:

<svg height="300" width="350">
        <defs>
            <clipPath id="circleView">
                <circle id="curve" cx="183" cy="156" r="25" fill="transparent" />
            </clipPath>
        </defs>
        <path id="forText" d="M 28,150, C 25,50, 180,-32,290,130" stroke="" fill="none"/>
            <text x="0" y="35" width="100">
                <textpath xlink:href="#forText">
                    <tspan font-weight="bold" fill="black">Place text here</tspan>
                </textpath>
            </text>
             <style>
    <![CDATA[
      text{
        dominant-baseline: hanging;
        font: 20px Arial;
      }
    ]]>
  </style>
</svg>

Now, let’s add in a background to bring more focus onto the title!

To add in the background, create another HTML item. We’ll again create the arc shape using an SVG element, so add the following code into the item’s source box:

<svg width="750" height="750" xmlns="http://www.w3.org/2000/svg">
  <path d="M 90 70
           A 56 56, 0, 0, 0, 133 140
           L 150 90 Z" fill="#414042" transform=" scale(2.1) rotate(68 150 150) " />/>
</svg>

(You can read more about SVG  curves and arcs paths over at MDN)

So there we go! These two techniques can help push your QGIS map creations further and make it easier to create beautiful cartography directly in QGIS itself. If you found these tips useful, keep an eye on this blog as we post more tips and tricks over the month of November. And don’t forget to follow the 30 day Map Challenge for a smorgasbord of absolutely stunning maps.

Amongst all the processing algorithms already available in QGIS, sometimes the one thing you need is missing. 

This happened not a long time ago, when we were asked to find a way to continuously visualise traffic on the Swiss motorway network (polylines) using frequently measured traffic volumes from discrete measurement stations (points) alongside the motorways. In order to keep working with the existing polylines, and be able to attribute more than one value of traffic to each feature, we chose to work with the M-values. M-values are a per-vertex attribute like X, Y or Z coordinates. They contain a measure value, which typically represents time or distance. But they can hold any numeric value.

In our example, traffic measurement values are provided on a separate point layer and should be attributed to the M-value of the nearest vertex of the motorway polylines. Of course, the motorway features should be of type LineStringM in order to hold an M-value. We then should interpolate the M-values for each feature over all vertices in order to get continuous values along the line (i.e. a value on every vertex). This last part is not yet existing as a processing algorithm in QGIS.

This article describes how to write a feature-based processing algorithm based on the example of M-value interpolation along LineStrings.

Feature-based processing algorithm

The pyqgis class QgsProcessingFeatureBasedAlgorithm is described as follows: “An abstract QgsProcessingAlgorithm base class for processing algorithms which operates “feature-by-feature”.  

Feature based algorithms are algorithms which operate on individual features in isolation. These are algorithms where one feature is output for each input feature, and the output feature result for each input feature is not dependent on any other features present in the source. […]

Using QgsProcessingFeatureBasedAlgorithm as the base class for feature based algorithms allows shortcutting much of the common algorithm code for handling iterating over sources and pushing features to output sinks. It also allows the algorithm execution to be optimised in future (for instance allowing automatic multi-thread processing of the algorithm, or use of the algorithm in “chains”, avoiding the need for temporary outputs in multi-step models).”

In other words, when connecting several processing algorithms one after the other – e.g. with the graphical modeller – these feature-based processing algorithms can easily be used to fill in the missing bits. 

Compared to the standard QgsProcessingAlgorithm the feature-based class implicitly iterates over each feature when executing and avoids writing wordy loops explicitly fetching and applying the algorithm to each feature. 

Just like for the QgsProcessingAlgorithm (a template can be found in the Processing Toolbar > Scripts > Create New Script from Template), there is quite some boilerplate code in the QgsProcessingFeatureBasedAlgorithm. The first part is identical to any QgsProcessingAlgorithm.

After the description of the algorithm (name, group, short help, etc.), the algorithm is initialised with def initAlgorithm, defining input and output. 

In our M-value example:

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr('Input layer'),
                [QgsProcessing.TypeVectorAnyGeometry]
            )
        )
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr('Output layer')
            )
        )

While in a regular processing algorithm now follows def processAlgorithm(self, parameters, context, feedback), in a feature-based algorithm we use def processFeature(self, feature, context, feedback). This implies applying the code in this block to each feature of the input layer. 

Interpolating M-values

This actual processing part can be copied and added almost 1:1 from any other independent python script, there is little specific syntax to make it a processing algorithm. Only the first line below really.

In our M-value example:

    def processFeature(self, feature, context, feedback):
        
        try:
            geom = feature.geometry()
            line = geom.constGet()
            vertex_iterator = QgsVertexIterator(line)
            vertex_m = []

            # Iterate over all vertices of the feature and extract M-value

            while vertex_iterator.hasNext():
                vertex = vertex_iterator.next()
                vertex_m.append(vertex.m())

            # Extract length of segments between vertices

            vertices_indices = range(len(vertex_m))
            length_segments = [sqrt(QgsPointXY(line[i]).sqrDist(QgsPointXY(line[j]))) 
                for i,j in itertools.combinations(vertices_indices, 2) 
                if (j - i) == 1]

            # Get all non-zero M-value indices as an array, where interpolations 
              have to start

            vertex_si = np.nonzero(vertex_m)[0]
            
            m_interpolated = np.copy(vertex_m)

            # Interpolate between all non-zero M-values - take segment lengths between 
              vertices into account

            for i in range(len(vertex_si)-1):
                first_nonzero = vertex_m[vertex_si[i]]
                next_nonzero = vertex_m[vertex_si[i+1]]
                accum_dist = itertools.accumulate(length_segments[vertex_si[i]
                                                                  :vertex_si[i+1]])
                sum_seg = sum(length_segments[vertex_si[i]:vertex_si[i+1]])
                interp_m = [round(((dist/sum_seg)*(next_nonzero-first_nonzero)) + 
                            first_nonzero,0) for dist in accum_dist]
                m_interpolated[vertex_si[i]:vertex_si[i+1]] = interp_m

            # Copy feature geometry and set interpolated M-values, 
              attribute new geometry to feature

            geom_new = QgsLineString(geom.constGet())
            
            for j in range(len(m_interpolated)):
                geom_new.setMAt(j,m_interpolated[j])
                
            attrs = feature.attributes()
            
            feat_new = QgsFeature()
            feat_new.setAttributes(attrs)
            feat_new.setGeometry(geom_new)

        except Exception:
            s = traceback.format_exc()
            feedback.pushInfo(s)
            self.num_bad += 1
            return []
        
        return [feat_new]

In our example, we get the feature’s geometry, iterate over all its vertices (using the QgsVertexIterator) and extract the M-values as an array. This allows us to assign interpolated values where we don’t have M-values available. Such missing values are initially set to a value of 0 (zero).

We also extract the length of the segments between the vertices. By gathering the indices of the non-zero M-values of the array, we can then interpolate between all non-zero M-values, considering the length that separates the zero-value vertex from the first and the next non-zero vertex.

For the iterations over the vertices to extract the length of the segments between them as well as for the actual interpolation between all non-zero M-value vertices we use the library itertools. This library provides different iterator building blocks that come in quite handy for our use case. 

Finally, we create a new geometry by copying the one which is being processed and setting the M-values to the newly interpolated ones.

And that’s all there is really!

Styling the result via M-values

The last step is styling our output layer in QGIS, based on the M-values (our traffic M-values are categorised from 1 [a lot of traffic -> dark red] to 6 [no traffic -> light green]). This can be achieved by using a Single Symbol symbology with a Marker Line type “on every vertex”. As a marker type, we use a simple round point. Stroke style is “no pen” and Stroke fill is based on an expression:

with_variable(

'm_value', m(point_n($geometry, @geometry_point_num)),

	CASE WHEN @m_value = 6
		THEN color_rgb(140, 255, 159)

		WHEN @m_value = 5
			THEN color_rgb(244, 252, 0)

		WHEN @m_value = 4
			THEN color_rgb(252, 176, 0)

		WHEN @m_value = 3
			THEN color_rgb(252, 134, 0)

		WHEN @m_value = 2
			THEN color_rgb(252, 29, 0)

		WHEN @m_value = 1
			THEN color_rgb(140, 255, 159)

		ELSE
			color_hsla(0,100,100,0)

	END
)

And voilà! Wherever we have enough measurements on one line feature, we get our motorway network continuously coloured according to the measured traffic volume.

One disclaimer at the end: We get this seemingly continuous styling only because of the combination of our “complex” polylines (containing many vertices) and the zoomed-out view of the motorway network. Because really, we’re styling many points and not directly the line itself. But in our case, this is working very well.

If you’d like to make your custom processing algorithm available through the processing toolbox in your QGIS, just put your script in the folder containing the files related to your user profile:

profiles > default > processing > scripts 

You can directly access this folder by clicking on Settings > User Profiles > Open Active Profile Folder in the QGIS menu.

That way, it’s also available for integration in the graphical modeller.

You can download the above-mentioned processing scripts (with numpy and without numpy) here.

Happy processing!

Starting with QField 2.2, users can fully rely on animation capabilities that have made their way into QGIS during its last development cycle. This can be a powerful mean to highlight key elements on a map that require special user attention.

The example below demonstrates a scenario where animated raster markers are used to highlight active fires within the visible map extent. Notice how the subtle fire animation helps draw viewers’ eyes to those important markers.

Animated raster markers is a new symbol layer type in QGIS 3.26 that was developed by Nyall Dawson. Supported image formats include GIF, WEBP, and APNG.

The second example below showcases more advanced animated symbology which relies on expressions to animate several symbol properties such as marker size, border width, and color opacity. While more complex than simply adding a GIF marker, the results achieved with data-defined properties animation can be very appealing and integrate perfectly with any type of project.

You’ll quickly notice how smooth the animation runs. That is thanks to OPENGIS.ch’s own ninjas having spent time improving the map canvas element’s handling of layers constantly refreshing. This includes automatic skipping of frames on older devices so the app remains responsive.

Oh, we couldn’t help ourselves but take the opportunity to demonstrate how nice the QField feature form layout is these days in the video above To know more about other new features in QField 2.2, go and read the release page.

Happy field mapping to all!

The lovely animal markers used in the zoo example above were made by Serbian artist Arsenije Vujovic.

Together with our partners at Lutra Consulting and Hobu, we have once again have collated your feature requests for even MORE point cloud and QGIS elevation improvements and are ready to start working on them!

Here’s a taster of what to expect if our latest crowd funding campaign is successful:

• Point cloud data management, transformation and analysis via incorporation of PDAL within the QGIS “Processing” toolbox
• Optimised handling of large datasets to speed up your 3D maps, including dynamic data loading for complex 3D scenes
• Combine multiple LAS/LAZ point cloud files into a single “virtual” point cloud for easy data management and display
• Elevation profile tool embedded into print layouts, totally customizable, more efficient and with the ability to export profile data to CSV and DXF
• General improvements to QGIS 3D map views, including an enhanced 3D measuring tool, additional camera controls and improving the configuration dialog options and functionality.

To find out more, check out the ‘Detailed proposal and deliverables’ section on the main crowd funding page. You may want to get some popcorn ready, it really is exciting!

We need your help to make this work possible. If you’d like to see these enhancements, pledge to the crowdfunding campaign before October 24, 2022.

Today’s post presents an experiment in modelling a common scenario in many IOT setups: time series of measurements at stationary sensors. The key idea I want to explore is to use MobilityDB’s temporal data types, in particular the tfloat_inst and tfloat_seq for instances and sequences of temporal float values, respectively.

For info on how to set up MobilityDB, please check my previous post.

Setting up our DB tables

As a toy example, let’s create two IOT devices (in table iot_devices) with three measurements each (in table iot_measurements) and join them to create the tfloat_seq (in table iot_joined):

CREATE TABLE iot_devices (
    id integer,
    geom geometry(Point, 4326)
);

INSERT INTO iot_devices (id, geom) VALUES
(1, ST_SetSRID(ST_MakePoint(1,1), 4326)),
(2, ST_SetSRID(ST_MakePoint(2,3), 4326));

CREATE TABLE iot_measurements (
    device_id integer,
    t timestamp,
    measurement float
);

INSERT INTO iot_measurements (device_id, t, measurement) VALUES
(1, '2022-10-01 12:00:00', 5.0),
(1, '2022-10-01 12:01:00', 6.0),
(1, '2022-10-01 12:02:00', 10.0),
(2, '2022-10-01 12:00:00', 9.0),
(2, '2022-10-01 12:01:00', 6.0),
(2, '2022-10-01 12:02:00', 1.5);

CREATE TABLE iot_joined AS
SELECT 
    dev.id, 
    dev.geom, 
    tfloat_seq(array_agg(
        tfloat_inst(m.measurement, m.t) ORDER BY t
    )) measurements
FROM iot_devices dev 
JOIN iot_measurements m
  ON dev.id = m.device_id
GROUP BY dev.id, dev.geom;

We can load the resulting layer in QGIS but QGIS won’t be happy about the measurements column because it does not recognize its data type:

Query layer with valueAtTimestamp

Instead, what we can do is create a query layer that fetches the measurement value at a specific timestamp:

SELECT id, geom, 
    valueAtTimestamp(measurements, '2022-10-01 12:02:00') 
FROM iot_joined

Which gives us a layer that QGIS is happy with:

Time for TemporalController

Now the tricky question is: how can we wire our query layer to the Temporal Controller so that we can control the timestamp and animate the layer?

I don’t have a GUI solution yet but here’s a way to do it with PyQGIS: whenever the Temporal Controller signal updateTemporalRange is emitted, our update_query_layer function gets the current time frame start time and replaces the datetime in the query layer’s data source with the current time:

l = iface.activeLayer()
tc = iface.mapCanvas().temporalController()

def update_query_layer():
    tct = tc.dateTimeRangeForFrameNumber(tc.currentFrameNumber()).begin().toPyDateTime()
    s = l.source()
    new = re.sub(r"(\d{4})-(\d{2})-(\d{2}) (\d{2}):(\d{2}):(\d{2})", str(tct), s)
    l.setDataSource(new, l.sourceName(), l.dataProvider().name())

tc.updateTemporalRange.connect(update_query_layer)

Future experiments will have to show how this approach performs on lager datasets but it’s exciting to see how MobilityDB’s temporal types may be visualized in QGIS without having to create tables/views that join a geometry to each and every individual measurement.

It’s been a while since we last talked about MobilityDB in 2019 and 2020. Since then, the project has come a long way. It joined OSGeo as a community project and formed a first PSC, including the project founders Mahmoud Sakr and Esteban Zimányi as well as Vicky Vergara (of pgRouting fame) and yours truly.

This post is a quick teaser tutorial from zero to computing closest points of approach (CPAs) between trajectories using MobilityDB.

Setting up MobilityDB with Docker

The easiest way to get started with MobilityDB is to use the ready-made Docker container provided by the project. I’m using Docker and WSL (Windows Subsystem Linux on Windows 10) here. Installing WLS/Docker is out of scope of this post. Please refer to the official documentation for your operating system.

Once Docker is ready, we can pull the official container and fire it up:

docker pull mobilitydb/mobilitydb
docker volume create mobilitydb_data
docker run --name "mobilitydb" -d -p 25432:5432 -v mobilitydb_data:/var/lib/postgresql mobilitydb/mobilitydb
psql -h localhost -p 25432 -d mobilitydb -U docker

Currently, the container provides PostGIS 3.2 and MobilityDB 1.0:

Loading movement data into MobilityDB

Once the container is running, we can already connect to it from QGIS. This is my preferred way to load data into MobilityDB because we can simply drag-and-drop any timestamped point layer into the database:

For this post, I’m using an AIS data sample in the region of Gothenburg, Sweden.

After loading this data into a new table called ais, it is necessary to remove duplicate and convert timestamps:

CREATE TABLE AISInputFiltered AS
SELECT DISTINCT ON("MMSI","Timestamp") *
FROM ais;

ALTER TABLE AISInputFiltered ADD COLUMN t timestamp;
UPDATE AISInputFiltered SET t = "Timestamp"::timestamp;

Afterwards, we can create the MobilityDB trajectories:

CREATE TABLE Ships AS
SELECT "MMSI" mmsi,
tgeompoint_seq(array_agg(tgeompoint_inst(Geom, t) ORDER BY t)) AS Trip,
tfloat_seq(array_agg(tfloat_inst("SOG", t) ORDER BY t) FILTER (WHERE "SOG" IS NOT NULL) ) AS SOG,
tfloat_seq(array_agg(tfloat_inst("COG", t) ORDER BY t) FILTER (WHERE "COG" IS NOT NULL) ) AS COG
FROM AISInputFiltered
GROUP BY "MMSI";

ALTER TABLE Ships ADD COLUMN Traj geometry;
UPDATE Ships SET Traj = trajectory(Trip);

Once this is done, we can load the resulting Ships layer and the trajectories will be loaded as lines:

Computing closest points of approach

To compute the closest point of approach between two moving objects, MobilityDB provides a shortestLine function. To be correct, this function computes the line connecting the nearest approach point between the two tgeompoint_seq. In addition, we can use the time-weighted average function twavg to compute representative average movement speeds and eliminate stationary or very slowly moving objects:

SELECT S1.MMSI mmsi1, S2.MMSI mmsi2, 
       shortestLine(S1.trip, S2.trip) Approach,
       ST_Length(shortestLine(S1.trip, S2.trip)) distance
FROM Ships S1, Ships S2
WHERE S1.MMSI > S2.MMSI AND
twavg(S1.SOG) > 1 AND twavg(S2.SOG) > 1 AND
dwithin(S1.trip, S2.trip, 0.003)

In the QGIS Browser panel, we can right-click the MobilityDB connection to bring up an SQL input using Execute SQL:

The resulting query layer shows where moving objects get close to each other:

To better see what’s going on, we’ll look at individual CPAs:

Having a closer look with the Temporal Controller

Since our filtered AIS layer has proper timestamps, we can animate it using the Temporal Controller. This enables us to replay the movement and see what was going on in a certain time frame.

I let the animation run and stopped it once I spotted a close encounter. Looking at the AIS points and the shortest line, we can see that MobilityDB computed the CPAs along the trajectories:

A more targeted way to investigate a specific CPA is to use the Temporal Controllers’ fixed temporal range mode to jump to a specific time frame. This is helpful if we already know the time frame we are interested in. For the CPA use case, this means that we can look up the timestamp of a nearby AIS position and set up the Temporal Controller accordingly:

More

I hope you enjoyed this quick dive into MobilityDB. For more details, including talks by the project founders, check out the project website.

This post is part of a series. Read more about movement data in GIS.