2

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
    layer_entry.append(raster_entry)
    num+=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)
        calc.processCalculation()
        # Reset Values
        num = 1
        layer_entry = []
    start_date+=delta

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?

closed as too broad by PolyGeo Jun 7 '18 at 6:51

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

1

In SQL only, you could try to compute the average raster with ST_Union():

SELECT date_part('year', fdate) year, 
       date_part('month', fdate) month, 
       ST_Union(rast, 'mean') rast
FROM model
WHERE layer = 1 
GROUP BY year, month
ORDER BY year, month;

Not the answer you're looking for? Browse other questions tagged or ask your own question.