[GIS] Reclassify a raster value to -9999 and set it to the nodata value using python and or gdal

gdalnumpypython-2.7rasterreclassify

I modified a script which reclassifies a raster. It is based on the script by Roger in his blog (the second script) and it works pretty well.The problem I am facing is as follows:
I would like to reclassify a value (in this case 100) to -9999 and then set -9999 as the nodata value. Below is the list of values to classify and the list to which the corresponding value should be reclassified. Note that value 100 should go to value -9999.

classification_values = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,100] #The interval values to classify
classification_output_values = [4,1,2,3,4,5,5,6,7,8,9,11,12,13,14,15,16,17,18,10,14,14,5,16,-9999] #The value assigned to each interval

However, when I view the resultant raster after running the reclass script, value 100 is not -9999, it defaults to 0, but -9999 is set as the nodata value.

In the script I added after line 33 in Roger's script

dst_ds.GetRasterBand(1).SetNoDataValue(-9999)

If I reclass value 100 to 0 and set 0 as the nodata value the resultant raster is correct and everywhere were value 100 used to be there is nodata.

Not sure if the followong info might be of use:
Input raster is 8bit unsigned (have tried 16bit signed with no luck);
I have changed the numpy array from 8bit to 16bit signed (line 47 in Roger's script);
Output raster is 16bit signed.

How would I go about solving this problem?

Best Answer

If you are looking for an easy way to change specific pixels to NoData (Nan) or -9999 please take a look at this alternative script I wrote: 

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 12 21:33:20 2015

@author: rutgerhofste

This script will replace values in a raster and save the result in geotiff format

"""

from osgeo import gdal 
import numpy as np

firstrun = 1


def readFile(filename):
    filehandle = gdal.Open(filename)
    band1 = filehandle.GetRasterBand(1)
    geotransform = filehandle.GetGeoTransform()
    geoproj = filehandle.GetProjection()
    Z = band1.ReadAsArray()
    xsize = filehandle.RasterXSize
    ysize = filehandle.RasterYSize
    return xsize,ysize,geotransform,geoproj,Z

def writeFile(filename,geotransform,geoprojection,data):
    (x,y) = data.shape
    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    # you can change the dataformat but be sure to be able to store negative values including -9999
    dst_datatype = gdal.GDT_Float32
    dst_ds = driver.Create(filename,y,x,1,dst_datatype)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(geoprojection)
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    return 1



pathname = "/yourpathhere/filename.tif"
writefilename = "/youroutputpathhere/filename.tif"


if firstrun == 1 :
    [xsize,ysize,geotransform,geoproj,Z] = readFile(pathname)


# Set large negative values to -9999
Z[Z<-9999]= -9999
# Choose your preference: (comment either rule)
Z[Z==-9999]= np.nan
# Or 
Z[np.isnan(Z)]= -9999



writeFile(writefilename,geotransform,geoproj,Z)

# Open your file in QGIS and set Nodata value if necessary for colormaps etc.

The current output datatype is set to float32.

Related Question