3

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, '')
zubro
  • 381
  • 2
  • 13
  • Try to build vrt with gdal in the processing tools. It has some parameters that's for no data. – SIGIS Nov 07 '17 at 06:42
  • Thanks. I can than work with virtual layer in pyqgis and save it ? – zubro Nov 07 '17 at 07:17
  • Yep definitely! – SIGIS Nov 07 '17 at 07: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 07 '17 at 08:20
  • Ok can you post the command that the menu tool (raster/merge) prompt you ? – SIGIS Nov 07 '17 at 09:17
  • QGIS/Raster/Merge: gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o ./mozaik.tif ./terrain_corrected01.tif ./terrain_corrected02.tif – zubro Nov 07 '17 at 09:48

1 Answers1

4

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.

SIGIS
  • 912
  • 7
  • 19
  • 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 07 '17 at 11:41
  • My bad I didn't declare the class error ProcessException. I edit the code. – SIGIS Nov 07 '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 07 '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 07 '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 07 '17 at 15:07
  • I have tried now gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /home/lukas/zdffgl /home/lukas/terrain_corrected05.tif /home/lukas/terrain_corrected06.tif and it was done correctly ! But I copied command to the script and there is error AttributeError: 'NoneType' object has no attribute 'messageBar'. But gdal alone works. – zubro Nov 07 '17 at 16:15
  • This error is normal because you don't have the iface instanciated as you run a standalone script. But it did not tell me why the exitCode is not 0 (everything is OK) but in shell it's OK. – SIGIS Nov 07 '17 at 20:54
  • Thank you for help and time. I have used gdal_merge without pyqgis after all, I did not know I it possible to use it without pyqgis. Could you just tell me how it works ? How it recognized gdal_merge.py in bash script just by name if the path to the gdal script is not declared anywhere and it is part of QGIS ? Also is it better to run gdal_merge in python scripts or it is ok to run it just as simple command in bash script. Thanks for all help it was really helpful. I have just thought I need to install GDAL standalone (not as part of QGIS) to use it without qgis as library in shell. – zubro Nov 08 '17 at 12:52
  • You are welcome, always a pleasure to help. I wanted to tell you to use gdal_merge without pyqgis, but it is not the question. I edit my answer about your questions. Which PPA have used ? (ubuntugis or qgis.org) – SIGIS Nov 09 '17 at 09:00
  • I have qgis.org. I will use it without QGIS. I use QGIS for another steps but I run whole process in docker and I have for every step in processing own container so I can choose for merging gdal without python. Many thanks. – zubro Nov 09 '17 at 09:23