6

I am doing NDVI analysis on a farm image. After NDVI, I need to convert the image into singleband pseudocolor image using following Python code.

def calculate():
file_name =  CUR_DIR + "\sampleodm.tif"
out_file = CUR_DIR + "\sample_output.tif"
out_file_ndvi = CUR_DIR + "\out_ndvi.tif"
try:

    file_info = QFileInfo(file_name)
    rasterName = file_info.baseName()
    raster = QgsRasterLayer(file_name, rasterName)

    if not raster.isValid(): print 'Invalid raster'

    ir = QgsRasterCalculatorEntry()
    r = QgsRasterCalculatorEntry()

    print ir, r

    ir.raster = raster
    r.raster = raster

    ir.bandNumber = 4
    r.bandNumber = 3

    ir.ref = rasterName + "@4"
    r.ref = rasterName + "@3"

    references = (ir.ref, r.ref, ir.ref, r.ref)
    exp = "1.0 * (%s - %s) / (%s + %s)" % references
    print exp
    print 'here2'

    output = out_file

    e = raster.extent()
    w = raster.width()
    h = raster.height()
    entries = [ir, r]

    ndvi = QgsRasterCalculator(exp, output, "GTiff", e, w, h, entries)

    ndvi_res = ndvi.processCalculation()
    print ndvi_res

    file_info_out = QFileInfo(output)
    rasterName_out = file_info_out.baseName()
    lyr = QgsRasterLayer(output, rasterName_out)
    if lyr.isValid(): print 'valid raster'

    min = 0.0341278
    max = 0.56844
    range = max-min

    algorithm = QgsContrastEnhancement.StretchToMinimumMaximum
    limits = QgsRaster.ContrastEnhancementMinMax
    lyr.setContrastEnhancement(algorithm, limits)


    i = []
    i.append(QgsColorRampShader.ColorRampItem(0, QColor(0,0,0,0), 'NODATA'))
    i.append(QgsColorRampShader.ColorRampItem(0.0341278, QColor(120,69,25,255), 'Lowest Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.168, QColor(255,178,74,255), 'Lower Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.244, QColor(255,237,166,255), 'Low Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.301, QColor(173,232,94,255), 'Moderate Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.435, QColor(135,181,64,255), 'High Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.494, QColor(3,156,0,255), 'Higher Biomass'))
    i.append(QgsColorRampShader.ColorRampItem(0.56844, QColor(1,100,0,255), 'Highest Biomass'))

    c = QgsColorRampShader()
    c.setColorRampType(QgsColorRampShader.INTERPOLATED)
    c.setColorRampItemList(i)

    s = QgsRasterShader()
    s.setRasterShaderFunction(c)
    print lyr.rasterType()
    ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), lyr.rasterType(), s)
    lyr.setRenderer(ps)

    QgsMapLayerRegistry.instance().addMapLayer(lyr)

    extent = lyr.extent()
    width, height = lyr.width(), lyr.height()
    print width, height
    renderer = lyr.renderer()
    print renderer.type()
    provider = lyr.dataProvider()
    crs = lyr.crs()     

    pipe = QgsRasterPipe()
    pipe.set(provider.clone())        
    pipe.set(renderer.clone())

    file_writer = QgsRasterFileWriter(out_file_ndvi)
    print file_writer
    print 'write'
    # upto here everything is working fine and program is generating black and white image.
    error = file_writer.writeRaster(pipe, width, height, provider.extent(), crs)
    # print error
    print 'Completed NDVI analysis!'


except Exception as e:
    print 'Exception: ', e

Consider every module is imported in program.
line error = file_writer.writeRaster(pipe, width, height, provider.extent(), crs) is causing Python.exe has stopped error.

I tried this program on QGIS Desktop application's python console and it is saving the singleband pseudocolor image perfectly.

How do I rectify the error?

Matt
  • 16,843
  • 3
  • 21
  • 52

1 Answers1

1

I would recommend adding QgsApplication at the beginning of the code when using Qgis

qgs = QgsApplication([], False)
QgsApplicationsetPrefixPath.setPrefixPath("{path to qgis}/qgis", True)
qgs.initQgis()

Also, Deleting the layer after writing the program is important to improve memory usage.

del lyr

Then at the end of the program, add:

qgs.exitQgis() # Should only run once per program.

Finally, If you are using Qgis 3.8 or higher then file_writer.writeRaster(pipe, width, height, provider.extent(), crs) is deprecated and You should use file_writer.writeRaster(pipe, width, height, provider.extent(), crs, transformContext) Here is more info about it.

Lastly, I would suggest taking a look on this discussed topics that may help to solve your problem:

Using exitQgis() in PyQGIS?
https://github.com/qgis/QGIS/issues/35354

Matt
  • 16,843
  • 3
  • 21
  • 52