PyQGIS: selecting WMS timestamp in a standalone script

Question

I'm trying to create a standalone python script that uses PyQGIS to create a map with a raster layer read from WMS.

My idea was to create a QgsRasterLayer with provider set to wms with the time argument encoded into the URL.

To test, I've just tried adding them in the python console within QGIS:

# default time
url1 = "layers=MCD43C3_M_BSA&styles=rgb&bbox=-180,-90,180,90&height=512&width=512&crs=CRS:84&format=image/png&url=https://neo.gsfc.nasa.gov/wms/wms?version%3D1.3.0"
selected time 2017-01-01
url2 = "layers=MCD43C3_M_BSA&styles=rgb&bbox=-180,-90,180,90&height=512&width=512&crs=CRS:84&format=image/png&time=2017-01-01&url=https://neo.gsfc.nasa.gov/wms/wms?version%3D1.3.0"
iface.addRasterLayer(url1, "default", "wms")
iface.addRasterLayer(url2, "selection", "wms")

However, I get the same raster both times, and the developer panel indicates the time parameter is not included in the query:

Is there a simple way I'm missing where I can select a date out of the available time steps from a WMS layer?

Edit:

til_b's excellent answer below has shown that one can use owslib to interact with the WMS server and with that easily select dates from the layer's time axis.

The issue is the intended use:

I want to use this WMS source to create a map using PyQGIS. I could indeed try to obtain the image some other way and then plug in into PyQGIS for mapping, but that seems overly complicated and has some gotchas:

1. I need to define the image shape beforehand (I suppose that could be calculated from extent, resolution, dpi etc.)
2. Some servers limit the max image shape to tiles (e.g. 512) so I'd need to get multiple tiles and stitch them together myself
3. I need to encode the geoinfo myself if I want to map it together with other geospatial data

Doing all of this myself feels a lot like re-inventing the wheel - particularly it's already being handled by PyQGIS with the notable issue that I can't select a date!

I really think there must be a way to do this with PyQGIS since you can clearly do it with the GUI. Either it's a bit more elaborate (maybe using the Temporal Controller perhaps), or I'm just missing the obvious! (or neither and it's a bug ... )

Answer 1

Use owslib per their instructions at https://geopython.github.io/OWSLib/usage.html#wms .

OWSlib passes any additional HTTP parameters directly to the server (checked with a simple netstat -lp 9999 listening socket).

Order of commands to test the behaviour:

In a linux shell: netstat -lp 9999. Leave running.

In python:

from owslib.wms import WebMapService
wms = WebMapService('http://localhost:9999/?bar=baz', version='1.3.0')

The request that arrives on the netstat socket is

GET /?bar=baz&service=WMS&request=GetCapabilities&version=1.3.0 HTTP/1.1
Host: localhost:9999
User-Agent: python-requests/2.22.0
Accept-Encoding: gzip, deflate
Accept: */*
Connection: keep-alive

It contains the bar=baz parameters.

You should then be able to produce an output image using "plain" OWSlib with the extra parameter arrriving at the server.

It the server supports image/geotiff as output format you can plug that directly into QGIS, if not you'll have to write a world file for the jpeg or png image with the following code (ripped from a real project, redacted, so may need some tweaking to run, but you can see the ideas):

def wf_params(ul,lr,size):
    """ul: upper left corner.
lr: lower right corner.
Only for straight-edge rectangles!"""
    ### World file writing from 
    ### https://gis.stackexchange.com/questions/269581/creating-tfw-world-files-from-tiff-coordinates-in-excel
    width, height = size
    xul,yul = ul  # -> coordinates of the upper left corner
    xlr,ylr = lr  # -> coordinates of the lower rigth corner
#print("%.3f %.3f" % (width,height))
#print("%.3f %.3f" % (xul,yul))
#print("%.3f %.3f" % (xlr,ylr))

A = (xlr-xul)/width
D = 0
B = 0
E = (ylr-yul)/height 
#center of the upper left pixel 
C = xul + (A * .5)
F = yul + (E * .5)
return (A,D,B,E,C,F)


def writeWF(params,fullpath):
    """Write a world file from the parameters produced by wf_param"""
    print("Writing worldfile to %s" % fullpath)
    f = open(fullpath,'w')
    for p in params:
        f.write("%.3f\n" % p);
    f.close()
...
bb = <the bounding box used for the owslib WMS request, given as [lowerleftX, lowerleftY, upperrightX, upperrightY>
size of the requested WMS image
outsize = 512
upperleft  = (bb[0],bb[3])
lowerright = (bb[2],bb[1])
parms = wf_params(upperleft,lowerright,outsize)
writeWF(parms ,full_outname_world)

Answer 2

Looking at the source code, it turns out to select a specific date you need

1. "time-aware" layer
2. select the time using the start/end format

Regarding 1), I've linked in the comments to a thread in the dev mailinglist which explains the required components.

Essentially it boils down to adding the following query arguments to your url:

• allowTemporalUpdates=true
• temporalSource=provider
• type=wmst
• timeDimensionExtent

with the timeDimensionExtent to be defined by the capabilities from the server.

Regarding 2), if all is needed is one date, we can just use the same date as both start and end, which results in a query time=2017-01-01/2017-01-01 (for this example).

Using the shortcut of copying the time aware URL from the GUI (pointed out by @til_b), everything seems to work:

# original non time-aware URL
url1 = "layers=MCD43C3_M_BSA&styles=rgb&bbox=-180,-90,180,90&height=512&width=512&crs=CRS:84&format=image/png&url=https://neo.gsfc.nasa.gov/wms/wms?version%3D1.3.0"
new time aware URL with time query 2017-01-01
url2 = "allowTemporalUpdates=true&contextualWMSLegend=0&crs=CRS:84&dpiMode=7&featureCount=10&format=image/png&layers=MCD43C3_M_BSA&styles&temporalSource=provider&timeDimensionExtent=2000-02-01/2016-04-01/P1M,2016-06-01/2017-02-01/P1M&type=wmst&time=2017-01-01/2017-01-01&url=https://neo.gsfc.nasa.gov/wms/wms?version%3D1.3.0"
iface.addRasterLayer(url1, "default", "wms")
iface.addRasterLayer(url2, "selection", "wms")

Not only are the two layers different now, the developer panel is also showing the TIME query as expected:

Answer 3

It's probably as simple as sticking the parameter urlencoded into the QGIS layer URI string.

url2 = "layers=MCD43C3_M_BSA&styles=rgb&bbox=-180,-90,180,90&height=512&width=512&crs=CRS:84&format=image/png&url=https://neo.gsfc.nasa.gov/wms/wms?version%3D1.3.0%3Dtime=2017-01-01"

Your original code had the parameter time as a part of the URL for QGIS addRasterLayer. You want it as part of the WMS URL inside the URL that gets to addRasterLayer.

From the "Add WMS Layer" dialog it works:

So what you can do is add the layer via the dialog and then get its URI from the python console. Select the layer, open the python console and run iface.activeLayer().dataProvider().uri(). I am not really sure how to get that URI in a format to stick back into pyqgis to create a URI programmatically. I have done that for MSSQL layers, but never for WMS layers...