Extract cells above certain values of raster file

Question

I have a very large raster layer and with the values for noise emission. I'd like to extract all cells which are above a certain threshold value (e.g. 50 db) and convert them to an ESRI Shapefile so that I can intersect then with a line shapefile at the next step. The files should be as small as possible. My TIF file is 600MB.

I'd prefer to solve the issue using R.

Answer 1

As previously commented, you may accomplish all your task using R exclusively; in the following code, a raster is created and then filtered to values above the threshold, all other values will become NA, then the pixels are masked to the lines. Raster works now in recent versions with both sf and sp objects, the code uses the latter kind.

library(raster)

r <- raster(ncol = 36, nrow = 18, vals = runif(648, min = 0, max = 70))
cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25))
cds2 <- rbind(c(80,20), c(140,60), c(160,0), c(140,-55))
lines <- spLines(cds1, cds2)


par(mfrow = c(1,3))
plot(r, main = "Original raster")
# to filter out values < 60
r[r[] < 50 ] = NA # check the use of braces to acces values of the raster
plot(r, zlim = c(0,70), main = "Filtered raster")

# to get the values of pixels that touch the lines
extract(r, lines)

# to filter the rasters to the lines
rf = mask(r, lines) 
plot(rf, main = "Extracted to lines")
plot(lines, add = T)

Answer 2

The following solution assumes that you have a GeoTiff with raster values, which you want to filter based on a threshold (e.g. 50) and convert the resulting pixels to an ESRI Shapefile for later use.

##load libraries
library(raster)
library(rgdal)

#read in raster
inputimage <- raster('input_raster.tif')

##check the image in R 
plot(inputimage) 

##mask the image (two steps)
tmpfilter <- inputimage < 50
filtered_image <- mask(inputimage, tmpfilter, maskvalue=1)

##check the filtered image 
plot(filtered_image)

##convert the output to a SpatialPointsDataFrame
p <- rasterToPoints(filtered_image)

##write the output to an ESRI Shapefile (writeOGR)
writeOGR(p, 'filtered_points.shp', driver='ESRI Shapefile', layer='filtered_points')