1

I am using c# GDAL bindings in attempts to read raster values and preform fairly simple calculations. In summary, I would like to get the mean of the values and then get the sum of values below the mean and the sum of values above the mean. Everything appears to work fine except for when I attempt to read the raster to a buffer. The result is producing an array in which every value is -3.402823E+38

The data is float data and i wish to keep that way for my calculations. Code is as follows.

class Program
{
    static void Main(string[] args)
    {
        Gdal.AllRegister();
        Ogr.RegisterAll();

        Dataset dataset = Gdal.Open(@"xxxxxxxxxxx.tif", Access.GA_ReadOnly);
        Band band = dataset.GetRasterBand(1);
        int width = band.XSize;
        int height = band.YSize;
        int size = width * height;
        double min = 0.00;
        double max = 0.00;
        double mean = 0.00;
        double stddev = 0.00;

        var stats = band.GetStatistics(1, 0, out min, out max, out mean, out stddev);

        //Console.WriteLine($"Statistics retrieved and returned a result of {stats}");
        Console.WriteLine($"X : {width} Y : {height} SIZE: {size}");
        Console.WriteLine($"MIN : {min} MAX : {max} MEAN : {mean} STDDEV : {stddev}");
        DataType type = band.DataType;
        Console.WriteLine($"Data Type : {type}");

        float gtMean = 0; //cut
        float ltMean = 0; //fill

        float[] data = new float[size];
        var dataArr = band.ReadRaster(0, 0, width, height, data, width, height, 0, 0);

        int i, j;
        for (i = 0; i < width; i++)
        {
            for (j = 0; j < height; j++)
            {
                float value = data[i + j * width];
                if (value > (float)mean)
                {
                    gtMean += value;
                }

                if (value < (float)mean)
                {
                    ltMean += value;
                }
            }
        }


        Console.WriteLine($"Sum of values above the mean {gtMean}");
        Console.WriteLine($"Sum of values below the mean {ltMean}");

        double pixelArea = 0.10763911106;

        Console.WriteLine("Press Any Key To Exit....");
        Console.ReadLine();
    }
}

enter image description here

EDIT / UPDATE As suggested in the comments the values i was observing represent the no data values. The below script accounts for these values and excludes them from the sums.

class Program
{
    static void Main(string[] args)
    {
        Gdal.AllRegister();
        Ogr.RegisterAll();

        Dataset dataset = Gdal.Open(@"\\fap2\mapping\GIS\Data\Third Party\AI_PRIME\Carbon_County\19N_91W\Calculated\Block_49_ToFeet_Clip.tif", Access.GA_ReadOnly);
        Band band = dataset.GetRasterBand(1);
        int width = band.XSize;
        int height = band.YSize;
        int size = width * height;
        double min = 0.00;
        double max = 0.00;
        double mean = 0.00;
        double stddev = 0.00;
        double noData = 0.00;
        int hasNodata = 0;

        band.GetNoDataValue(out noData, out hasNodata);

        var stats = band.GetStatistics(1, 0, out min, out max, out mean, out stddev);

        //Console.WriteLine($"Statistics retrieved and returned a result of {stats}");
        Console.WriteLine($"X : {width} Y : {height} SIZE: {size}");
        Console.WriteLine($"MIN : {min} MAX : {max} MEAN : {mean} STDDEV : {stddev}");
        Console.WriteLine($"HAS NO DATA : {hasNodata} Value : {noData}");
        DataType type = band.DataType;
        Console.WriteLine($"Data Type : {type}");


        float gtMean = 0; //cut
        float ltMean = 0; //fill

        float[] data = new float[size];

        var dataArr = band.ReadRaster(0, 0, width, height, data, width, height, 0, 0);

        float fMean = (float)mean;
        float fNoData = (float)noData;
        int i, j;
        for (i = 0; i < width; i++)
        {
            for (j = 0; j < height; j++)
            {
                float value = data[i + j * width];
                if (value != fNoData)
                {
                    if (value > fMean)
                    {
                        gtMean += value;
                    }

                    else if (value < fMean)
                    {
                        ltMean += value;
                    }
                }

            }
        }

        Console.WriteLine($"Sum of values above the mean {gtMean}");
        Console.WriteLine($"Sum of values below the mean {ltMean}");
    }
}

enter image description here

  • 1
    Are you sure it's every value? the number -3.402823E+38 is the min Float value, looks like NoData, it wouldn't be uncommon to have a couple of rows of NoData before you hit some data. I can't see anything intrinsically wrong with your code, perhaps run in debug with a break in the less/greater than if (note I'd make the 2nd an else if (value < (float)mean) but cast the mean as float_mean = (float)mean; and use that rather than casting on each iteration - space is cheap, time is expensive) also be aware than with very large rasters it will fall over when the raster is more than available RAM. – Michael Stimson Dec 3 '18 at 23:11
  • That makes total sense! Ive added a section to highlight the changes. In you opinion are the nested loops the appropriate way to parse through the data? I tested with a single loop on the whole data buffer and got a significantly different answer. I do not quite understand why. – LCaraway Dec 4 '18 at 2:32
  • I usually go line by line, reading into a buffer that is int width = band.XSize; by 1 and loop for each row.. it doesn't affect the sum but you also need to keep track of the number of cells for each category (greater than mean, less than mean) so divide the sum by (float)gtCount for example to get the mean of each class. – Michael Stimson Dec 4 '18 at 2:58
  • Any chance you have an example? I cant help but think more numbers are being added together than there should. The results just do not feel right if that makes since. My resolution on the raster is 10cm thus using 0.10763911106 as the pixel area (converted to sq ft) I attempt to get a rough volume of above and below the mean. These volumes are in the billions which i am having trouble accepting considering the area is only 49000 sq ft and the diff between min and max is only 23ft or so. May be a separate issue entirely. – LCaraway Dec 4 '18 at 3:29
  • 1
    You're getting the overall level, to get a cut/fill it should be gtMean += ( value - fMean); and ltMean += ( fMean - value ); - the values should be difference from the mean. – Michael Stimson Dec 4 '18 at 3:59
2

As Michael Stimson pointed out in his comments, -3.402823E+38 is the minimum value for a float and should be considered as the no data value in this instance. The inclusion of these values drastically effects the values in which we are attempting to solve for.

Below is a working example. For context, the program bellow attempts a simple cut/fill analysis on a DEM initiated at the mean elevation. In other words it is analyzing values above and below the mean. The tif being used was clipped to a specific area and the z values where converted to feet in preparation for this analysis. The raster's resolution is 10 cm. Furthermore, the code posted in the initial question was not properly summing the values above and below the mean. This has also been corrected below.

class Program
{
    static void Main(string[] args)
    {
        Gdal.AllRegister();
        Ogr.RegisterAll();

        Dataset dataset = Gdal.Open(@"path\to\raster.tif", Access.GA_ReadOnly);
        Band band = dataset.GetRasterBand(1);
        int width = band.XSize;
        int height = band.YSize;
        int size = width * height;
        double min = 0.00;
        double max = 0.00;
        double mean = 0.00;
        double stddev = 0.00;
        double noData = 0.00;
        int hasNodata = 0;

        band.GetNoDataValue(out noData, out hasNodata);


        var stats = band.GetStatistics(1, 0, out min, out max, out mean, out stddev);

        //Console.WriteLine($"Statistics retrieved and returned a result of {stats}");
        Console.WriteLine($"X : {width} Y : {height} SIZE: {size}");
        Console.WriteLine($"MIN : {min} MAX : {max} MEAN : {mean} STDDEV : {stddev}");
        Console.WriteLine($"HAS NO DATA : {hasNodata} Value : {noData}");
        DataType type = band.DataType;
        Console.WriteLine($"Data Type : {type}");


        double gtMean = 0; //cut
        double ltMean = 0; //fill
        int noDataCount = 0;

        float[] data = new float[size];

        var dataArr = band.ReadRaster(0, 0, width, height, data, width, height, 0, 0);

        float fMean = (float)mean;
        float fNoData = (float)noData;
        int i, j;
        for (i = 0; i < width; i++)
        {
            for (j = 0; j < height; j++)
            {
                float value = data[i + j * width];
                if (value != fNoData)
                {
                    if (value > fMean)
                    {
                        gtMean += value - fMean;
                    }

                    else if (value < fMean)
                    {
                        ltMean +=  fMean - value;
                    }
                }
                else
                {
                    noDataCount++;
                }

            }
        }

        Console.WriteLine($"Sum of values above the mean {gtMean}");
        Console.WriteLine($"Sum of values below the mean {ltMean}");
        Console.WriteLine($"Cells W/ Data : {size - noDataCount} Cells W/ No Data : {noDataCount}");

        double pixelArea = 0.10763911106; //cubic feet raster is (10cm x 10cm pixel)
        double dLtMean = (double)ltMean;
        double dGtMean = (double)gtMean;

        double ltMeanFt3 = dLtMean * pixelArea; //fill in cubic feet
        double gtMeanFt3 = dGtMean * pixelArea; //cut in cubic feet

        Console.WriteLine($"Cut : {gtMeanFt3} cubic feet");
        Console.WriteLine($"Fill : {ltMeanFt3} cubic feet");

        var meanDelta = gtMeanFt3 - ltMeanFt3;

        if(meanDelta > 0)
        {
            Console.WriteLine($"Excess Cut: {meanDelta} cubic feet");
        }
        else if(meanDelta < 0)
        {
            Console.WriteLine($"Insufficient Fill: {meanDelta} cubic feet");
        }
        else
        {
            Console.WriteLine("cut & fill are balanced");
        }

        Console.WriteLine("Press Any Key To Exit....");
        Console.ReadLine();
    }
} 

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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