[GIS] PyQGIS save raster as rendered image, then use GDAL tools on it

gdalgdal-translategdalwarppyqgisqgis

I have tried to use the code from QGis Save Raster as Rendered Image. i is a raster layer:

pipelayer = i
pipeextent = pipelayer.extent()
pipewidth, pipeheight = (pipelayer.width(),
                         pipelayer.height())
piperenderer = pipelayer.renderer()
pipeprovider = pipelayer.dataProvider()
crs = pipelayer.crs().toWkt()
pipe = QgsRasterPipe()
pipe.set(pipeprovider.clone())
pipe.set(piperenderer.clone())
pipedFile = os.path.join(tempfile.gettempdir(),
                         safeLayerName + '_pipe.tif')
print pipedFile
file_writer = QgsRasterFileWriter(pipedFile)
file_writer.writeRaster(pipe,
                        pipewidth,
                        pipeheight,
                        pipeextent,
                        pipelayer.crs())

in_raster = pipedFile
prov_raster = os.path.join(tempfile.gettempdir(),
                           'json_' + safeLayerName +
                           '_prov.tif')
out_raster = dataPath + '.png'
crsSrc = i.crs()
crsDest = QgsCoordinateReferenceSystem(4326)
xform = QgsCoordinateTransform(crsSrc, crsDest)
extentRep = xform.transform(i.extent())
extentRepNew = ','.join([unicode(extentRep.xMinimum()),
                         unicode(extentRep.xMaximum()),
                         unicode(extentRep.yMinimum()),
                         unicode(extentRep.yMaximum())])
processing.runalg("gdalogr:warpreproject", in_raster,
                  i.crs().authid(), "EPSG:4326", "", 0, 1,
                  5, 2, 75, 6, 1, False, 0, False, "",
                  prov_raster)
processing.runalg("gdalogr:translate", prov_raster, 100,
                  True, "", 0, "", extentRepNew, False, 0,
                  0, 75, 6, 1, False, 0, False, "",
                  out_raster)

It doesn't work, instead giving me the unstyled raster as final output. I really don't understand the problem, because the output of writeRaster (the file pipedFile) is styled, so the process from the answer linked to above is working. It's just that when I try to use the output image to run through GDAL's warpreproject and translate, it somehow reverts to the unstyled raster.

What have I done wrong?

UPDATE: The output of QgsRasterFileWriter is styled. The output of warpreproject is not styled, but the .tif is accompanied by a .aux.xml file which appears to have colour information in it:

<PAMDataset>
  <PAMRasterBand band="1">
    <Histograms>
      <HistItem>
        <HistMin>-0.498046875</HistMin>
        <HistMax>255.498046875</HistMax>
        <BucketCount>256</BucketCount>
        <IncludeOutOfRange>0</IncludeOutOfRange>
        <Approximate>1</Approximate>
        <HistCounts>314|4|3|2|3|5|4|4|5|8|10|9|16|16|23|26|37|46|58|62|69|77|77|94|127|94|131|136|133|162|169|172|171|184|200|186|207|196|186|177|179|194|182|182|178|177|181|195|196|196|195|195|179|226|206|215|193|197|202|212|207|206|232|204|229|253|233|240|248|234|234|268|240|238|259|286|268|287|259|264|255|267|268|271|256|321|277|284|317|286|290|277|312|334|319|325|323|310|349|357|347|326|337|338|336|383|374|349|411|382|382|417|406|405|414|427|404|434|447|430|468|424|446|442|450|461|448|458|457|459|412|414|469|466|443|475|464|472|481|519|504|458|473|481|514|523|522|494|542|580|604|623|686|648|707|763|808|799|819|853|932|1046|1062|1118|1145|1215|1293|1393|1335|1392|1390|1317|1232|1247|1215|1088|1032|896|847|731|708|666|608|634|556|500|484|456|433|408|406|401|386|388|414|414|387|382|340|346|382|385|368|348|334|376|359|318|339|326|357|315|319|357|333|328|329|310|344|295|300|351|347|324|358|349|363|361|347|396|386|433|395|434|479|473|523|550|644|655|672|778|882|1014|1183|1400|1681|2075|2531|3529|4851|7235|10793|60773|7297|1999</HistCounts>
      </HistItem>
    </Histograms>
    <Metadata>
      <MDI key="STATISTICS_MAXIMUM">255</MDI>
      <MDI key="STATISTICS_MEAN">203.31482680479</MDI>
      <MDI key="STATISTICS_MINIMUM">0</MDI>
      <MDI key="STATISTICS_STDDEV">63.158702013579</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>

The output of translate is not styled, and has no additional file. However, since it's a PNG not a GeoTIFF, that's perhaps the explanation.

I've tried looking at the expand option of translate, but that gives me image files which do not display in a web browser – I've not investigated why.

UPDATE 2: No, using expand results in no output file at all. Wondering about trying pct2rgb.

UPDATE 3: I'm wondering if file locking is causing the problems. Could it be?

Best Answer

You should understand that when you save the file from the QgsRenderer as in your code:

piperenderer = pipelayer.renderer()
pipe.set(piperenderer.clone())
file_writer.writeRaster(pipe...

You are not saving the raster styled, instead you are saving the actual RGB file that represents that visualization, so it is not classified or styled, but it is like you've taken a print screen of the rendered image.

So gdalwarp and translate can't change this, as it is hardcoded in the image's bands, what must be happening is that they are clipping by a single band or so, I have to check it thorughly.

EDIT

I can see in the xml file you provided that the raster has now only one band

<PAMRasterBand band="1">

There should be the tags:

<PAMRasterBand band="2">
<PAMRasterBand band="3">

So the problem is in gdalogr:warpreproject. I used the same code you provided in a raster I have but the output is RGB, I don't know what's going on, maybe you should give a sample of your raster file.

NOTE

The xml has no style data, as I said the raster you firstly saved is not a styled layer, but an actual RGB image representing the styled raster.

Related Question