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)
array([[-2147483648,6645],
[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')
