8

I would like to reclassify a raster file from a raster with 10 classes to a raster with 8 classes using pyhton, gdal and/or numpy. The classes are represented as integers. I have tried following the steps from this post Reclassify rasters using GDAL and Python , the numpy.equal doc and also gdal_calc doc. However, to no avail.

The raster file to be reclassified has integer values ranging from 0 to 11 and also include values 100 and 255. The following show the reclass (from value : to value):

nodata : 4, 0 : 4, 1 : 1, 2 : 2, 3 : 3, 4 : 3, 5 : 4, 6 : 5, 7 : 5, 8 : 6, 9 : 7, 10 : 8, 100 : nodata, 255 : nodata,

What I have been able to do is select the raster file to be reclassified using tkinter.FileDialog and get the raster info such as geotransform, and pixel size with reclass = gdal.Open(raster, GA_ReadOnly).

How do I go about solving the above.

It might be worth mentioning that the rasters to be reclassified can be fairly large in some cases (500mb to 5gb).

  • There is another example on the GeoExamples Blog – bennos Sep 16 '15 at 7:53
  • @bennos, tried the script on the blog but it returns a memory error when unpacking the array. – PyMapr Sep 17 '15 at 8:25
  • I suggest you discuss this problem with Roger Veciana i Rovira, the author of the post, as he knows his code better than I do and maybe knows how to resolve the issue – bennos Sep 17 '15 at 8:36
  • Changing the input raster from 16Bit unsigned to 8Bit unsigned solved the memory issue. However, it takes about the same amount of time to reclassify as the dmh126 script below. – PyMapr Sep 17 '15 at 9:52
6

Here you have a simple python script for reclassification, I wrote it and it works for me:

from osgeo import gdal

driver = gdal.GetDriverByName('GTiff')
file = gdal.Open('/home/user/workspace/raster.tif')
band = file.GetRasterBand(1)
lista = band.ReadAsArray()

# reclassification
for j in  range(file.RasterXSize):
    for i in  range(file.RasterYSize):
        if lista[i,j] < 200:
            lista[i,j] = 1
        elif 200 < lista[i,j] < 400:
            lista[i,j] = 2
        elif 400 < lista[i,j] < 600:
            lista[i,j] = 3
        elif 600 < lista[i,j] < 800:
            lista[i,j] = 4
        else:
            lista[i,j] = 5

# create new file
file2 = driver.Create( 'raster2.tif', file.RasterXSize , file.RasterYSize , 1)
file2.GetRasterBand(1).WriteArray(lista)

# spatial ref system
proj = file.GetProjection()
georef = file.GetGeoTransform()
file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.FlushCache()

Just change the ranges.

I hope it will help.

  • 3
    You should close file2 either with del file2 or file2 = None to make sure it gets written to disk. .FlushCache() only influences GDALs internal tile cache. – Kersten Sep 15 '15 at 13:30
  • @dmh126, I modified the ranges and the script works. However, it is not very quick (quick being debatable). The input raster was about 120mb and took about 15 min to complete.With the aid of a commercial remote sensing package it takes seconds. Any recommendations on decreasing processing time? – PyMapr Sep 17 '15 at 8:23
  • I think that multithreading can help. You can try to use all of your cores, there is a question gis.stackexchange.com/questions/162978/… – dmh126 Sep 17 '15 at 9:07
  • Doesn't make sense to use a double for loop, see answer below – Mattijn Jan 21 '16 at 8:43
  • Right, the double loop and per-element reclassification is the slowest of all possible ways to do this. Use the powerful parts of numpy like ufuncs: docs.scipy.org/doc/numpy-1.10.1/reference/ufuncs.html. – sgillies May 27 '16 at 15:57
15

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where:

# reclassification    
lista[np.where( lista < 200 )] = 1
lista[np.where((200 < lista) & (lista < 400)) ] = 2
lista[np.where((400 < lista) & (lista < 600)) ] = 3
lista[np.where((600 < lista) & (lista < 800)) ] = 4
lista[np.where( lista > 800 )] = 5

On an array of 6163 by 3537 pixels (41.6mb) the classification is done in 1.59 seconds, where it takes 12min 41s using the double for loop. In total just a speedup of 478x.

Bottomline, never use a double for loop using numpy

  • Thanks for the hint, but I think that will lead into a problem if the input classes overlap with the output classes. I don't want my new value to be changed by the next rule. – etrimaille May 11 '16 at 6:20
  • @Gustry - Just ran into that problem here. – relima Aug 29 '16 at 23:20
  • So check my answer below : gis.stackexchange.com/questions/163007/… – etrimaille Aug 31 '16 at 6:50
6

Here's a basic example using rasterio and numpy:

import rasterio as rio
import numpy as np


with rio.open('~/rasterio/tests/data/rgb1.tif') as src:
    # Read the raster into a (rows, cols, depth) array,
    # dstack this into a (depth, rows, cols) array,
    # the sum along the last axis (~= grayscale)
    grey = np.mean(np.dstack(src.read()), axis=2)

    # Read the file profile
    srcprof = src.profile.copy()

classes = 5
# Breaks is an array of the class breaks: [   0.   51.  102.  153.  204.]
breaks = (np.arange(classes) / float(classes)) * grey.max()

# classify the raster
classified = np.sum(np.dstack([(grey < b) for b in breaks]), axis=2).reshape(1, 400, 400).astype(np.int32)

# Update the file opts to one band
srcprof.update(count=1, nodata=None, dtype=classified.dtype)

with rio.open('/tmp/output.tif', 'w', **srcprof) as dst:
    # Write the output
    dst.write(classified)
2

Just to complete the answer from @Mattijn, I think that will lead into a problem if the input classes overlap with the output classes. I don't want my new value to be changed by the next rule.

I don't know if I loose speed but I should do a deep copy :

list_dest = lista.copy()

list_dest[np.where( lista < 0 )] = 0
list_dest[np.where((0 <= lista) & (lista <= 1)) ] = 1
list_dest[np.where((1 < lista) & (lista <= 5)) ] = 2
list_dest[np.where( 5 < lista )] = 3
1

Here is another Rasterio approach that I hacked together using the Rasterio Cookbook and @Mattijn's answer.

import rasterio
import numpy as np

with rasterio.open('input_raster.tif') as src:    
    # Read as numpy array
    array = src.read()
    profile = src.profile

    # Reclassify
    array[np.where(array == 0)] = 4 
    array[np.where(array == 2)] = 1
    # and so on ...  

with rasterio.open('output_raster.tif', 'w', **profile) as dst:
    # Write to disk
    dst.write(array)
0

In some cases, numpy digitize can be useful to quickly go from ranges to bins.

import rasterio
import numpy as np

with rasterio.open('my_raster.tif') as src:    
    array = src.read()
    profile = src.profile
    bins = np.array([-1.,-0.7,-0.4, 0.2, 1]) 
    inds = np.digitize(array, bins)

with rasterio.open('output_raster.tif', 'w', **profile) as dst:
    dst.write(inds)

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.