AnsweredAssumed Answered

Cost raster creation ignores NoData areas???!!!

Question asked by geonetadmin on Dec 11, 2013
Original User: MikeEber

We are running 10.1 SP1 and trying to run routines to calculate a distance between two points.  This should be so easy to do yet even with ESRI help we have not gotten a functional model.  So here is where we stand:

Everything works as far as returning a line and getting a value from it.  However when we pass in our raster (hyrdrography of the east coast) then publish the resulting cost raster onto that map, we have large chunks of land completely ignored.  We selected a point on the hyrdrography layer and used Identify and it indicates that the value at that location is INDEED NoData.

The raster showing that the island has a value of NoData:
[ATTACH=CONFIG]29776[/ATTACH]

The resulting path between points using the cost raster:
[ATTACH=CONFIG]29777[/ATTACH]

The displayed cost layer and it's value:
[ATTACH=CONFIG]29778[/ATTACH]

So here is our code:

                // First get the origin point for computing distances
                IFeature originFeature = stationClass.Search(null, false).NextFeature();
                IPoint origin = (IPoint)originFeature.ShapeCopy;

                messages.AddMessage("Station feature defined.");
                System.Diagnostics.Debugger.Launch();

                // set spatial reference and define limit to raster area around the origin point
                Type factoryType = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
                ISpatialReferenceFactory3 spatialReferenceFactory = (ISpatialReferenceFactory3)Activator.CreateInstance(factoryType);
                ISpatialReferenceInfo spatialReference;
                int bytesRead;
                spatialReferenceFactory.CreateESRISpatialReferenceInfo(prjData, out spatialReference, out bytesRead);
                ITopologicalOperator pTopoOp = (ITopologicalOperator)origin;
                IEnvelope pExtent = pTopoOp.Buffer(buffer).Envelope;
                IRasterAnalysisEnvironment rasterAnalysis = new RasterDistanceOpClass();
                rasterAnalysis.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, pExtent);
                rasterAnalysis.SetAsNewDefaultEnvironment();

                // Empty the output class since ESRI now lets previous run hang around!
                if (outputClass.FeatureCount(null) > 0)
                {
                    IFeatureCursor cursor = outputClass.Search(null, false);
                    IFeature nextLine = cursor.NextFeature();
                    while (nextLine != null)
                    {
                        nextLine.Delete();
                        nextLine = cursor.NextFeature();
                    }
                    Marshal.FinalReleaseComObject(cursor);
                }

                // Get the gdb file then get the raster set
                //factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory");
                //IWorkspaceFactory pWorkspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
                //Gis.Libraries.DllConfigurationManager.Properties.Settings settings =
                //new Gis.Libraries.DllConfigurationManager.Properties.Settings();
                //string fileName = settings.ConfigExists() ?
                //  settings["RMW_Geodatabase"] : Properties.Settings.Default.RMW_Geodatabase;
                //IRasterWorkspaceEx workspaceExecutable = (IRasterWorkspaceEx)pWorkspaceFactory.OpenFromFile(fileName, 0);
                //IRasterDataset rasterDataset = workspaceExecutable.OpenRasterDataset("HydrographyRasterDS");
                IRasterDataset rasterDataset = GetActiveHydrography(origin);
                IRasterLayer holdLayer = new RasterLayerClass();
                holdLayer.CreateFromDataset(rasterDataset);

                bool rasterObjectId = GPUtilities.CheckCoincidence(originFeature.ShapeCopy, holdLayer, messages);
                if (!rasterObjectId)
                {
                    messages.AddMessage("Station is not coincident with the waterways and will be relocated.");
                    if (!RelocateToWaterway(ref originFeature, rasterDataset, messages))
                    {
                        messages.AddMessage("Station could not be relocated to water within 4 nautical miles.");
                        return;
                    }

                }

                messages.AddMessage("Station is coincident with the hydrography map and IRasterLayer is initialized.");

                IDistanceOp2 distanceComputer = new RasterDistanceOpClass();
                messages.AddMessage(string.Format("Beginning process of {0} target(s)", targets.Count));

                for (int index = 0; index < targets.Count; index++)
                {
                    IFeatureClass targetClass = (IFeatureClass) targets[index].Class;
                    IPoint targetPoint = (IPoint) targets[index].ShapeCopy;
                    if (!GPUtilities.CheckCoincidence(targetPoint, holdLayer, messages))
                    {
                        messages.AddMessage("Target is not coincident with the waterways and will be relocated.");
                        IFeature targetFeature = targets[index];
                        if (!RelocateToWaterway(ref targetFeature, rasterDataset, messages))
                        {
                            CreatePlaceholder(outputCursor, outputClass, messages);
                            continue;
                        }
                        targets[index] = targetFeature;
                        targetClass = (IFeatureClass)targetFeature.Class;
                        targetPoint = (IPoint) targetClass;
                    }
                    messages.AddMessage("Target is coincident with the hydrography map.");

                    IGeoDataset[] costSets = ComputeFullRaster(
                        targetClass, (RasterDataset) rasterDataset, messages);

                    // **  Temporary code for adding cost array to the map TOC
                    IApplication baseApp = (IApplication)Activator.CreateInstance(Type.GetTypeFromProgID("esriFramework.AppRef"));
                    IMxDocument doc = (IMxDocument) (baseApp.Document);

                    IMap activeMap = doc.ActiveView.FocusMap;
                    RasterLayer layer = new RasterLayerClass();
                    IRaster ds = (IRaster) costSets[0];
                    layer.CreateFromRaster(ds);
                    layer.Name = "Cost Array";
                    layer.Visible = true;
                    activeMap.AddLayer(layer);


                    IPointCollection points = new MultipointClass();
                    points.AddPoint(origin);
                    points.AddPoint(targetPoint);

                    messages.AddMessage(string.Format("Computing distance from origin to target {0}", index + 1));
                    IGeometryCollection result = distanceComputer.CostPathAsPolyline(points, costSets[0], costSets[1]);             
                    if (result == null)
                    {
                        messages.AddMessage("No result was returned for this target.");
                        CreatePlaceholder(outputCursor, outputClass, messages);
                        continue;
                    }
                    for (int subindex = 0; subindex < result.GeometryCount; subindex++)
                    {
                        Polyline data = (Polyline) result.Geometry[subindex];
                        IPolyline idata = (IPolyline) data;
                        
                        messages.AddMessage(string.Format("Their are {0} points on geometry {1} with length of {2}", data.PointCount, subindex, idata.Length));
                        messages.AddMessage("Saving geometry as new feature.");
                        IFeatureBuffer featureBuffer = outputClass.CreateFeatureBuffer();
                        IFeatureCursor cursor = outputClass.Insert(true);
                        featureBuffer.Shape = result.Geometry[subindex];
                        cursor.InsertFeature(featureBuffer);
                        cursor.Flush();
                        Marshal.FinalReleaseComObject(cursor);
                    }
                    
                }

                Marshal.ReleaseComObject(outputCursor);
            }
            catch (Exception ex)
            {
                messages.AddMessage(string.Format("Process failed: {0}", ex.Message));
                Marshal.ReleaseComObject(outputCursor);
                throw;
            }
        }


        private static IGeoDataset[] ComputeFullRaster(IFeatureClass origin, RasterDataset costArray, IGPMessages messages)
        {
            List<IGeoDataset> results = new List<IGeoDataset>();
            IDistanceOp distanceOp = new RasterDistanceOpClass();

            // Create the cost distance
            messages.AddMessage("Creating the distance cost raster...");
            object objectMissing = Type.Missing;

            results.Add( distanceOp.CostDistanceFull(
                (IGeoDataset)origin, (IGeoDataset)costArray, 
                true, false,
                false, ref objectMissing, ref objectMissing) );
            messages.AddMessage("Creating the backlink cost raster...");
            results.Add(distanceOp.CostDistanceFull(
                (IGeoDataset)origin, (IGeoDataset)costArray,
                false, true,
                false, ref objectMissing, ref objectMissing));
            return results.ToArray();
        }

Attachments

Outcomes