1

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 "<input>", line 1, in <module>
  File "/home/nikan/Untitled-2.py", line 35, in <module>
    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.
2

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:

enter image description here

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

enter image description here

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.