AdrianK

Problem generating IRasterBand statistics

Discussion created by AdrianK on Jun 29, 2011
Latest reply on Jan 11, 2012 by AdrianK
I have an issue with the generated statistics when creating a new raster. I find that I get a different result when I generate statistics through the IRasterBand as opposed to through ArcMap Toolkit. The ArcMap Toolkit generates the correct stats which displays the raster correctly. However, when I generate the stats through IRasterBand the raster doesn't display correctly until I regenerate them through ArcMap Toolkit. Very annoying after a batch process. I've included the code below in case it sheds some light on the issue. The code is converting an ENVI file to grid. The envi class, which isn't included, is simply a BinaryStream obtaining the values for the GRID pixels.

Another possibly related issue occurs after producing a large GRID. A 'Function evaluation is disabled because a previous function evaluation timed out' error causes the program to hang when trying to display a 'Finished' message box. The raster is actually produced fine (except for previous stats issue). This apparently has something to do with the debugger but I've yet to resolve it.

If anyone has any ideas about these two issues I'd like to hear it.

Cheers

Adrian

        private bool Output_ENVI2GRID(string directoryName, string fileName, ref ENVI_Reader envi)
        {
            IPixelBlock3 pixelBlock3 = null;
            IPnt pixelBlockOrigin = new PntClass();
            IPnt pixelBlockSize = new PntClass();
            IRasterBandCollection rasterBandCollection = null;
            IRasterProps rasterProps = null;
            IRasterDataset rasterDataset = null;
            IPoint originPoint = new PointClass();
            IRasterWorkspace2 rasterWorkspace2;
            ISpatialReference spatialRef;
            int blockCount;
            int lastBlockLineCount = 0;
            int lineCount = 128; ;
            IRasterEdit rasterEdit = null;
            IRaster raster = null;
            IRasterBand rasterBand;
            string outPath = @"C:\TEMP\"
            System.Array pixelData;
            double pixVal = 0;
            double min = 999999999999;
            double max = -999999999999;
            double sum = 0;
            double count = 0;
            double sumSq = 0;


            try
            {
                if (envi.mapInfo_Projection_Info == null || envi.mapInfo_Projection_Info == "")
                {
                    spatialRef = new UnknownCoordinateSystemClass();
                }
                else
                {
                    spatialRef = GetSpatialReference(GetGeoProjectionType(envi.mapInfo_Projection_Info));
                }

                originPoint.PutCoords(envi.mapInfo_MapX, envi.mapInfo_MapY - ((double)envi.header_Lines * envi.mapInfo_PixelX));

                // Create the dataset
                rasterWorkspace2 = null;
                rasterWorkspace2 = SetRasterWorkspace(outPath);

                //create empty raster dataset
                rasterDataset = rasterWorkspace2.CreateRasterDataset(fileName, "GRID", originPoint, (int)envi.header_Samples, (int)envi.header_Lines, envi.mapInfo_PixelX, envi.mapInfo_PixelY, 1, GetGridPixelType(envi.header_DataType), spatialRef, true);

                rasterBandCollection = (IRasterBandCollection)rasterDataset;
                rasterBand = rasterBandCollection.Item(0);
                rasterProps = (IRasterProps)rasterBand;
                envi.OutputNoDataValue = rasterProps.NoDataValue;
                raster = rasterDataset.CreateDefaultRaster();
                
                rasterEdit = (IRasterEdit)raster;
                
                //work out the size of the envi matrix and determine if needs 128 line pixelblocks
                blockCount = (int)Math.Ceiling((double)(envi.header_Lines / (double)lineCount));

                if (blockCount == 1)
                {
                    pixelBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
                    lineCount = rasterProps.Height;
                }
                else
                {
                    pixelBlockSize.SetCoords(rasterProps.Width, lineCount);
                    lastBlockLineCount = (int)(((double)(envi.header_Lines / (double)lineCount) % 1) * lineCount);
                }
                
                //loop through pixel block and fill pixel values
                for (int i = 1; i <= blockCount; i++)
                {
                    
                    //reconfigure pixel block for last block dimensions
                    if (i == (blockCount) && blockCount > 1)
                    {
                        lineCount = lastBlockLineCount;
                        pixelBlockSize.SetCoords(rasterProps.Width, lastBlockLineCount);
                    }

                    pixelBlock3 = (IPixelBlock3)raster.CreatePixelBlock(pixelBlockSize);

                    pixelBlockOrigin.SetCoords(0, (i - 1) * 128);

                    // Get pixeldata array 
                    pixelData = (System.Array)pixelBlock3.get_PixelData(0);

                    // Loop through all the pixels and assign value  
                    for (int k = 0; k < lineCount; k++)
                    {
                        for (int j = 0; j < rasterProps.Width; j++)
                        {
                            pixelData.SetValue(envi.Read_Pixel(), j, k);
                            
                            //gather info for stats
                            if (!pixelData.GetValue(j, k).Equals(envi.OutputNoDataValue))
                            {
                                pixVal = Convert.ToDouble(pixelData.GetValue(j, k));
                                count++;
                                if (pixVal < min) { min = pixVal; }
                                if (pixVal > max) { max = (pixVal); }
                                sum += pixVal;
                                sumSq += Math.Pow(pixVal,2);
                            }
                        }
                    }

                    pixelBlock3.set_PixelData(0, (System.Array)pixelData);

                    //Write the pixel block.
                    rasterEdit.Write(pixelBlockOrigin, (IPixelBlock)pixelBlock3);

                    pixelData = null;
                    pixelBlock3 = null;

                }

                GC.Collect();

                if (!rasterBand.Statistics.IsValid)
                {

                    //rasterBand.Statistics.SkipFactorX = 1;
                    //rasterBand.Statistics.SkipFactorY = 1;
                    rasterBand.Statistics.Recalculate();

                    //rasterBand.Statistics.Minimum = min;
                    //rasterBand.Statistics.Maximum = max;
                    //rasterBand.Statistics.Mean = sum / count;
                    //rasterBand.Statistics.StandardDeviation = Math.Sqrt((sumSq - Math.Pow(sum, 2) / count) / (count - 1));
                }
                return true;

            }
            catch (Exception ex)
            {
                 throw new Exception(ex.Message + "\nError in Output_ENVI2GRID");
            }
            finally
            {
                pixelBlock3 = null;
                pixelBlockOrigin = null;
                pixelBlockSize = null;
                rasterBandCollection = null;
                spatialRef = null;
                rasterProps = null;
                rasterEdit = null;
                Common.ReleaseObject(raster);
                raster = null;
                rasterDataset = null;
                originPoint = null;
                rasterWorkspace2 = null;
                blockCount = 0;
                lastBlockLineCount = 0;
                lineCount = 0;
            }
        }
        
        public IRasterWorkspace2 SetRasterWorkspace(string sName)
        {
            // Given a pathname, returns the raster workspace object for that path
            //If not a valid raster workspace returns nothing. Caller must test

            IWorkspaceFactory pWKSF = null;
            IRasterWorkspace2 pWKS = null;

            try
            {

                pWKSF = new RasterWorkspaceFactory();
                pWKS = (IRasterWorkspace2)pWKSF.OpenFromFile(sName, 0);

                return pWKS;
            }
            catch (Exception ex)
            {
                throw new Exception("Multiple Raster Bands: " + ex.Message);
                //        Set SetRasterWorkspace = Nothing
            }
            finally
            {
                pWKSF = null;
                pWKS = null;
            }

        } 
        
        public ISpatialReference GetSpatialReference(GeoProjectionType GeoCorrdSys)
        {
            ISpatialReferenceFactory3 spatialReferenceFactory = new SpatialReferenceEnvironmentClass();
            ISpatialReference spatialReference;
            try
            {
                spatialReference = spatialReferenceFactory.CreateSpatialReference((int)GeoCorrdSys);
                return spatialReference;
            }
            finally
            {
                spatialReferenceFactory = null;
                spatialReference = null;

            }
        }        

Outcomes