Sum the values of a single raster and output the value as a number?

Question

I have a floating point raster and have applied a conditional to it in order to only extract the positive values, so there are 'no data' "holes" in the raster. I simply need to sum the values from this raster and output it as a number.

I came across this thread, which I'm sure some of you will recognize. Has anyone developed a working solution? I'm working in Arc10 and making use of ArcObjects if necessary.

EDIT: I've included a VB.NET version of @Kirk's C# code. It's provided as is, but it should be a good starting point if it doesn't already completely satisfy its purpose.

    Public Shared Sub TestSumRaster()
        Dim path As String = "D:\topo_modeling\topo_modeling\gis\usgs\ned_10m.img"
        Dim rds = OpenRaster(path)
        Dim sum As Double = SumRaster(rds)
        Debug.Print("sum : {0}", sum)
    End Sub

    Public Shared Function SumRaster(rasDs As IRasterDataset) As Double
        ' sum up each band in a raster dataset
        Dim rasBandCol As IRasterBandCollection = DirectCast(rasDs, IRasterBandCollection)
        Dim sum As Double = 0.0

        ' Iterate through each band of the dataset.
        For i As Integer = 0 To rasBandCol.Count - 1
            sum += SumBand(rasBandCol.Item(i))
        Next
        Return sum
    End Function

    Private Shared Function SumBand(rBand As IRasterBand) As Double

' sum up each block in a band
        Dim sum As Double = 0.0
        ' QI to IRawBlocks from IRasterBandCollection
        Dim rawBlocks = DirectCast(rBand, IRawBlocks)

        ' Determine the tiling scheme for the raster dataset.
        Dim tile = New Tile(rawBlocks.RasterInfo)

        ' NoData is array, but usually just has one value in it.
        Debug.Print(rawBlocks.RasterInfo.NoData.[GetType]().ToString())

        ' Create the pixel block.
        Dim pb = rawBlocks.CreatePixelBlock()

        ' Iterate through the pixel blocks.
        For pbYcursor As Integer = tile.startY To tile.endY - 1
            For pbXcursor As Integer = tile.startX To tile.endX - 1
                ' Get the pixel block.
                rawBlocks.ReadBlock(pbXcursor, pbYcursor, 0, pb)
                sum += SumBlock(pb, pbXcursor, pbYcursor, rawBlocks.RasterInfo.NoData)
            Next
        Next
        Return sum
    End Function

    Private Shared Function SumBlock(pb As IPixelBlock, pbXcursor As Integer, pbYcursor As Integer, noDataValues As Object) As Double
        ' sum up each pixel in the pixelblock
        Dim sum As Double = 0.0

        ' Put the pixel block into a SafeArray for manipulation.
        Dim safeArray = TryCast(DirectCast(pb.get_SafeArray(0), System.Array), System.Array)

        Dim noDataList = New List(Of [Single])(DirectCast(noDataValues, [Single]()))

        ' Iterate through the pixels in the pixel block.
        For row As Integer = 0 To pb.Height - 1
            For col As Integer = 0 To pb.Width - 1
                Dim val As Object = safeArray.GetValue(row, col)
                Dim singVal = DirectCast(val, [Single])

                'Debug.Print(val.GetType().Name);
                If Not noDataList.Contains(singVal) Then
                    sum += singVal
                End If
            Next
        Next
        Return sum
    End Function

Answer 1

Could you experiment with a known grid...here is a suggestion for implementation using Python, numpy and arcpy.

>>> import arcpy
>>> npArray = arcpy.RasterToNumPyArray("c:/temp/x")
>>> import numpy
>>> print numpy.sum(npArray)
45763093
>>> 

EDIT Additional comments on the above.

You have to be careful when dealing with grids which contain nodata. Consider a small grid with the following values:

1 2 3

-- 5 --

7 8 --

where -- represents nodata. For example, if the above grid is converted to a numpy array, then summed:

>>> z = arcpy.RasterToNumPyArray("c:/temp/small_grid")
>>> np.sum(z)
nan

the result is obviously unexpected.
If array is converted to a masked array, then the summation is correct.

>>> z1 = np.ma.masked_array(z, np.isnan(z))
>>> np.sum(z1)
26.0

Answer 2

Create a single value raster (CreateConstantRaster).

Calculate the zonal sum (ZonalStatisticsAsTable).

Read "SUM" column of resulting single row table.

Answer 3

Here is code adapted from this sample. You'll need to modify it if your raster is double (instead of Single). Tested with a DEM downloaded from here. I heard somewhere that C# is now better at handling contravariance. If so, I'd be curious how this code could be generalized to support either single or double precision rasters.

public static void TestSumRaster()
{
    string path = @"D:\topo_modeling\topo_modeling\gis\usgs\ned_10m.img";
    var rds = OpenRaster(path);
    double sum = SumRaster(rds);
    Debug.Print("sum : {0}", sum);
}

public static double SumRaster(IRasterDataset rasDs)
{
    // sum up each band in a raster dataset
    IRasterBandCollection rasBandCol = (IRasterBandCollection)rasDs;
    double sum = 0.0;

    // (can't think of a case where you'd need to sum across more than one band)

    // Iterate through each band of the dataset.
    for (int i = 0; i < rasBandCol.Count; i++)
    {
        sum += SumBand(rasBandCol.Item(i));
    }
    return sum;
}

private static double SumBand(IRasterBand rBand)
{
    // sum up each block in a band
    double sum = 0.0;
    // QI to IRawBlocks from IRasterBandCollection.
    var rawBlocks = (IRawBlocks)rBand;

    // Determine the tiling scheme for the raster dataset.
    var tile = new Tile(rawBlocks.RasterInfo);

    // NoData is array, but usually just has one value in it.
    Debug.Print(rawBlocks.RasterInfo.NoData.GetType().ToString());

    // Create the pixel block.
    var pb = rawBlocks.CreatePixelBlock();

    // Iterate through the pixel blocks.
    for (int pbYcursor = tile.startY; pbYcursor < tile.endY; pbYcursor++)
    {
        for (int pbXcursor = tile.startX; pbXcursor < tile.endX; pbXcursor++)
        {
            // Get the pixel block.
            rawBlocks.ReadBlock(pbXcursor, pbYcursor, 0, pb);
            sum += SumBlock(pb, pbXcursor, pbYcursor,rawBlocks.RasterInfo.NoData);
        }
    }
    return sum;
}

private static double SumBlock(IPixelBlock pb, int pbXcursor, int pbYcursor, object noDataValues)
{
    // sum up each pixel in the pixelblock
    double sum = 0.0;

    // Put the pixel block into a SafeArray for manipulation.
    var safeArray = (System.Array)pb.get_SafeArray(0) as System.Array;

    var noDataList = new List<Single>((Single[])noDataValues);

    // Iterate through the pixels in the pixel block.
    for (int row = 0; row < pb.Height;row++)
    {
        for (int col = 0; col < pb.Width;col++)
        {
            object val =safeArray.GetValue(row, col);
            var singVal = (Single) val;

            //Debug.Print(val.GetType().Name);
            if(!noDataList.Contains(singVal))
                sum += singVal;
        }
    }
    return sum;
}

public static IRasterDataset OpenRaster(string path)
{
    IGPUtilities3 gpUtil = new GPUtilitiesClass();
    var rds = gpUtil.OpenFromString(path) as IRasterDataset;
    return rds;
}
public class Tile
{
    public int startX;
    public int endX;
    public int startY;
    public int endY;
    public Tile(IRasterInfo rasInfo)
    {
        this.startX = (int)Math.Floor((rasInfo.Extent.Envelope.XMin - rasInfo.Origin.X)
            / (rasInfo.BlockWidth * rasInfo.CellSize.X));
        this.endX = (int)Math.Ceiling((rasInfo.Extent.Envelope.XMax - rasInfo.Origin.X)
            / (rasInfo.BlockWidth * rasInfo.CellSize.X));

        this.startY = (int)Math.Floor((rasInfo.Origin.Y - rasInfo.Extent.Envelope.YMax)
            / (rasInfo.BlockHeight * rasInfo.CellSize.Y));
        this.endY = (int)Math.Ceiling((rasInfo.Origin.Y - rasInfo.Extent.Envelope.YMin)
            / (rasInfo.BlockHeight * rasInfo.CellSize.Y));
    }
}

Answer 4

I have worked out a solution using standard geoprocessing that can be implemented using ArcObjects.

Input: Floating-point raster (can include no-data holes)

1) Convert the input raster to integer

2) Convert the integer raster to a polygon

3) Perform Zonal Statistics on the input raster (using SUM or ALL)

4) Using ArcObjects read the SUM column and add all the value to get the raster sum

You can store this value as a double. In my case I was calculating raster value so I simply multiplied the sum of the raster values by the raster cell resolution.

Thanks for all your input - it really helped me to figure out this task.