Raster Calculator for monthly input rasters from postgis

by Biplov   Last Updated July 17, 2017 10:22 AM

I have a PostGIS table (soil_moist) with columns as

id (integer), rid (integer), fdate (date), rast(raster), layer(integer)

The table is updated every day for 3 layers (layer value as 1,2 and 3). Hence, there are 3 rows added every day for the same area. The table contains the record from 2009-01-01 up to 2016-12-31. My task is to calculate the monthly average raster for layer 1 from 2009-01.

So, I created table as soil_moist_avg and then tried to accomplish this using Python and QGIS. My code snippet is as follows:

# -*- coding: utf-8 -*-

from datetime import timedelta, date
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

config = {
    "dbname"  : "model",
    "user"    : "postgres",
    "host"    : "localhost",
    "password": "******"

start_date = date(2009, 1, 1)
end_date = date(2017, 1, 1)
delta = timedelta(1)

initial_date = start_date
layer_entry = []
num = 1
while start_date < end_date:
    connString_ = "PG: dbname=%(dbname)s host=%(host)s user=%(user)s password=%(password)s port=5432 mode=2 column=rast table=soil_moist where=%(where)s" % \
                      {"dbname": config["dbname"],
                       "host": config["host"],
                       "user": config["user"],
                       "password": config["password"],
                       "where": 'fdate=\'%s\' AND layer=1' % start_date
    layer = QgsRasterLayer(connString_, "%s_%s" % (str(start_date), str(num)))
    raster_entry = QgsRasterCalculatorEntry()
    raster_entry.ref = 'layer@%s' % num
    raster_entry.raster = layer
    raster_entry.bandNumber = 1
    if ((start_date + delta).month > start_date.month) or ((start_date + delta).year > initial_date.year):
        output_path = "D:/output/%s.tif" % start_date
        expression = "("
        for idx, layer_ in enumerate(layer_entry):
            expression += "layer@%s*" % str(idx+1)
        expression += "1)"
        expression += "/%s" % str(len(layer_entry))
        calc = QgsRasterCalculator(expression, output_path, "GTiff", layer.extent(), layer.width(), layer.height(), layer_entry)
        # Reset Values
        num = 1
        layer_entry = []

Here, I am trying to dump the raster as tif for the monthly average and later thought of dumping into PostGIS database. However, the layer returned from the connection string is empty (I checked with layer.width() and layer.height()).

What am I doing wrong here?

Is it because of the where clause?

Also, I would like to know if there is a better way to accomplish this task, to directly dump the output raster into PostGIS raster data type?

Related Questions

How to sum multiple rasters in pyqgis?

Updated December 29, 2016 08:09 AM

Why QgsRasterCalculator crashes Qgis unexpectedly?

Updated December 30, 2016 08:09 AM

QGIS Script for Raster Equation

Updated January 26, 2017 14:09 PM