QGIS Server and the WMS protocol with Time dimension

October 06, 2020

QGIS Server, WMS, Time, raster, vector, Python plugin

Paul Blottiere

Introduction

QGIS Server is a map server-based on the QGIS core library and rendering engine-which provides numerous classical services like WMS, WFS, WCS, WMTS and lastly OGC API Features (implemented by the QCooperative) since QGIS 3.10. Furthermore, the certification process has been reached for WMS 1.3.0 since QGIS 3.0 and has also been renewed for the LTR version 3.10.

However, there are several kind of extensions for the WMS protocol, especially for dimensions. This way, some specific optional parameters may be supported like ELEVATION or TIME. Today, we’re going to take a closer look to the last point: the WMS protocol with Time dimension.

QGIS Server 3.10 is shipped with a basic dimension support (implemented by 3Liz) for vector layers. Unfortunately, raster layers are left behind… but nothing is impossible with Python plugins :).

OGC Badge

In this article, we’re going to take a look on the time dimension configuration in case of vector layers. And last but not least, we’ll see how to write a basic Python plugin to support the TIME parameter for raster layers in case of simple use cases.

WMS with Time dimension and vector layers

As its name suggests, the WMS protocol is useful as soon as the underlying data has a temporal characteristic. In our case, we’re going to use population changes in Brittany thanks to a GeoJSON data.

To use the concept of dimension in QGIS Server, we need to have a dedicated field. So we have to preprocess the data above in order to create a new layer with a sepicific time field (which may be named differently). During this stage, we only keep the population for 1999, 2010, 2014 and 2017.

preprocess.py

vl = QgsVectorLayer("Polygon", "temp", "memory")
vl.startEditing()
vl.addAttribute(QgsField("code", QVariant.String))
vl.addAttribute(QgsField("time", QVariant.String))
vl.addAttribute(QgsField("pop", QVariant.Int))
vl.updateFields()

pop_layer = iface.activeLayer()
idx_level = pop_layer.fields().indexOf("level")
idx_code = pop_layer.fields().indexOf("code_geo")
idx_1999 = pop_layer.fields().indexOf("p_pop1999")
idx_2010 = pop_layer.fields().indexOf("p_pop2010")
idx_2014 = pop_layer.fields().indexOf("p_pop2014")
idx_2017 = pop_layer.fields().indexOf("p_pop2017")

for feature in pop_layer.getFeatures():
    attrs = feature.attributes()

    if attrs[idx_level] != "Commune":
        continue

    code = attrs[idx_code]
    pop_1999 = ("1999", attrs[idx_1999])
    pop_2010 = ("2010", attrs[idx_2010])
    pop_2014 = ("2014", attrs[idx_2014])
    pop_2017 = ("2017", attrs[idx_2017])

    for pop in [pop_1999, pop_2010, pop_2014, pop_2017]:
        f = QgsFeature()
        f.setGeometry(feature.geometry())
        f.setAttributes([code, pop[0], pop[1]])
        vl.addFeature(f)

crs = QgsCoordinateReferenceSystem("epsg:4326")
QgsVectorFileWriter.writeAsVectorFormat(vl, "population_time", "UTF-8", crs, "GeoJSON")

Finally we just have to adjust the symbology and configure a Time dimension pointing to the time field of the layer in the QGIS Server properties. Some more parameters are also availble for a finer configuration (like End attribute).

OGC Badge

Once QGIS Server is started, we’re able to run a GetMap request with a TIME parameter. For example, we can retrieve the population map of Brittany in 1999 thanks to the next request:

http://qgisserver? \
SERVICE=WMS \
&REQUEST=GetMap \
&WIDTH=800 \
&HEIGHT=400 \
&CRS=EPSG:4329 \
&BBOX=46.9917,-5.6523,49.1485,-0.1043 \
&VERSION=1.3.0 \
&LAYER=population_time \
&TRANSPARENT=TRUE \
&TIME=1999

And by playing with the redlining mechanism of QGIS Server (concretely HIGHLIGHT_ parameters), we can add the year on top of the map.

OGC Badge

WMS with Time dimension and raster layers

In case of raster layers, there isn’t an in-core solution but we can implement a very simple Python plugin based on the QgsServerFilter class. The idea is to catch and interpret the TIME parameter (sent by the client) with the Python plugin and update the LAYER parameter accordingly to retrieve the desired layers. This way the QGIS Server process itself isn’t even aware of the TIME parameter and only use the usual LAYER information.

Firstly, we use the sentinelsat Python API to download some S2 rasters on a specific area by using the tile number and the relative orbit number (unfortunately you need to create an account).

download_s2.py

from collections import OrderedDict
from sentinelsat import SentinelAPI

api = SentinelAPI(USER, PASSWORD)

query_kwargs = {
        'platformname': 'Sentinel-2',
        'producttype': 'S2MSI1C',
        'tileid': '30UUU',
        'relativeorbitnumber': 37,
        'date': ('NOW-25DAYS', 'NOW')}

products = OrderedDict()
kw = query_kwargs.copy()
pp = api.query(**kw)
products.update(pp)

api.download_all(products)

Once the download has finished, we can open the next true color rasters in QGIS Desktop:

T30UUU_20200901T112119_TCI.jp2
T30UUU_20200906T112121_TCI.jp2
T30UUU_20200911T112119_TCI.jp2
T30UUU_20200916T112121_TCI.jp2
T30UUU_20200921T112119_TCI.jp2

OGC Badge

The naming convention for S2 images is very simple and contains the datetime. So if we don’t want to use custom names for our layers, a very basic implementation for our plugin is to use this characteristic. Indeed, we can configure our web client to use a layer named T30UUU. Then, if the Python plugin receives TIME=2020-09-01 11:21:19, the layer name is updated accordingly: LAYER=T30UUU_20200901T112119_TCI.

The Python implementation for such a basic need is very simple. Of course the implementation may be improved to be more robust and generic, especially to support several layers.

timefilter.py

class TimeFilter(QgsServerFilter):

    def __init__(self, serverIface):
        super().__init__(serverIface)

    def requestReady(self):
        request = self.serverInterface().requestHandler()
        params = request.parameterMap()

        dt = QDateTime.fromString(params['TIME'], "yyyy-MM-dd hh:mm:ss")
        name = "{}_{}_TCI".format(params['LAYER'], dt.toString("yyyyMMddThhmmss"))

        request.setParameter('LAYER', name)

In other cases, we don’t want to have the naming constraint. Then even if the underlying idea is the same, we have to adjust the mechanism allowing to retrieve the correct layer’s name. So the first thing to do is to set up the temporal property on each layer. A script can easily do it automatically in case of a lot of images.

OGC Badge

To gain some speed, a group is created for each tile (named T30UUU in this case) and the web client is configured to use the group as the layer name. This way the Python plugin will iterate over each layer of the group instead of all layers of the project. The temporal setting is finally used to retrieve the correct layer.

timefilter_v2.py

class TimeFilter(QgsServerFilter):

    def __init__(self, serverIface):
        super().__init__(serverIface)

    def requestReady(self):
        request = self.serverInterface().requestHandler()
        params = request.parameterMap()

        dt = QDateTime.fromString(params['TIME'], "yyyy-MM-dd hh:mm:ss").date()
        project = QgsConfigCache.instance().project(params['MAP'])

        gp = project.layerTreeRoot().findGroup(params['LAYER'])
        if not gp:
            return

        for tree_layer in gp.findLayers():
            layer = tree_layer.layer()
            if layer.type() != QgsMapLayer.RasterLayer:
                continue

            prop = layer.temporalProperties()
            begin = prop.fixedTemporalRange().begin().date()

            if begin.year() != dt.year() or begin.month() != dt.month():
               continue

            request.setParameter('LAYER', layer.name())

Conclusion

As of now, QGIS Server is lacking an in-core solution for rasters and the WMS protocol Time dimension but in case of simple scenarios, a Python plugin may be implemented. But it may quickly become inefficient in case of a huge number of raster layers, especially with the second approach. Actually we can still create some groups and subgroups to optimize the research of the good layer, but the whole mechanism doesn’t scale-up. Moreover, the GetCapabilities document is not up to date with the time dimension, but it’s not necessarily an issue for simple use cases.

Some discussions are currently ongoing in the QGIS underground world, but visibility (and funds) are still missing. Let us now cross our fingers for the future :).