Monday, 2 October 2017

raster - Finding point with highest elevation in area using ArcObjects?


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.



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
};
}

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...