0

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.

1

1 Answer 1

6

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:

4
  • 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.