NDWI Calculation with Surface Reflectance Products

Question

I have downloaded L8 OLI/TIRS C1 Higher Level Data from USGS Earth Explorer to map Chlorophyll content in lake. The file contains original images and surface reflectance images. When I calculate NDWI using the original images, the lakes and surrounding area can be distinguised clearly. However, when I use the surface reflectance images, they produce a map of single colour. I have used the same NDWI formula for both the images. I have used the following code to calculate NDWI:

LS_ndwi <- (tmp3-tmp5)/(tmp3+tmp5) 
plot(LS_ndwi) 

Here tmp3 is the Green band and tmp5 is the NIR band. I wanted to know, should I use the original images or the surface reflectance images for my analysis? I have to retrieve the reflectance data from the images for my analysis, so wanted to use the surface reflectance imageries. I checked my data type and it is in "FLT8S". So the problem of integer is not valid here I think.

Answer 1

There is no real difference in applying the ratio to integer or float data. In evaluating results, you should first, be looking at the distributional characteristics of your data and not the visual. If the data is in fact "flat" you should see this in a plot of the probability density function (PDF).

Here is a worked example, using data in the RStoolbox package, that illustrates the expected results from NDWI. Please note that, although the USGS uses the Gao method, I prefer the McFeeters ratio of NIR and Green.

###############################################################
####  Normalized difference water index (NDWI) 
# parameters:
# nir - (630 to 690nm) landsat 5&7 band 4 
# swir - short-wave infrared 1 (1,550 to 1,750nm), landsat 5&7 band 5
# green - Green (2,090 to 2,350nm), landsat 5&7 band 2
# s – scaling factor (default is NULL)
# method – ("mcfeeters", "gao") Gao is the common method used 
#             by USGS in level 1 processing 
# scale – apply data scaling (default is FALSE) 
###############################################################
ndwi <- function(nir = NULL, swir = NULL, green = NULL, s = NULL, 
                 method = c("mcfeeters", "gao"), scale = FALSE) {
  if(is.null(nir)) stop("Must define NIR band")
    if(class(nir) != "RasterLayer")
      stop("Data must be raster class objects")  
      if(method[1] == "gao" ) {
        if(is.null(swir)) stop("Must define SWIR band")
          cat("Calculating NDWI using NIR and SWIR")      
        if(class(swir) != "RasterLayer")
          stop("Data must be raster class objects")  
        i <- (nir - swir) / (nir + swir)    
      } else if(method[1] == "mcfeeters" ) {
        if(is.null(green)) stop("Must define Green band")
          cat("Calculating NDWI using NIR and Green")     
        if(class(green) != "RasterLayer")
          stop("Data must be raster class objects")  
        i <- (green - nir) / (green + nir)  
      } else {
        stop("Not a supported method")
      }     
  if(scale == TRUE) {
   i <- (i - cellStats(i, min, asSample=FALSE)) * (1 - -1) /
        (cellStats(i, max, asSample=FALSE) - 
         cellStats(i, min, asSample=FALSE)) + -1
  } 
    if( !is.null(s) ) { i <- s * i }
  return( i )
}

Example with at-sensor reflectance correction (float) and without (8-bit integer).

library(raster)
library(RStoolbox)

data(lsat)

metaData <- readMeta(system.file("external/landsat/LT52240631988227CUB02_MTL.txt",package="RStoolbox"))
lsat.ref <- radCor(lsat, metaData = metaData, method = "apref")

ndwi.int <- ndwi(nir=lsat[[4]] , green=lsat[[2]])
ndwi.float <- ndwi(nir=lsat.ref[[4]] , green=lsat.ref[[2]])

par(mfrow=c(2,2))
  plot(density(ndwi.int[]),main="distribution of ndwi int")
  plot(ndwi.int, main="ndwi from integer")
  plot(density(ndwi.float[]),main="distribution of ndwi float")
  plot(ndwi.float, main="ndwi from float")

As you can see there is not a notable difference in the index or the resulting distributions. But there is some variation. We can see the specific differences in the distributional moments.

summary(ndwi.int[])
summary(ndwi.float[])

This can likely be attributed to the fact that at-sensor reflectance for level 1 processing also includes dark-object subtraction and changes atmospheric attenuation thus, changing some the characteristics of the data. This is why it is at-sensor reflectance recommended for calculating metrics such as this.

I am not sure why you are seeing a uniform color but, my guess would be a stretch issue in the plot function. You may want to try RStoolbox::ggR plot function, with stretch = "lin", "sqrt" or "log", to see if you get the same result. Also, first check the distribution of the data. If it is flat then there is, in fact, something going sideways in the metric. I included a function for NDWI so, you may want to try it but, that said, your calculation is correct.

Answer 2

I applied the following to successfully create a NDWI GeoTiff after I clipped a region. You need to clip a region otherwise the processing area might be too large and you'll run into memory issues.

library(raster)
library(rgdal)

#Load the landsat
landsat8 = brick("C:/Test.tif")

# Review the raster dimensions and print details for inspection
dim(landsat8)
print(landsat8)

# plot RGB
plotRGB(landsat8, r=4, g=4, b=3, axes = TRUE, stretch="lin")

#Calculate NDWI
band3 = raster(landsat8, layer = 3)

band5 = raster(landsat8, layer = 5)

NDWI = (band3-band5)/(band3+band5)

plot(NDWI)

#Write the GeoTiff for use in QGIS or ArcGIS
writeRaster(NDWI,'C:/NDWI.tif', options = c('TFW=YES'))

The following provides guidance on "surface reflectance": https://lta.cr.usgs.gov/L8Level2SR https://landsat.usgs.gov/sites/default/files/documents/lasrc_product_guide.pdf and https://landsat.usgs.gov/sites/default/files/documents/si_product_guide.pdf

There are no issues in calculating NDVI or NDMI indices using reflectance data according to the final guide linked above. However, you should do some research on your methods and terminology so that we can understand what you mean by "original images". The following paper can be used as a guide to communicate more clearly what you are doing:

Young et al. (2017). A survival guide to Landsat preprocessing. Ecology, 98(4):920–932

Answer 3

NDMI is calculated using the near-infrared (NIR) and the short-wave infrared (SWIR) reflectance:

NDMI = (NIR – SWIR) / (NIR + SWIR)

Sentinel-2: NDMI = (B08 – B11)/(B08 + B11)

Landsat 4-5 TM: NDMI = (B04 – B05) / (B04 + B05)

Landsat 7 ETM+: NDMI = (B04 – B05) / (B04 + B05)

Landsat 8: NDMI = (B05 – B06) / (B05 + B06)

MODIS: NDMI = (B02 – B06) / (B02 + B06)

I found this article here: https://eos.com/make-an-analysis/ndwi/