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

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.
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

-

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.
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));
}
}

-

Create a single value raster (CreateConstantRaster).

Calculate the zonal sum (ZonalStatisticsAsTable).

Read "SUM" column of resulting single row table.

-

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
>>>


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

-
Doesn't ArcGIS 10 now expose the grid statistics? That would likely be the most efficient way to read the sum out of an entire grid. –  whuber Jul 9 '11 at 0:08
Mean and Standard Deviation, min ana max are reported only. –  Dan Patterson Jul 9 '11 at 0:13
Mean (times) #rows (times) #columns = sum. –  whuber Jul 9 '11 at 11:26
In my example above, grid statistics yielded mean = 4.3333 (sum 26 / 6 data cells) but the sum would be incorrect if multiplied by #rows and # cols. –  Dan Patterson Jul 9 '11 at 13:47
That's a good point, but there's a simple workaround: divide the apparent sum by the mean of the data-indicator grid. In your example the data-indicator grid has 6 ones and 3 zeros for a mean of 2/3. Multiplying 13/3 by 9 cells gives 39; adjusting by 2/3 reduces this to 26, the correct answer. –  whuber Jul 11 '11 at 13:48

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

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

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.

-
Sure seems like a lot of disk i/o. –  Kirk Kuykendall Jul 12 '11 at 20:48
@Kirk - You're right - it is. It's by no means perfect, but it does the job in a straight-forward fashion (relying on geoprocessing). I'm bound by some requirements on this task, but I will take a shot at converting your above code sample to VB (also task-bound) and see if it can cut the processing time along with the i/o. –  Radar Jul 12 '11 at 20:52
@Radar - Why are you converting to polygon in step 2? ArcGIS will convert it back to a raster behind the scenes anyway... BTW - my answer above does exactly what you want using standard geoprocessing and is simpler than your answer as you don't have to to add multiple values from the zonal stats SUM column. With my method, you simply read the first (and only) row. –  Luke Jul 13 '11 at 5:02