4

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?

  • Did you by chance change line 47 to r = zeros((rows, cols), numpy.uint16)? That would be an unsigned integer which would explain why any negative value becomes 0 instead. If that is the case try numpy.int16. – Kersten Sep 21 '15 at 13:12
  • @Kersten, had a quick look at the script to make sure and it is numpy.int16. – PyMapr Sep 21 '15 at 13:24
4

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.

  • 1
    The script works great thus far, thanks. I made some small changes to the script to suit my specific requirements. At the bottom of the writeFile definition I added dst_ds.GetRasterBand(1).SetNoDataValue(-9999). Which removes the extra step of setting the nodata value in QGIS. Also, I changed Z[Z<-9999]= -9999 to Z[Z==100]= -9999. – PyMapr Sep 22 '15 at 7:49
  • thats a great addition actually. much appreciated – RutgerH Sep 22 '15 at 12:38

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.