1. Introduction

In this notebook I explore functionalities provided by the imager package developed by Simon Barthelmé as well as the imagerExtra package developed by Shota Ochi.

I applied here these packages’ functions to improve remote sensing image visualization.

In a earlier notebook, we saw that the terra package provide a few functions to improve image visualization, that it, linear stretch and histogram equalization. Some times, these types of stretch enhancement do not meet users’ needs.

I will illustrate the topic using a spatial subset of a Landsat 8 scene collected in 2013. The subset covers the area known as Montes de Maria, in the Colombian Caribbean region.

2. Reading the image

Let’s use the raster library to create a raster object from a Landsat multispectral TIF file obtained in a previous notebook.

library(raster)
## Loading required package: sp
landsat13 <- brick("./landsat/montes-landsat.TIF")

Let’s check what is landsat13:

landsat13
## class      : RasterBrick 
## dimensions : 3392, 3409, 11563328, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 422745, 525015, 1019775, 1121535  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : /Users/ivanlizarazo/Documents/ivan/PR/datos/landsat/montes-landsat.TIF 
## names      : ultra.blue,  blue, green,   red,   NIR, SWIR1, SWIR2 
## min values :          0,     0,     0,     0,     0,     0,     0 
## max values :      20788, 21783, 23294, 27555, 43270, 65535, 65535

We can see the spatial resolution, extent, number of layers, coordinate reference system and more.

Note that digital values range from 0 to 65535. In order to obtain apparent reflectance values it is necessary to make a conversion:

Reflectance conversion factors

3. Image conversion

Now, let’s use the imager package to convert a single band of a raster object into a cimg object:

library(imager)
img5 <- round(as.cimg(landsat13[[5]])*2E-03)

What is img5?

img5
## Image. Width: 3170 pix Height: 3154 pix Depth: 1 Colour channels: 1
img6 <- round(as.cimg(landsat13[[6]])*2E-03)
img4 <- round(as.cimg(landsat13[[4]])*2E-03)

4. Image display

Let’s plot the three bands:

layout(t(1:3))
plot(img5, main= "nir band")
plot(img6, main= "swir1 band")
plot(img4, main= "red band")

Bands appear dark because we are displaying them using default parameters (i.e., there is not contrast stretch).

Let’s append the individual bands to create a three-band image:

ilandsat564 <- imappend(list(img5,img6,img4),"c")
ilandsat564
## Image. Width: 3170 pix Height: 3154 pix Depth: 1 Colour channels: 3
plot(ilandsat564, main="Landsat RGB564 -default visualization") #recombine and plot

Not too many things to see, isn’t it?

5. Image histogram

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
landsat.df <- as.data.frame(ilandsat564)
landsat.df <- mutate(landsat.df,channel=factor(cc,labels=c('NIR','SWIR1','RED')))
ggplot(landsat.df,aes(value,col=channel))+geom_histogram(bins=30)+facet_wrap(~ channel)

What we immediately see from these histograms is that the digital values are located to the left part, that is, most of the middle and high values are in a sense under-used: there’s very few pixels with middle or high values.

Contrast stretch solves the problem by making histograms wider.

6. Contrast enhancement

It is said that contrast enhancement improve the perceptibility of objects in the scene by enhancing the brightness difference between objects and their backgrounds. Read a complete explanation here.

As an illustration of what contrast enhancement means, let’s do a simple example:

library(imagerExtra)

BalanceSimplest saturates a percentage sleft % of the pixels on the left side of the histogram, and a percentage sright % of the pixels on the right side of the histogram.

See Limare et al. (2011) for detail.

layout(matrix(1:2, 1, 2))
plot(img5, main = "Original NIR band")
BalanceSimplest(img5, 2, 2) %>% plot(main = "sleft = 2.5, sright = 2.5")

layout(t(1:2))
hist(img5,main="Histogram of band5")
b5 <- BalanceSimplest(img5, 2, 2)
hist(b5,main="Histogram of balanced band5")

Let’s do this for the three-band image:

b6 <- BalanceSimplest(img6, 2, 2)
b4 <- BalanceSimplest(img4, 2, 2)
n564 <- imappend(list(b5,b6,b4),"c")
layout(t(1:2))
plot(ilandsat564, main="L8 RGB564 - No stretch") 
plot(n564, main="L8 RGB564 - Balanced")

The simple balanced false-color composite reveals much more about the landscape than the earlier no stretch composite. This explains why the contrast stretch techniques are most used to support visual interpretation of satellite images.

7. Further processing

Other functions for contrast enhancement provided by imagerExtra are:

You may be interested in applying some of these functions to your own image(s).

7. Reproducibility

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] imagerExtra_1.3.2 dplyr_1.0.4       ggplot2_3.3.3     imager_0.42.7    
## [5] magrittr_2.0.1    raster_3.4-5      sp_1.4-5         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        compiler_4.0.4    pillar_1.4.7      highr_0.8        
##  [5] tools_4.0.4       digest_0.6.27     evaluate_0.14     lifecycle_1.0.0  
##  [9] tibble_3.0.6      gtable_0.3.0      lattice_0.20-41   pkgconfig_2.0.3  
## [13] png_0.1-7         rlang_0.4.10      igraph_1.2.6      DBI_1.1.1        
## [17] yaml_2.2.1        rgdal_1.5-23      xfun_0.21         readbitmap_0.1.5 
## [21] bmp_0.3           withr_2.4.1       stringr_1.4.0     knitr_1.31       
## [25] fftwtools_0.9-11  generics_0.1.0    vctrs_0.3.6       tidyselect_1.1.0 
## [29] grid_4.0.4        glue_1.4.2        R6_2.5.0          jpeg_0.1-8.1     
## [33] rmarkdown_2.7     farver_2.0.3      purrr_0.3.4       scales_1.1.1     
## [37] codetools_0.2-18  htmltools_0.5.1.1 ellipsis_0.3.1    assertthat_0.2.1 
## [41] colorspace_2.0-0  tiff_0.1-7        labeling_0.4.2    stringi_1.5.3    
## [45] munsell_0.5.0     crayon_1.4.1