Here is an example using the QgsZonalStatistics class to find the highest elevation from a DEM within a given distance from the clicked coordinate. Specifically, we can use the calculateStatistics() static method which takes a QgsGeometry argument rather than a vector layer.
It's also worth bearing in mind that the clicked point will be in the canvas CRS. I assume that you are getting your user-defined radius distance in meters, so this example includes a check to see if the canvas CRS map units are not meters and uses the QgsDistanceArea class to convert to the appropriate unit if necessary. This answer also shows an example of using the QgsCoordinateTransform class to transform the buffer geometry to match that of the raster layer if those are different.
To also extract the locations of the pixels with the maximum value (which of course there may be more than one), we can read the raster into a numpy array and use the numpy.where() method to get the indices of all pixels with the maximum value, then calculate the coordinates of the pixel centroids based on the raster origin & pixel dimensions (georeferenced rasters do not inherently store the location of every pixel). To make this process a bit more efficient, this example makes a processing.run() call to "gdal:cliprasterbyextent" to clip the raster to the bounding box of the buffered click point, before running the zonal stats calculation.
from osgeo import gdal
import numpy as np
class CustomMapTool(QgsMapToolEmitPoint):
def __init__(self, iface):
self.iface = iface
self.canvas = self.iface.mapCanvas()
super(CustomMapTool, self).__init__(self.canvas)
self.buffer_distance = 5000.0# Hard coded for this example (units are meters)
self.da = QgsDistanceArea()
def canvasReleaseEvent(self, e):
clicked_point = e.mapPoint()
lyr = self.iface.activeLayer()
if lyr.type() != QgsMapLayerType.RasterLayer:
return
max_val, coords = self.get_zonal_max(lyr, clicked_point)
print(max_val, coords)
def get_zonal_max(self, rl, click_pt):
canvas_crs = self.canvas.mapSettings().destinationCrs()
if canvas_crs.mapUnits() != QgsUnitTypes.DistanceMeters:
self.da.setSourceCrs(canvas_crs, self.canvas.mapSettings().transformContext())
self.da.setEllipsoid(canvas_crs.ellipsoidAcronym())
radius = self.da.convertLengthMeasurement(self.buffer_distance, canvas_crs.mapUnits())
else:
radius = self.buffer_distance
pt_geom = QgsGeometry.fromPointXY(click_pt)
buffer = pt_geom.buffer(radius, 50)
if rl.crs() != canvas_crs:
buffer = self.transformed_geom(buffer, canvas_crs, rl.crs())
bb = buffer.boundingBox()
ext = f'{bb.xMinimum()},{bb.xMaximum()},{bb.yMinimum()},{bb.yMaximum()}'
params = {'INPUT':rl,
'PROJWIN':ext,
'OVERCRS':False,
'NODATA':-999,
'OPTIONS':'',
'DATA_TYPE':0,
'EXTRA':'',
'OUTPUT':'TEMPORARY_OUTPUT'}
temp_rl = QgsRasterLayer(processing.run("gdal:cliprasterbyextent", params)['OUTPUT'], '', 'gdal')
pixel_size_X = temp_rl.rasterUnitsPerPixelX()
pixel_size_Y = temp_rl.rasterUnitsPerPixelY()
res = QgsZonalStatistics.calculateStatistics(temp_rl.dataProvider(),
buffer,
pixel_size_X,
pixel_size_Y,
1,
QgsZonalStatistics.Statistic.Max)
if not list(res):
return NULL, []
max_val = res[list(res)[0]]
max_coords = []
highest_pts = self.get_pixel_coords_from_val(temp_rl, max_val)
for pt in highest_pts:
pt_geom = QgsGeometry.fromPointXY(pt)
if pt_geom.within(buffer):
max_coords.append((pt.x(), pt.y()))
return max_val, max_coords
def transformed_geom(self, g, src_crs, dest_crs):
xform = QgsCoordinateTransform(src_crs, dest_crs, QgsProject.instance())
g.transform(xform)
return g
def get_pixel_coords_from_val(self, r_lyr, val):
path = r_lyr.source()
ds = gdal.Open(path)
geotransform = ds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
arr = ds.ReadAsArray()
value_indices = list(zip(*np.where(arr == val)))
all_points = []
for pp in value_indices:
x = (originX+(pp[1]*pixelWidth))+pixelWidth/2
y = (originY+(pp[0]*pixelHeight))+pixelHeight/2# Pixel height is negative
all_points.append(QgsPointXY(x, y))
ds = None
return all_points
map_tool = CustomMapTool(iface)
iface.mapCanvas().setMapTool(map_tool)
This script will print the maximum value plus a list containing tuples of coordinate pairs for each of the pixels having the maximum elevation within the buffer distance of the click point, which you can then do with whatever you like.
Just as an example, in the screencast below, I added some logic to the script to create features and add them to a temporary layer.
QgsGeometryobject from the clicked point and buffer by the radius distance instead of creating aQgsCircle. Don't forget you may need to convert the buffer distance to canvas CRS units and then transform to the raster CRS if different. Let me know if this would suit your needs and I will add a full answer tomorrow.