[GIS] Finding point with highest elevation in area using ArcObjects

arcobjectsdemelevationnetraster

Is there a way to find a point with highest elevation in specified area by using ArcObjects? Let's say that I have a raster with DEM and I would like to establish some areas of interest in this raster. In each AOI, I would need to find a point where terrain has the highest elevation and return coordinates of this point.

I looked into ArcObjects API and I found only a manual way to discover such point, that is:

  1. Establish AOIs
  2. For each AOI, find corresponding pixels in the raster
  3. Having sets of pixels defined for AOIs, iterate over each set and return a point with highest elevation

This process is doable but absolutely impractical due to performance concerns.

Best Answer

Thanks to Vince's comments I was able to perform succesfully the whole process.

As I wrote, first I wanted to establish a number of AOIs on the raster. Nothing fancy here, I got raster properties by using this small function below:

public static IRasterProps GetRasterProperties(IRasterDataset rasterDataset, int rasterBandIndex)
{
     IRasterBandCollection rasterBands = (IRasterBandCollection)rasterDataset;
     var rasterBand = rasterBands.Item(rasterBandIndex);
     return (IRasterProps)rasterBand;
}

Given the raster properties, I created a number of areas each defined by an IEnvelope. In order to split the input raster to AOIs, I used the ExtractByRectangle Geoprocessing tool. 

ExtractByRectangle extract = new ExtractByRectangle(inputRaster, envelope, path);
IGeoProcessorResult2 result = gp.Execute(extract, null) as IGeoProcessorResult2;

parameters:
inputRaster - IRasterDataset (input raster)
envelope - IEnvelope (definition of our AOI)
path - file path to resulting raster

I won't describe all Geoprocessing shenanigans, but that is where ArcObjects get very cryptic in my opinion. A number of error codes that have no obvious explanation doesn't help too. I found great code by Kirk Kuykendall that helps a lot with debugging here: Avoiding fails from ArcObjects geoprocessing with .NET?

Now, we have our small raster AOIs. I need to get the point of maximal elevation for each of them. Hence, I start with computing statistics for every AOI using the function below:

public void ComputeRasterStatistics(IRasterDataset rasterDataset, int rasterBandIndex)
{
     IRasterBandCollection rasterBands = (IRasterBandCollection)rasterDataset;
     var rasterBand = rasterBands.Item(rasterBandIndex);
     RasterStatistics = rasterBand.Statistics;
}

Statistics deliver information concerning values of raster extremes, inclusing the maximum point value of a raster. Next I convert the raster to points using the function below.

public static IFeatureClass RasterToPoints(IRasterDataset raster)
{
     IConversionOp conversionOp = new RasterConversionOpClass();
     IWorkspace shapeWS = FeatureWorkspaceHelper.CreateInMemoryWorkspace();
     var featClass = conversionOp.RasterDataToPointFeatureData((IGeoDataset)raster, shapeWS, Guid.NewGuid().ToString());
     return (IFeatureClass)featClass;
}

FeatureWorkspaceHelper.CreateInMemoryWorkspace(); is a helper function of mine, that creates an empty in-memory workspace. Have in mind this can be obviously done with 'classic' file workspace too.

Now I just need to find a point that has the maximal elevation and return it. The function below (some hardcodes there!) does that.

private IPoint GetRasterMaxElevationPoint(IFeatureClass featureClass, double val, int elevationIndex)
{
    IQueryFilter queryFilter=new QueryFilterClass();
    queryFilter.WhereClause = "GRID_CODE >= " + (val - 0.01).ToString();

    var cursor = featureClass.Search(queryFilter, true);
    IFeature feature = null;
    IGeometry shape = null;
    double maxValue = double.MinValue;
    while ((feature = cursor.NextFeature()) != null)
    {
        if ((double) feature.Value[elevationIndex] > maxValue)
        {
            shape = feature.Shape;
            maxValue = (double) feature.Value[elevationIndex];
        }
    }

    return new PointClass()
    {
        X = shape.Envelope.LowerLeft.X,
        Y = shape.Envelope.LowerLeft.Y,
        SpatialReference = shape.SpatialReference
    };
}
Related Question