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:
- Establish AOIs
- For each AOI, find corresponding pixels in the raster
- 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