Sample Multiple Points from Raster Layer Python QGIS 3.26.1

Question

I have the coordinates for 30000 locations I would like to to sample from raster layers. These locations are not located next to one another. This is currently done with an external python script and takes hours to loop through the array of coordinates, creating a QgsPointXY object and then sampling the raster layer at that single location.

values = [0.0] * len(x_coords)
for i in range(len(x_coords)):
 layer = QgsRasterLayer(path)
 values[i] = layer.dataProvider().sample(QgsPointXY(x_coord[i], y_coord[i]), 1)[0]

Is there a way to sample multiple points from a raster layer in QGIS 3.26.1 with python?

If this is possible, is it faster than iterating through each point separately?

I have looked through the QGIS Python API documentation for QGIS and cannot find anything

• Erik: Do the points already exist? And what exactly do you mean by "sample points"?
• Ethan: @Erik At present I have the x and y coordinates, then create the QgsPointXY right before I sample it.
• Erik: And "sampling" is defined as "attaching raster values to points" for you?
• Ethan: @Erik I would like to obtain the value of the raster layer at that point. The QGIS Python API do this for a single location and specific band with the sample() method and use the language "sample"
• Erik: Have you had a look at the native rastersampling tool?

Answer 1

Short answer:

Yes, Yes (way faster), https://docs.qgis.org/3.22/en/docs/user_manual/processing_algs/qgis/rasteranalysis.html#sample-raster-values

Long answer:

I would also prefer using the algorithm ('Sample raster values') from the processing toolbox, but note well that all processing algorithms can be executed from the python console. Accessing the processing toolbox gives you useful information as well:

1. Mouseover the algorithm reveals the algorithm ID, here: native:rastersampling

enter image description here

2. Execute the algorithm from the gui:

enter image description here

3. Have a look at the Log:

enter image description here

Note first the execution time, 0.70sec. for 30k points (way faster than a couple hours)

Note also the Input parameters and the Results, both python dictionaries.

With this informations gathered, the python code is almost self-explanatory:

point_lyr = QgsProject.instance().mapLayersByName('my_points')[0]
raster_lyr = QgsProject.instance().mapLayersByName('my_raster')[0]
params = {'COLUMN_PREFIX':'SAMPLE_',
'INPUT' : point_lyr,
'OUTPUT' : 'TEMPORARY_OUTPUT',
'RASTERCOPY' : raster_lyr}
result = processing.run('native:rastersampling', params)
result_lyr = result['OUTPUT']
QgsProject.instance().addMapLayer(result_lyr)

Results in

enter image description here

... and welcome to GIS.SE btw!