Identify the coordinates of the max values resulting from using the Zonal Statistics tool

Geographic Information Systems Asked by JAT86 on December 13, 2020

I used the Zonal Statistics tool to extract the max values of a raster/DEM layer. The attribute table shows the result, but the problem now is to identify the coordinates of the said maximum elevation. I tried looking for solutions online and saw Find the coordinates of the max n values of a raster layer but it does not apply to QGIS.

How do I identify these max values?

2 Answers

So, as output of the zonal stat, you have the max value of the raster in the attribute table for each polygon.

Convert the polygons to raster (raster > conversion > rasterize) with the same pixel size and extend than you raster dataset (in the parameter "output extent (xmim,ymin,xmax,ymax)" press the button on the right and "use layer/canvas extent") and based on the value of your "_max" field in the parameter "field for use for a burn in value [optional]" .

Use raster calculator to get all values where your input_raster = converted_zonal_stat_raster

Convert this result to points, either with "Processing > SAGA > Raster values to points" (which gives you the XY coordinate) or with "vector creation > raster pixels to points" (faster, but then you need to compute the coordinate values in the attribute table)

Answered by radouxju on December 13, 2020

If you know a little bit of using python in Qgis, you can easily use this code and get the coordinates. If you don't know how to use the Python console, follow the steps below:

  1. Press CTRL+ALT+P on keyboard. This will open the console.
  2. Click 'Show Editor' in console panel. Look picture.
  3. Paste the following code to the editor. And replace 'your_raster_name' variable with your own raster name.(Look to the code)
  4. And run it.(Icon: Green Triangle)
  5. After running, Coordinates and elevation will be displayed as text in the console
  6. Note: Raster Crs and Project Crs must be same.

raster = QgsProject.instance().mapLayersByName("your_raster_name")[0]
xmin = raster.extent().xMinimum()
ymin = raster.extent().yMinimum()
xunit = raster.rasterUnitsPerPixelX()
yunit = raster.rasterUnitsPerPixelY()
vals = list()
for i in range(0,int(raster.width())):
    for j in range(0,int(raster.height())):
        val, res = raster.dataProvider().sample(QgsPointXY(xmin+i*xunit+xunit/2,ymin+j*yunit+yunit/2), 1)
result = max(vals, key=lambda x: x[2])
print(f'Max Elevation is {result[2]}')
print(f'Coordinates are {result[0]} by {result[1]}')

enter image description here

Answered by Xi Jin on December 13, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP