2

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, '')
  • Try to build vrt with gdal in the processing tools. It has some parameters that's for no data. – SIGIS Nov 7 '17 at 6:42
  • Thanks. I can than work with virtual layer in pyqgis and save it ? – zubro Nov 7 '17 at 7:17
  • Yep definitely! – SIGIS Nov 7 '17 at 7:25
  • @SIGIS I have tried to build virtual layer but there is also a problem, If I choose option stacklayer it creates mosaic correctly but it creates a band for each raster. When I polygonize the virtual raster, only first band is polygonized. If I do not choose stacklayer, it creates one band but same wrong output as GDAL. – zubro Nov 7 '17 at 8:20
  • Ok can you post the command that the menu tool (raster/merge) prompt you ? – SIGIS Nov 7 '17 at 9:17
3

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.

  • Thank you for help :) I have a error using your script: Traceback (most recent call last): ` File "<input>", line 1, in <module> File "/tmp/tmpaEv06m.py", line 20, in execute raise ProcessException(command, exitCode, output) NameError: global name 'ProcessException' is not defined`, I have copied your script to the editor in qgis and ran and then I executed the command: execute(u"gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o ./mozaik.tif ./terrain_corrected05.tif . Have I done it wrong ? – zubro Nov 7 '17 at 11:41
  • My bad I didn't declare the class error ProcessException. I edit the code. – SIGIS Nov 7 '17 at 11:47
  • This is my error know: Traceback (most recent call last): File "<input>", line 1, in <module> File "/tmp/tmp8SKgrh.py", line 24, in execute raise ProcessException(command, exitCode, output) ProcessException: ('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/terrain_corrected05.tif /media/sf_dpz/dpz_data/outputs/snap/terrain_correction/terrain_corrected06.tif', 1, ''). It is same command as used in tool window. I added prtscr of what I tried to my post. I think I am missing smtng. – zubro Nov 7 '17 at 12:03
  • What path do you mean ? Files paths or wrong system path to the gdal_merge ? Files paths are correct because I used this paths in tool. Maybe I need to set path to the gdal_merge.py in the script ? – zubro Nov 7 '17 at 14:20
  • I have used full paths as can be seen in the error in my comment above. I have tried bash shell but there is same error, I have added the error to the post. – zubro Nov 7 '17 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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