Monday 23 October 2017

pyQGIS Merge raster layers


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
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, '')

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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...