Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
add comment

4 Answers 4

up vote 1 down vote accepted

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));
    }
}
share|improve this answer
add comment

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.

share|improve this answer
    
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. –  Nick Ochoski 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
add comment

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
share|improve this answer
    
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
add comment

Create a single value raster (CreateConstantRaster).

Calculate the zonal sum (ZonalStatisticsAsTable).

Read "SUM" column of resulting single row table.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.