pyqgis - How to reclassify a raster based on quantile and without no data in QGIS?

I want to reclassify quantile a raster layer without negative values (nodata values) so I used Reclassify a raster file with quantiles and How to filter no data value with GDAL? in below code:

import numpy as np
from osgeo import gdal, gdal_array


# open the dataset and retrieve raster data as an array
dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()
lst = []
for v in array:
    for element in v:
        if element > 0:
            lst.append(element)

# create an array of zeros the same shape as the input array

output = np.zeros_like(lst).astype(np.uint8)

# use the numpy percentile function to calculate percentile thresholds
percentile_80 = np.percentile(lst, 80)
print percentile_80
percentile_60 = np.percentile(lst, 60)
print percentile_60
percentile_40 = np.percentile(lst, 40)
print percentile_40
percentile_20 = np.percentile(lst, 20)

print percentile_20
percentile_0 = np.percentile(lst, 0)
print percentile_0

output = np.where((lst > percentile_0), 1, output)
output = np.where((lst > percentile_20), 2, output)
output = np.where((lst > percentile_40), 3, output)
output = np.where((lst > percentile_60), 4, output)
output = np.where((lst > percentile_80), 5, output)


outname = "/path/to/newfile.tif"
gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)

But I received below error:

Traceback (most recent call last):
  File "", line 1, in 
  File "/home/nikan/Untitled-2.py", line 35, in 
    gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal_array.py", line 239, in SaveArray
    return driver.CreateCopy( filename, OpenArray(src_array,prototype) )

  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 1538, in CreateCopy
    return _gdal.Driver_CreateCopy(self, *args, **kwargs)
ValueError: Received a NULL pointer.

Answer

I tried out approach in Reclassify a raster file with quantiles and it produces bad results because it is necessary to change values in a loop to avoid self reference (produced with 'where' numpy method). Instead, you can use following code.

import numpy as np
from osgeo import gdal, gdal_array, osr

dataset = gdal.Open("/home/zeito/pyqgis_data/test_raster_nodata.tif")
band = dataset.GetRasterBand(1)

nodata = band.GetNoDataValue()
array = dataset.ReadAsArray()

new_array = array
nan_array = array

nan_array[array == nodata] = np.nan

percentile_80 = np.nanpercentile(nan_array, 80)
percentile_60 = np.nanpercentile(nan_array, 60)

percentile_40 = np.nanpercentile(nan_array, 40)
percentile_20 = np.nanpercentile(nan_array, 20)
percentile_0 = np.nanpercentile(nan_array, 0)

for i, v in enumerate(new_array):
    for j, element in enumerate(v):
        if element <= percentile_0:
            new_array[i,j] = 1
        if element > percentile_0 and element <= percentile_20:
            new_array[i,j] = 2

        if element > percentile_20 and element <= percentile_40:
            new_array[i,j] = 3
        if element > percentile_40 and element <= percentile_60:
            new_array[i,j] = 4
        if element > percentile_60 and element <= percentile_80:
            new_array[i,j] = 5
        if element > percentile_80:
            new_array[i,j] = 6

new_array[ new_array != new_array ] = nodata


geotransform = dataset.GetGeoTransform()
wkt = dataset.GetProjection()

# Create gtif file
driver = gdal.GetDriverByName("GTiff")
output_file = "/home/zeito/pyqgis_data/test_raster_nodata_reclass.tif"

dst_ds = driver.Create(output_file,
                       band.XSize,

                       band.YSize,
                       1,
                       gdal.GDT_Int16)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(nodata)
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution

dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

dataset = None
dst_ds = None

After running above code with raster of following image:

I got reclassified version below. Value Tool QGIS plugin helps me to corroborate that reclassification was as expected.