Finding maximum value on raster transect between two points using PyQGIS

Question

I have a raster layer and two points with their coordinates, and I want to find the maximum raster value on the line between these two points.

Is there a built-in function in PyQGIS to find the pixel with the maximum value between those points (and optionally retrieve its coordinates) ?

I'm extracting the pixel value for a point with raster.dataProvider().sample(QgsPointXY(x, y), 1), so I was thinking maybe raster.dataProvider().sample(QgsLineString(QgsPointXY(x1, y1), QgsPointXY(x2, y2)), 1) would do it, but I get an error :

TypeError: QgsRasterDataProvider.sample(): argument 1 has unexpected type 'QgsLineString'

Is there an easy way to do this?

Edit: I'm doing this by dividing the transect in x points and calculating the raster value for each, but a PyQGIS function could do it better/faster I suppose.

Answer 1

You could do something like this:

# get reference to project
p = QgsProject.instance()

# get raster layer
raster = p.mapLayersByName('dem')[0]

# using the linestring from your example
linestring = QgsLineString(QgsPointXY(x1, y1), QgsPointXY(x2, y2))

# make QgsGeometry from the linestring
linestring_geom = QgsGeometry.fromPolyline(linestring)

# create source CRS (change EPSG code as needed)
source_crs = QgsCoordinateReferenceSystem(4326)
# get target CRS directly from the raster layer 
target_crs = raster.crs()

# create a QgsCoordinateTransform instance
tr = QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())

# transform the geometry to the projected CRS by passing the `QgsCoordinateTransform` instance to the `transform` method of the geometry
linestring_geom.transform(tr)

# densify (add vertices to) the linestring at a specified distance smaller than your raster cell size (I used 1 because my data is high resolution)
dist = 1
dense = linestring_geom.densifyByDistance(dist)

# make QgsPointXYs from the dense vertices
dense_verts = [QgsPointXY(vert) for vert in list(dense.vertices())]

# sample raster for every dense point, return a tuple of the QgsPointXY and the raster value 
# points that sample NoData cells are discarded (raster.dataProvider().sample(vert, 1)[1] returns False)
smp = [(vert, raster.dataProvider().sample(vert, 1)[0]) for vert in dense_verts if raster.dataProvider().sample(vert, 1)[1]]

# sort the list of tuples by the raster value
smp_sorted = sorted(smp, key = lambda x: x[1])

# get the last element (largest raster value)
max_val = smp_sorted[-1]

# you can get the various values from the result using indexing  (\n is the newline character for the print statement)
# max_val[0] is a QgsPointXY which has .x() and .y() methods
# max_val[1] is the largest raster value
print(f'The maximum value is {max_val[1]}\nx: {max_val[0].x()},\ny:{max_val[0].y()}')

# make point layer of highest raster value and add to map
lyr = QgsVectorLayer(f'?query=SELECT {max_val[1]} AS max_val, SetSRID(ST_GeomFromText(\'{max_val[0].asWkt()}\'), {int(target_crs.authid().replace("EPSG:",""))})', 'highest raster value', 'virtual')
p.addMapLayer(lyr)

Output:

The maximum value is 1501.905029296875
x: 584795.5171731369,
y:7305533.836840233

Note: The point coordinates will not be the pixel centroid, but they will fall inside the pixel with the maximum value.

---

References:

• https://stackoverflow.com/a/36955637/4227297
• https://gis.stackexchange.com/a/314864/128665
• https://gis.stackexchange.com/a/428574/128665