2

I have a problem with the grass function r.mfilter outputting an empty file. I'm running the tool from qgis processing standalone.

My system:

OS: Ubuntu 16.04.2 LTS
Saga: 2.2.3
Grass: 7.0.3
QGIS: 2.18.4

Saga installed as per instructions in this question

The script:

from qgis.core import *
from qgis.gui import *
import os
import sys
from PyQt4.QtGui import *
from PyQt4.QtCore import Qt

app = QgsApplication(sys.argv, True)
QgsApplication.initQgis();
QgsApplication.setPrefixPath('/usr', True)

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
import processing.tools.general as processing

Processing.initialize()

# Processing directory
cwd = os.getcwd()
procdir = cwd + '/LST-trad/hojd/'
tmpdir = cwd + '/LST-trad/Test2/'
outdir = cwd + '/LST-trad/treecover2/'
filter = cwd + "/LST-trad/qgis5x5filter_1or0_div1.txt"
files = os.listdir(procdir)
# Filter files
filefilt =[]
for f in files:
    if f.endswith(".tif"):
        filefilt.append(f)

print "help:"
#processing.alghelp("saga:rastercalculator")
processing.alghelp("grass:r.mfilter")

# Coordinates
arglist = []
for file in filefilt:
    x1 = int(file[7] + file[11:13]+"000")
    x2 = x1 + 25000
    y1 = int(file[4:6] + file[9:11] +"000")
    y2 = y1 + 25000
    coords =  str(x1) + "," + str(x2) + "," + str( y1) + "," + str(y2)
    outfile = "COV_" + file[4:]
    arglist.append([procdir+file, coords, outdir+outfile, outfile])

#print arglist

formula1 =('ifelse(a > 50, 1, 0)')
formula2 = ('a * 4')

for arg in arglist:
    tmp1 = tmpdir + "t1_" + arg[3]
    tmp2 = tmpdir + "t2_" + arg[3]
    print "saga"
    processing.runalg("saga:rastercalculator", arg[0], None, formula1, True, 1, tmp1)
    print "grass"
    processing.runalg("grass:r.mfilter", tmp1, filter,1,False,arg[1], 0, tmp2)
    print "rascal"
    processing.runalg("saga:rastercalculator", tmp2, None, formula2, True, 1, arg[2])

QgsApplication.exitQgis();

The filter file:

MATRIX    5
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
DIVISOR   1
TYPE      P 

output form alghelp(r.mfilter):

ALGORITHM: r.mfilter - Performs raster map matrix filter.
    input <ParameterRaster>
    filter <ParameterFile>
    repeat <ParameterNumber>
    -z <ParameterBoolean>
    GRASS_REGION_PARAMETER <ParameterExtent>
    GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber>
    output <OutputRaster>

One thing that confuses me is that the parameters listed here does not match the parameters listed in the grass documentation

2 Answers 2

1

I can't test your code, but you need to follow the syntax returned by processing.alghelp().

However, it's unclear which doesn't work with your code: maybe you need to add the outputs from the processing operations as objects (otherwise QGIS will interpret them as strings).

Try replacing the last for-loop of your code with something like this:

for arg in arglist:
    tmp1 = tmpdir + "t1_" + arg[3]
    tmp2 = tmpdir + "t2_" + arg[3]
    print "saga"
    first_op = processing.getObject(arg[0])
    processing.runalg("saga:rastercalculator", first_op, None, formula1, True, 1, tmp1)
    print "grass"
    second_op = processing.getObject(tmp1)
    processing.runalg("grass:r.mfilter", second_op, filter,1,False,arg[1], 0, tmp2)
    print "rascal"
    third_op = processing.getObject(tmp2)
    processing.runalg("saga:rastercalculator", third_op, None, formula2, True, 1, arg[2])

Note: I don't understand why you are running a loop without looping inside arglist: in fact, you directly call arg[0], arg[1] and arg[2] (I suppose you know what are you doing).

8
  • I've tried that now, the first operation runs fine, but r.mfilter still outputs an empty file. I don't understand the GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber> parameter. And since it's not in the grass documentation for the tool I don't know where too find an explanation. Re: Note: arglist should really be called inputfiles or something. Each element is a list with filename, extent etc.
    – Iamlukesky
    Commented Mar 27, 2017 at 12:04
  • The GRASS_REGION_PARAMETER should define the size of the cell for the output raster: if you leave 0 as value, it will use the size of the input raster, otherwise it will assign a custom value specified by the user. Is there a creation of an empty file or do you get an error without writing anything on the disl? Furthermore, did you try using the algorithm from GRASS7 instead? I mean, the "grass7:r.mfilter" algorithm.
    – mgri
    Commented Mar 27, 2017 at 12:51
  • I just tried running the tool within QGIS, works fine, output is just what want. I took the exact snippet used for running the tool from the Processing History: processing.runalg("grass7:r.mfilter","/home/johnnie/GiB/Projekt/qgisHeadless/LST-trad/Test2/t1_COV_66_6_0025.tif","/home/johnnie/GiB/Projekt/qgisHeadless/LST-trad/qgis5x5filter_1or0_div1.txt",1,False,"625000.0,650000.0,6600000.0,6625000.0",0,None). So specifying the input as strings seems to be just fine. But when I substitue the line in the script with this line I still get an output of an empty file, no error message.
    – Iamlukesky
    Commented Mar 27, 2017 at 13:31
  • 1
    The GRASS_REGION_PARAMETER is the extent, not the size. See for an extent() example gis.stackexchange.com/questions/129639/…
    – markusN
    Commented Mar 29, 2017 at 18:06
  • 1
    @markusN you are obviously right, I made some confusion! Thanks for the clarification! =)
    – mgri
    Commented Mar 29, 2017 at 19:25
1

This question turned out to have the answer. It seems like the Grass tools are really picky about how you set the extent. It's somewhat annoying to not be able to set it to None, but at least this works.

I've also set the output type to None for the intermediate steps to avoid creating temporary middle-files as recommended in the comments.

from qgis.core import *
from qgis.gui import *
import os
import sys
from PyQt4.QtGui import *
from PyQt4.QtCore import Qt, QFileInfo

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
#import processing.tools.general as processing
from processing.tools.dataobjects import *
from processing.tools.general import *
from processing.tools.vector import *
from processing.tools.raster import *
from processing.tools.system import *

app = QgsApplication(sys.argv, True)
QgsApplication.initQgis();
QgsApplication.setPrefixPath('/usr', True)
Processing.initialize()

# Processing directory
cwd = os.getcwd()
procdir = cwd + '/LST-trad/hojd/'
outdir = cwd + '/LST-trad/treecover2/'
filterfile = cwd + "/LST-trad/qgis5x5filter_1or0_div1.txt"
files = os.listdir(procdir)
#files = [f for f in os.listdir(procdir) if os.path.isfile(f)]
# Filter files
#filefilt =[]
#for f in files:
    #if f.endswith(".tif"):
        #filefilt.append(f)

print "startar"

#print "help:"
#alghelp("saga:rastercalculator")
#alghelp("grass:r.mfilter")

inputraster = procdir + "/THL_66_6_0025.tif"
outputraster = outdir + "/uttest.tif"

fileInfo = QFileInfo(inputraster)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(inputraster, baseName)

extent = rlayer.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

print "rascal 1"
outputs_SAGARASTERCALCULATOR_1=runalg('saga:rastercalculator', inputraster,[],'ifelse(a > 50, 1, 0)',True,1,None)

print "mfilter"
outputs_GRASS7R_MFILTER_1=runalg('grass7:r.mfilter', outputs_SAGARASTERCALCULATOR_1['RESULT'],filterfile,1.0,False,"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),0.0,None)

print "rascal 2"
outputs_SAGARASTERCALCULATOR_2=runalg('saga:rastercalculator', outputs_GRASS7R_MFILTER_1['output'],[],'a * 4',True,1,outputraster)

print "done"

QgsApplication.exitQgis();

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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