1

I'm pretty new to QGIS and Python but I've been making a number of processing scripts that automate GIS processes for my job. I'm working on one now that requires a raster reclassification step - the raster was produced by interpolating a point shapefile based on elevation data and then clipping it using a polygon shapefile. I want to make this code universally applicable for the data I will be working with so I wrote something that takes the range of the pixels and breaks them up into even steps to create 7 logical classes for whatever data I may be working with. My code snippet to perform the reclass is below:

        #Now we must reclassify the raster
        
        temporary_output = 'C:/temp/QGIS_Temp/ReclassRas.sdat'
        
        raster_provider = clippedLyr.dataProvider()
        stats = raster_provider.bandStatistics(1, QgsRasterBandStats.All)
        range = stats.maximumValue - stats.minimumValue
        step = ceil((range/7)*100)/100
        pos_one = floor(stats.minimumValue)
        pos_two = ceil(stats.minimumValue + step)
        pos_three = ceil(pos_two + step)
        pos_four = ceil(pos_three + step)
        pos_five = ceil(pos_four + step)
        pos_six = ceil(pos_five + step) 
        pos_seven = ceil(pos_six + step)
        pos_eight = ceil(pos_seven + step)
        
        table_list = [pos_one,pos_two,1,pos_two,pos_three,2,pos_three,pos_four,3,pos_four,pos_five,4,pos_five,pos_six,5,pos_six,pos_seven,6,pos_seven,pos_eight,7]
                
        alg_params = {
            'INPUT_RASTER': clippedLyr,
            'RASTER_BAND': 1,
            'TABLE': table_list,
            'NO_DATA':-9999,
            'RANGE_BOUNDARIES':0,
            'NODATA_FOR_MISSING':False,
            'DATA_TYPE':5,
            'OUTPUT': temporary_output
        }
        
        processing.run("native:reclassifybytable", alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        
        reclassLyr = QgsRasterLayer(temporary_output)
        reclassLyr.isValid()
        
        QgsProject.instance().addMapLayer(reclassLyr)

When I print the table_list variable I made using the PyConsole, I get the correct values in a list. When I use these values in the reclassify by table GUI, I get the correct reclassified output raster that I want. But when I try to run my code, the output raster has only two visible classes: white (7) and black (-9999) and it also extends beyond my clipped raster boundary for some reason.

Has anyone been able to correctly execute the reclassify by table processing tool in a python script?

1 Answer 1

1

I found a way to accomplish what I wanted using the saga tool: reclassify values (simple). I can use the exact same list I made for the lookup table and everything works fine. I'm not sure if there is an advantage to using the native tool: reclassify by table, but at least for my purposes the saga tool is great.It's possible that the saga tool works better because it is more specific to grid files (.sdat) but I'm not sure.

        temporary_output = 'C:/temp/QGIS_Temp/ReclassRas_{}.sdat'.format(time_stamp)
        
        raster_provider = clippedLyr.dataProvider()
        stats = raster_provider.bandStatistics(1, QgsRasterBandStats.All)
        range = stats.maximumValue - stats.minimumValue
        step = ceil((range/7)*100)/100
        pos_one = floor(stats.minimumValue)
        pos_two = (stats.minimumValue + step)
        pos_three = (pos_two + step)
        pos_four = (pos_three + step)
        pos_five = (pos_four + step)
        pos_six = (pos_five + step) 
        pos_seven = (pos_six + step)
        pos_eight = ceil(pos_seven + step)
        
        table_list = [pos_one,pos_two,1,pos_two,pos_three,2,pos_three,pos_four,3,pos_four,pos_five,4,pos_five,pos_six,5,pos_six,pos_seven,6,pos_seven,pos_eight,7]
        
        alg_params_3 = {
            'GRID_IN': clippedLyr,
            'METHOD': 2,
            'LOOKUP': table_list,
            'GRID_OUT': temporary_output,
        }
        
        processing.run("saga:reclassifyvaluessimple", alg_params_3, context=context, feedback=feedback, is_child_algorithm=True)
        
        reclassLyr = QgsRasterLayer(temporary_output)
        reclassLyr.isValid()
        
        QgsProject.instance().addMapLayer(reclassLyr)


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.