I have a shp with more than 1700 polygons on a elevation raster.
I need to plot a point that represents the highest altitude location within the polygons, and then export those points as a new shp.
[GIS] How to plot a point on the highest elevation of a raster on a polygon. QGIS
altitudepolygonqgisraster
Related Solutions
Is Python an option?
Use RasterIO (a Python GDAL/ numpy bridge) to load the raster to NumPy array, then use numpy.amax() to find the maximum value, followed by numpy.where() to find the row/column indices, then calculate the lat and lon from the raster extents.
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
};
}
Best Answer
I first thought of recommending an approach here Extract highest point in raster and convert to point vector, but the large number of your polygons (1700) gives another challenge, so let me suggest a three-steps way as below.
(1) Run SAGA
Raster values to pointstool (in Processing Toolbox | SAGA | Vector <-> Raster). Select your elevation raster layer (let's call itelevation) as Grids. You may of may not assign your polygon layer toPolygonsoption. This will give you an extracted point layer (Shapes).(2) Start
Join attributes by locationtool. SelectShapeslayer as the Target vector layer and polygon layer as Join vector layer. (Any predicate will do: e.g. intersects). This will return aJoined layerpoint shapefile which has polygon id (let's sayfid) field added onto theShapeslayer.(3) Start
Extract by expressiontool (in Processing Toolbox | QGIS geoalgorithms | Vector selection tools if it is QGIS 2.18). Expression would be"elevation" = maximum("elevation", group_by:= "fid"). This will return a new point shapefile layerExtracted(expression)which is the highest point within each polygon.Please note the above expression can extract multiple "highest" points if there are cells of the same highest elevation.