# Finding maximum value on raster transect between two points using PyQGIS

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.

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')
``````

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:

• 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...
– S.E.
Apr 14, 2022 at 19:35
• 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.
– Matt
Apr 14, 2022 at 20:34
• I just meant I'm not familiar with all these functions and methods... with all the added details it's more digest now !
– S.E.
Apr 14, 2022 at 23:04
• I see :). Let me know if you need anything made clearer
– Matt
Apr 14, 2022 at 23:08