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.

• Taras: Please, do not forget about "What should I do when someone answers my question?"

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

enter image description here

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

• S.E.: Thanks a lot, that's a concise, sophisticated way... I just can't manipulate these Qgs tuples easily, if I apply a CRS transformation to the initial points before your 1st line, should the result be the same as if I applied it just before ` raster.dataProvider().sample(vert, 1)[0]` I'm asking because my results are 1 meter off compared to this method...
• Matt: I'm not sure what you mean by can't manipulate these Qgs tuples easily. Do you need the output in a different format? It sound slike both those transformations should be the same, indeed. I added a coordinate transformation to my answer.
• S.E.: I just meant I'm not familiar with all these functions and methods... with all the added details it's more digest now !
• Matt: I see :). Let me know if you need anything made clearer