[GIS] pyQGIS Merge raster layers

mergepyqgisraster

I want to merge raster layers I have created in different software (SNAP). I know about merge tools in QGIS (qgis merge and gdal merge). I have problem with null values in the rasters which were created.

When I use merge tool which is in the QGIS menu I can set to ignore null values and result is correct.

enter image description here

When I use gdal merge tool (processing tool), null values are not ignored and result is not correct because there are stripes with null values in the result raster. There is no parameter to ignore null values (null data parameter does not ignore null values it just set selected value as null value)

enter image description here

These are my parameters in tools.

enter image description here

And my problem is there is no qgis merge tool in the pyqgis. Only GDAL merge tool. I checked by processing.alglist() in QGIS. I think I can use vector merge which is implemented (qgis:mergevectorlayers -> polygonize every raster, select and save only features with value 1 and merge) but I would like to "get rid of" lot of rasters (vectors) sooner and use raster merge and not to poligonize lots of rasters. Yes, the final step is polygonization but I want to polygonized merged raster.

These are my original rasters, there are terrain corrected so they were rotated and null values are around.

enter image description here

EDIT @SIGI script I have used

#!/usr/bin/env python
# -*- coding: utf-8 -*-

#import libraries
import sys, os
import qgis
from qgis.core import *
from PyQt4.QtCore import QFileInfo
from subprocess import Popen, PIPE, STDOUT

# supply path to qgis install location
QgsApplication.setPrefixPath("/usr", True)

# create a reference to the QgsApplication
qgs = QgsApplication([],True, None) 
#app.setPrefixPath("/usr", True)

# load providers
qgs.initQgis()

# import processing and set system path to processing
sys.path.append('/usr/share/qgis/python/plugins')   
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

class ProcessException(Exception):
    pass
# function to run the command
def execute(command):

    process = Popen(command.split(" "), shell=True, stdout=PIPE, stderr=STDOUT, stdin=PIPE)
    # Poll process for new output until finished
    while True:
        nextline = process.stdout.readline()
        if nextline == '' and process.poll() is not None:
            break
        # You can add here some progress infos

    output = process.communicate()[0]
    exitCode = process.returncode

    if (exitCode == 0):
        return output
    else:
        raise ProcessException(command, exitCode, output)



execute(u"gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /media/sf_dpz/dpz_data/outputs/qgis/mosaica100.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra‌​in_corrected05.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra‌​in_corrected06.tif")

Error in bash shell:

raceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 196, in qgis_excepthook
    showException(type, value, tb, None, messagebar=True)
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 107, in showException
    open_stack_dialog(type, value, tb, msg)
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 142, in open_stack_dialog
    iface.messageBar().popWidget()
AttributeError: 'NoneType' object has no attribute 'messageBar'

Original exception was:
Traceback (most recent call last):
  File "/home/lukas/mosaic_gdal.py", line 49, in <module>
    execute(u"gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /media/sf_dpz/dpz_data/outputs/qgis/mosaica100.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra‌​in_corrected05.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra‌​in_corrected06.tif")
  File "/home/lukas/mosaic_gdal.py", line 45, in execute
    raise ProcessException(command, exitCode, output)
__main__.ProcessException: (u'gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /media/sf_dpz/dpz_data/outputs/qgis/mosaica100.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra\u200c\u200bin_corrected05.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terra\u200c\u200bin_corrected06.tif', 1, '')

Best Answer

You don't need processing to execute gdal_tools you can use subprocess module to execute your command. Even Processing use the subprocess so why don't you ? here a function to use it :

import sys
from subprocess import Popen, PIPE, STDOUT
# Create a class herit from Exception class
class ProcessException(Exception):
    pass
# function to run the command
def execute(command):

    process = Popen(command.split(" "), shell=True, stdout=PIPE, stderr=STDOUT, stdin=PIPE)
    # Poll process for new output until finished
    while True:
        nextline = process.stdout.readline()
        if nextline == '' and process.poll() is not None:
            break
        # You can add here some progress infos

    output = process.communicate()[0]
    exitCode = process.returncode

    if (exitCode == 0):
        return output
    else:
        raise ProcessException(command, exitCode, output)

To use it :

try:
    execute(u"gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /path/of/your/output/mozaik.tif /path/of/your/input/terrain_corrected01.tif /path/of/your/input/terrain_corrected02.tif")
except ProcessException as p:
    print(u"something goes wrong : {}".format(p))

--EDIT-- Where is those command

Gdal comes with QGIS as you install QGIS. In Windows they are in the :

C:\path\of\OSGEO4W\bin\

In linux it's in the :

/usr/bin/
# for the manual
/usr/share/man/manl/gdal_merge.l.gz

In my opinion if you don't need the QGIS API it will be faster and easy to use it in a bash script. But it depends on your need.

Related Question