Selecting features from a .shp file based on rasters pixel values

Question

I am on windows 10, using QGIS 3.18.

I have a shapefile with multiple entries (polygons) and a raster with a few pixel values. I would like to use a "select by ... " to select the entries within the '.shp' that intersect with specific pixel values.

For example, try answering the question: Which entries (lines in my attribute table) intersect with "1"s in raster 'raster_name' (raster with 0 and 1s). In example 2, how to do the same but with rasters with a few different values (let's say 1 to 10) when I would like to 'select by ...' intersect with pixels with values equals '9'.

I got as far as:

'select by expression'

intersects($geometry, layer_property('raster_name', 'extent'))

With that, I could select every entry that intersects with a raster extent.

I also tried:

intersects($geometry, raster_name@band = value)

And it did not cam as a valid expression.

Answer 1

Here is an example of the workflow.

1. Use gdal_translate to convert the raster to XYZ format.
2. Filter out nodata values from the resulting file.
3. Add a header to each set ("X Y value").
4. Load the XYZ as a CSV file of Points. This converts the raster to a vector.
5. Do a spatial join using Join Attributes by Location against the Polygon layer.

You could read all valid data in one layer and set a query for each data value of you prefer.

Filtering

gawk '{if($3>-1)print($0)}' test.xyz >data.csv

may help if you have awk or gawk, and a negative nodata value.

Example

Here are some polygons and a raster. Only some points in the raster have values. The rest are -999 nodata values.

After converting to XYZ, filtering out the nodata values and reloading as CSV Points, the loaded points overlay the raster correctly.

Running the tool Join Attributes by Location:

Viewing the resulting table, ID is the polygon ID, and field3 is the point value.

Answer 2

Simply use this expression to get the raster value at a certain point (and convert lines etc. to points, see below). Be sure your vector- and raster-layers have the same CRS (projection):

raster_value('layer', 1,$geometry)

• 'layer' is the name of the raster-layer
• 1 is the number of the raster-band
• $geometry is the point where the raster value should be returned. You can define any other point-geometry: based on QGIS expression, you could e.g. use the start-point of a line: start_point ($geometry)

This last point means that you should convert your line to points and sample the raster value at each point. Thus densify your line (create points in a regular interval) in a distance slightly shorter than the size of your pixels (in case the line just passes through a corner of the pixel). Do this using Menu Processing / Toolbox / Points along geometry - see screenshot 2.

Than query the pixel value for each point to see if any fulfills the condition.

Screenshot: the raster-value at the pixel where the red points are located is calculated, see preview at the bottom (red rectangle):

Screenshot 2: line (red), interpolated Points along the line (blue) and selected points where raster-value is > 3000 (yellow points):