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.

Taras
  • 32,823
  • 4
  • 66
  • 137
S.E.
  • 336
  • 1
  • 8

1 Answers1

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:

Taras
  • 32,823
  • 4
  • 66
  • 137
Matt
  • 16,843
  • 3
  • 21
  • 52
  • 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 '22 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 '22 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 '22 at 23:04
  • I see :). Let me know if you need anything made clearer – Matt Apr 14 '22 at 23:08