Simple calculation for multiple rasters

[StephenFricke]

Hello, I am working on a project that chronicles the snow covered area in the Clearwater Basin in Idaho.  I have hundreds of rasters which all have pixels that are classified as either snow or no snow.  I want to calculate the snow covered area for each raster, which is the number of pixels in the raster classified as snow divided by the total number of pixels in the raster.  Obviously it would be a huge pain to have to do this manually for each raster so I was hoping that there would be a way to write a program to do this.  Let me know if you have any suggestions.  Thanks very much!

[JakeSkinner]

Hi Peter,

Are you rasters classified by pixel value (i.e. value > 3 equals snow covered) or by an actual field attribute?  What pixel type are your rasters (i.e. floating point, signed integer)?

Below is an example on how you can do this.  In the example, the rasters I am working with have a pixel type of floating point.  What the script does is convert the floating point rasters to integer by multiplying the raster by the number of decimal places I would like to preserve.  For example, I multiply the rasters by 100 to maintain 2 decimal values.  Next I add an attribute table to the raster and calculate the total number of pixels.

Next I want to find the percentage of pixels with a greater value of 3.25 (which represents snow cover).  I multiply this value by 100, sum all the pixels with a value greater than this and divide by the total number of pixels to find the percentage.

import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

env.workspace = r"C:\temp\python\rasters"
env.overwriteOutput = True

lstRasters = arcpy.ListRasters("*")
for raster in lstRasters:
    outTimes = Times(raster, 100)
    outTimes.save(r"C:\temp\python\outTimes.tif")
    intRas = Int(r"C:\temp\python\outTimes.tif")
    intRas.save(r"C:\temp\python\intRas.tif")
    arcpy.BuildRasterAttributeTable_management(r"C:\temp\python\intRas.tif", "Overwrite")
    desc = arcpy.Describe(r"C:\temp\python\intRas.tif")
    totalpixels = desc.height * desc.width

    list = []
    pixValue = 3.25 * 100

    rows = arcpy.SearchCursor(intRas, "VALUE > " + str(pixValue))
    for row in rows:
        list.append(row.Count)

    print float(sum(list)) / float(totalpixels)

del row, rows
arcpy.Delete_management(r"C:\temp\python\intRas.tif")
arcpy.Delete_management(r"C:\temp\python\outTimes.tif")

[StephenFricke]

JSkinn-
Thank you so much!  That script is very helpful.  The only problem is I am not getting the total pixels correctly.  Each raster has an attribute table that looks like this....

Rowid VALUE COUNT
0        1        16456
1        2        50594


The values for no snow are 1 and snow is 2.  I am able to get the sum(list) of snow values correctly, but unable to calculate the total pixels correctly (which is simply the two values in the count field added together).  How would I make the "totalpixels" value equal to the sum of the COUNT field.  Thanks again!

[JakeSkinner]

Hi Peter,

It appears that if the raster has any rotation, the columns and rows will include the 'NoDATA' areas which will give you more pixels when multiplying the columns and rows together using the Describe function.  You can sum the COUNT field using a Search Cursor and appending the values to another list.  The below should work:

import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

env.workspace = r"C:\temp\python\rasters"
env.overwriteOutput = True

lstRasters = arcpy.ListRasters("*")
for raster in lstRasters:
    outTimes = Times(raster, 100)
    outTimes.save(r"C:\temp\python\outTimes.tif")
    intRas = Int(r"C:\temp\python\outTimes.tif")
    intRas.save(r"C:\temp\python\intRas.tif")
    arcpy.BuildRasterAttributeTable_management(r"C:\temp\python\intRas.tif", "Overwrite")
    list2 = []
    rows = arcpy.SearchCursor(r"C:\temp\python\intRas.tif")
    for row in rows:
        list2.append(row.getValue("COUNT"))
    totalpixels = sum(list2)

    list = []
    pixValue = 3.25 * 100

    rows = arcpy.SearchCursor(intRas, "VALUE > " + str(pixValue))
    for row in rows:
        list.append(row.Count)

    print float(sum(list)) / float(totalpixels)

del row, rows
arcpy.Delete_management(r"C:\temp\python\intRas.tif")
arcpy.Delete_management(r"C:\temp\python\outTimes.tif")

[StephenFricke]

Thanks alot JSkinn, I really appreciate it!