Landsat 8 Image Processing & Visualization

Natural Color Image:

p1<-ggRGB(s, r=4, g=3, b=2, stretch = "lin")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  ggtitle("Natural Color\n (R= Red, G= Green, B= Blue)")
p1

False Color Image (Vegetation):

p2<-ggRGB(s, r=5, g=4, b=3, stretch = "lin")+
theme(axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank())+
  ggtitle("False Color Vegetation\n (R= NIR, G= Red,  B=Green)")
p2

False Color Image (Urban):

p3<-ggRGB(s, r=7, g=6, b=4, stretch = "lin")+
theme(axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank())+
  ggtitle("False Color Urban\n (R= SWIR1, G= SWIR2,  B=Red)")

Let’s plot all maps side by side, for this, you have to use following function

grid.arrange(p1, p2, p3, nrow = 1)

Conversion of DN values at spectral radiance at the sensor

At-sensor spectral radiance (also called top-of-atmosphere radiance) is the amount of energy received by the sensor per second per steradian (sr) per square meter. To convert the processing raw digital number (DN) to spectral radiance we will use radCor() function in RStoolbox. There are four methods available for doing this conversion: rad, apref, sdos, dos, costz.

• rad - Radiance
• apref - Apparent reflectance (top-of-atmosphere reflectance)
• dos - Dark object subtraction following Chavez (1989)
• costz - Dark object subtraction following Chavez (1996)
• sdos - Simple dark object subtraction. Classical DOS, Lhaze must be estimated for each band separately.

The gain (G) and bias or intercept (I ) are different for each band and are provided in the Landsat header file. For Landsat 8, a linear equation is used to convert from DN to radiance (L):

Clipping

lsat_tra <- radCor(lsat,    # landsat stack
    metaData = metaData,      # metadata (MTL file) 
    method = "rad")           # convert DN to radiance
lsat_tra

## class       : RasterStack 
## dimensions  : 504, 541, 272664, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 318195, 334425, 4737735, 4752855  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :      B1_tra,      B2_tra,      B3_tra,      B4_tra,      B5_tra,      B6_tra,      B7_tra,      B9_tra,     B10_tra,     B11_tra 
## min values  : 43.66201600, 33.50591100, 20.96619600, 10.20694100,  1.76545000,  0.27992000,  0.07246402,  0.00000000,  7.56201760,  7.00657720 
## max values  : 176.0562830, 191.2132650, 195.7889600, 176.2474060, 137.5278500,  49.9815120,  17.4347412,   0.2496565,  10.1630962,   9.1691854

Conversion of spectral radiance to reflectance at sensor

When spectral radiance values are compared within and between sensors, variations occur frequently because of difference in sun-target-sensor configurations. Hence, these data are often converted to Top of Atmosphere (TOA) or at sensor reflectance. It eliminates effect of different solar zenith angles at the time of image acquisition, minimizes the solar radiance band differences occurring between different sensors and accounts for differences in the Earth-Sun distance.

lsat_tre<- radCor(lsat,        # landsat stack
    metaData = metaData,         # metadata 
    method = "apref")              # convert DN top-of-atmosphere reflectance
lsat_tre

## class       : RasterStack 
## dimensions  : 504, 541, 272664, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 318195, 334425, 4737735, 4752855  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :       B1_tre,       B2_tre,       B3_tre,       B4_tre,       B5_tre,       B6_tre,       B7_tre,       B9_tre,       B10_bt,       B11_bt 
## min values  : 1.036013e-01, 7.764161e-02, 5.272034e-02, 3.043955e-02, 8.603770e-03, 5.488612e-03, 4.212881e-03, 0.000000e+00, 1.036556e+02, 2.290230e+02 
## max values  : 4.177576e-01, 4.430942e-01, 4.923433e-01, 5.256014e-01, 6.702040e-01, 9.794354e-01, 1.000000e+00, 3.115158e-03, 1.106258e+02, 2.413258e+02

Atmospheric correction and Conversion of sensor reflectance to surface reflectance

TOA reflectance does not take into account any atmospheric effects and it is simply a ratio of the at-sensor radiance with the incoming radiance from the sun, where the distance and angle of the sensor to the target are considered.
It is assumed that any radiation originating from dark pixels is due to atmospheric haze and not the reflectance of the surface. The haze values are estimates of path radiance, which can be subtracted in a dark object subtraction (DOS) procedure. Atmospheric haze affects almost exclusively the visible wavelength range. Therefore, typically, you’d only want to estimate haze in blue, green and red bands, but occasionally also in the NIR band. The atmospheric correction methods (sdos, dos and costz) apply to the optical (solar) region of the spectrum and do not affect the thermal band. The dark object subtraction approach relies on the estimation of atmospheric haze based on dark pixels. Dark pixels are assumed to have zero reflectance, hence the name. It is then assumed further that any radiation originating from such dark pixels is due to atmospheric haze and not to the reflectance of the surface itself. For Dark Object Subtraction (DOS), we use the estimateHaze() function to estimate the digital number (DN) pixel value for dark objects in the visible wavelength range, and then use the values to convert sensor reflectance to surface reflectance.

For Dark Object Subtraction (DOS), first we have to estimate the digital number (DN) pixel value of dark objects for the visible wavelength range using estimateHaze() function, then use these values convert sensor reflectance to surface reflectance

# Automatic haze estimation
hazeDN <- estimateHaze(lsat, # raster object
    hazeBands = 2:5,               # haze band B2,B3, B4 and B5  
    plot = TRUE,                   # display hisgram and haze values
    darkProp=0.01)               # Proportion of pixels estimate haze

The following methods are available in radCor() function:

• dos - Dark object subtraction following Chavez (1989)
• costz - Dark object subtraction following Chavez (1996)
• sdos - Simple dark object subtraction. Classical DOS, Lhaze must be estimated for each band separately.

We will convert the DN to values to surface reflectance using using the function radCor() with estimated haze and the simple dark object subtraction, sdos, argument.

lsat_sref <- radCor(lsat,       # raster object  
    metaData = metaData,          # metadata 
    method = "sdos",                # Simple dark object subtraction
    hazeValues = hazeDN,          # haze value 
    hazeBands = 2:5,                # Bands corresponding to hazeValues
    atmosphere = "veryClear")     # Atmospheric characteristics
lsat_sref

## class       : RasterStack 
## dimensions  : 504, 541, 272664, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 318195, 334425, 4737735, 4752855  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :       B1_sre,       B2_sre,       B3_sre,       B4_sre,       B5_sre,       B6_sre,       B7_sre,      B9_sre,       B10_bt,       B11_bt 
## min values  : 1.077987e-01, 6.414577e-03, 2.357145e-03, 3.929587e-03, 0.000000e+00, 5.283789e-03, 3.948982e-03,           ?, 1.036556e+02, 2.290230e+02 
## max values  :    0.4346718,    0.3744582,    0.4517692,    0.5054899,    0.5213277,    0.9434545,    0.9501195,           ?,  110.6257801,  241.3257968

writeRaster(lsat_sref, filename=(paste0(dataFolder,".//Landsat8//Landsat8_PROCESSED//multi.tif")), bylayer=F, overwrite=TRUE, format="GTiff")

rm(list = ls())