I have two raster re-sampled and same extent one is elevation class raster from DEM reclassification classified in 10 zone with 500m interval and another is snow cover fraction raster has snow fraction value 0-100 and some reassigned classes like 201 for cloud 205 for water etc.

I want to have a table of elevation class-wise snow cover fraction. How to calculate that in Python? How to sort the array and make query?

I do not want to convert the input raster format to vector format because in the later stage I need the output file in raster format actually the array to read in to build contingency matrix.

I have tried with from rasterstats import zonal_stats but it needs to have one file in vector format however I have two raster data.

1 Answer 1


Solution part 1 with gdal tool

  • Multiply your 10 DEM class raster (i.e from 1 to 10) by 1000 (i.e from 1000 to 10000).

gdal_calc.py -A dem_reclass.tif --calc=A*1000 --outfile=dem_reclass_1000.tif

  • Add this to your 100 snow class raster (i.e from 0 to 100).

gdal_calc.py -A dem_reclass_1000.tif -B snow_reclass.tif --outfile=composite.tif --calc="A+B"

You will obtain a raster with values from 0 to 10100, where the thousand number is the dem class and the units number < 101 is the snow class.

Example 1 : for a pixel where dem class is 5 and snow class is 45 the final value will be 5045. Example 2 : for a pixel where dem class is 10 and snow class is 201 the final value will be 10201.

Now you can use some count method to see the pixel values :

Solution part 2 with osgeo tools

from osgeo import gdal_array
import numpy as np

# Read raster data as numeric array from file
rasterArray = gdal_array.LoadFile('composite.tif')

#convert the array to integer
x = rasterArray.astype(int)

#create an array with the pixel value and the number of pixel
np.array(np.unique(x, return_counts=True)).T

I get a 64bit data type array like below on a test monoband image where it appears there's 3355 pixels of value = 1000 and 6645 pixels of nodata (-2147483648)

[1000,3355]], dtype=int64)

To get the data, i use the following

result = np.array(np.unique(x, return_counts=True)).T
result[0,0] # returns -2147483648 (no data)
result[0,1] # returns 6645 (nb of pixels of nodata)
result[1,0] # returns 1000 (value 1000)
result[1,1] # returns 3355 (nb of pixels of value 1000)

you should be able to export the array with

np.savetxt('result_file.csv', result, delimiter=',', fmt='%s')
  • is it necessary to multiply the DEM with 1000 ? I have tried to run the 1st command but got the error ` 1000' does not exist in the file system, and is not recognised as a supported dataset name.
    – Tua
    Commented Jan 29, 2020 at 15:06
  • I modified the formula. The B was for a dataset name indeed and i mistook x instead of * for multiplication. Sorry for the error. The 1000 is meant to push your dem classes in the thousands position to be able to find them in the resultant calculated raster
    – gisnside
    Commented Jan 29, 2020 at 15:13
  • I have simply tried with 2nd command it worked as you suggested however I have huge number of every day raster which I can do like this by creating a for loop but how to get the table statistics like in 0-500m (class-1) has snow cover count (pixel) from counter??
    – Tua
    Commented Jan 29, 2020 at 15:23
  • My example is from another post, as you have 2 questions in one. I need to check the second one now :)
    – gisnside
    Commented Jan 29, 2020 at 16:33
  • I changed the code, added an import of numpy, removed Counter couldn't make it work properly + had a try on a monoband TIF, it seems to work. Please tell me if it works for you. The result is an array you should be able to use.
    – gisnside
    Commented Jan 29, 2020 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.