Geospatial Raster reclassification

Often in image processing, one wishes to assign a single unique number to a particular range of raster values, for example when creating different types of masks. In remote sensing, this method of converting pixel values into classes is called threshold classification. An example of this is a water mask from satellite imagery, where pixels representing reservoirs take values of 1. The others have values of 0.

How do you do this in R? Nothing simpler than that. First we load an image using the raster function in the raster library. We will use a Sentinel satellite image taken in band 8 (infrared). We will display the loaded image:

plot(b8,col=grey.colors(255))

In a GIS software we have previously vectorized four areas where water bodies occur. We will use this layer to determine the threshold for the mask. We load the created layer with the function readOGR. With the function extract we read the pixel values for the areas with waters and write them into the variable ext:

ext = extract(b8, obszary)

The ext variable contains a list of four vectors with the read pixel values for each area:

ext

[[1]]

  [1] 274 274 276 273 266 267 269 266 275 260 273 270 266 267 266 270 272

 [18] 258 262 266 259 259 253 260 274 267 267 267 268 265 266 264 258 265

…

[[2]]

  [1] 264 262 261 256 258 264 261 256 256 260 262 255 258 256 256 256 256

 [18] 258 258 259 261 257 257 256 257 256 256 256 256 256 256 257 259 256

…

[[3]]

…

[[4]]

…

The unlist function converts the list into a single vector containing the pixel values for all areas:

ext = unlist(ext)

ext

[1] 290 296 290 288 290 286 281 276 274 276 283 285 282 281 274 268

  [17] 260 268 276 275 282 281 277 270 269 264 263 264 273 270 270 275

…

Let’s calculate the average of the values in ext:

avg = mean(ext)

The average value will be our threshold. Let’s reclassify the image pixels assuming that all values below the average take the value 1. The others will take the value 0. We will use the reclassify function. For this function, we need to define the raster that we want to reclassify and the vector of change in value, which is of the form:

c(left_threshold, right_threshold, new_value, left_threshold, right_threshold, new_value, …)

Thresholds can be infinite values defined as -Inf and Inf. Finally, we create our mask using the expression:

mask = reclassify(b8,c(-Inf,avg,1,avg,Inf,0))

Show result:

plot(mask)

We have obtained a mask of water bodies based on near infrared satellite imagery. We can use this tool to create arbitrary masks for different types of raster data.