# How can I overlay two raster without turning them in vector format in Python to get zone wise statistics?

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.

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

#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')
``````
• 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 Jan 29 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 Jan 29 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 Jan 29 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 Jan 29 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 Jan 29 at 16:56