3

I'm trying to calculate the NDVI for a scene, and I've downloaded Landsat 8 OLI/TIRS Collection 2 Level 2 data, which means it has already been processed to surface reflectance.

The scaling factors for Collection 2, however, are new, and include an "additive" factor in addition to a multiplication factor. Following the specification then, I multiplied all the cell values in my scene by 0.0000275 and subtracted 0.2. Here is a simplified example (I'm working in R and using the raster and terra libraries):

library(raster)
library(terra)
library(tidyverse)

# Getting the names of all the Surface reflectance band files
filenames <- paste0("landsat/LC08_L2SP_048026_20200805_20200916_02_T1/LC08_L2SP_048026_20200805_20200916_02_T1_SR_B", 1:7,".TIF")

# Loading them as a raster stack
aug20 <- terra::rast(filenames)

# Setting band names
names(aug20) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

# Applying scaling factors (these operations are applied to every layer in the stack)
aug20 <- (aug20*0.0000275) - 0.2

I assumed that this manipulation should be enough but I ran into 2 pretty major issues. First, the resulting "surface reflectance" values have a range that far exceeds 0 and 1 on both sides:

# Summary of the stack, showing that all bands are beyond the 0 and 1 range
print(aug20)
> aug20
class       : SpatRaster 
dimensions  : 8051, 7951, 7  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 283485, 522015, 5292285, 5533815  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
source      : spat_jvKP7k9AhCUT5B3.tif 
names       : ultra-blue,       blue,      green,        red,        NIR,      SWIR1, ... 
min values  : -0.1999725, -0.1999725, -0.1999725, -0.1999725, -0.1156025, -0.0204800, ... 
max values  :   1.243200,   1.257198,   1.274935,   1.317257,   1.480910,   1.527193, ... 

To account for this I tried to set everything greater than 1 to 1 and everything less than 0 to 0, but I'm not sure if this is appropriate since I don't know how to interpret these out of range values.

Second, and more importantly, when I try to calculate the NDVI, the results are completely wrong for some pixels. For example, some areas that I know to be water have "unscaled" values that are as low as 7300-7500 across all bands. The image shows a sample of pixel values inside this known water area using QGIS (this image is a true-colour composite):

Highlighted areas known to be water and the band values for 1 pixel inside it

The issue is that since the scaling factor involves both multiplication and subtraction, the resulting SR values end up being extremely small for all bands over water, which in turn means that calculating an NDVI ratio between then returns a REALLY high value which signals dense vegetation (for pixels that I know to be water!)

# Using the sample cell values from above as an example, applying the rescaling factor:
b5_scaled <- (7414*0.0000275) - 0.2
[1] 0.003885
b4_scaled <- (7273*0.0000275) - 0.2
[1] 0.0000075

# NDVI for this pixel (WAY too high for water):
(b5_scaled - b4_scaled)/(b5_scaled + b4_scaled)
[1] 0.9961464

This is a simplified example to illustrate where the issue is, but it is affecting the scene overall where an NDVI calculated on the "scaled" data totally obscures water while calculating it on the "unscaled" data manages to capture it just fine:

# Plot of unscaled NDVI (left)
plot((aug20_unscaled$NIR - aug20_unscaled$red) / (aug20_unscaled$NIR + aug20_unscaled$red))

# Plot of scaled (right)
plot((aug20_scaled$NIR - aug20_scaled$red) / (aug20_scaled$NIR + aug20_scaled$red))

NDVI comparison

I'm working with several tiles, and this seems to be an issue with all of them so I don't think it is a feature of just this single tile. I'm really not sure what to do. I will likely be reverting to just using Collection 1 data for the time being but I'd still really like to know if I'm just missing something obvious or there is another way to resolve this.

2

Thanks for making the example file available. Clearly you need to apply the scale (gain) and offset values to get reflectance values (I first suggested that this may not be the case). With the example data I get the same as you do

library(terra)
filenames <- paste0("LC08_L2SP_048026_20200805_20200916_02_T1_SR_B", 1:7,".TIF")
aug_20 <- terra::rast(filenames)
names(aug_20) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')  
aug20 <- aug_20 * 0.0000275 - 0.2

But my map is different --- I get low values for water. Did your NDVI computation go wrong?

aug20 <- clamp(aug20, 0, 1)
ndvi <- function(red, nir) (nir - red) / (nir + red)
x <- lapp(aug20[[c("red", "NIR")]], ndvi)
plot(x)

enter image description here

The reason for getting high NDVI in the water is that these pixels have very low red and NIR values. You could consider them zero. So you could do

aug20 <- clamp(aug20, 0.02, 1)

And the problem goes away. But it is a stop-gap. The basic problem is that you get negative and very small numbers after adjusting.

If all you are doing is computing NDVI, then there is no reason for applying the scale/offset. Applying gain/offset only introduces noise in this case. Just use the raw data you get good NDVI values. (as you showed yourself).

y <- lapp(aug_20[[c("red", "NIR")]], ndvi)
plot(y)
4
  • This is interesting. I repeated your steps exactly and got the same map as you did (didn't know about clamp, thanks for that). However, when I crop this ndvi to my study area I see that the problem persists. The inland water features that I pointed to in the example above are still classified as dense vegetation. I've added the study-area shapefile to the folder, so you can see for yourself. I've also added the map I made following your steps and the same map cropped to the study area.
    – Soosh98
    Apr 13 at 14:48
  • Yes, I see that, the problem is with the very small values. Just do not rescale and you are good. Apr 14 at 3:54
  • It seems odd that applying the specified rescaling factor should alter the data this much... Isn't the point of rescaling just to optimize file-size/storage (since integers require less space than floating points)? Regardless, it helps to know that I can just run the NDVI without it. Thank you for all your help!
    – Soosh98
    Apr 14 at 13:45
  • I agree. But negative reflectance makes no sense either. Apr 14 at 19:40
1

Landsat atmospheric correction and surface reflectance retrieval algorithms are not ideal for water bodies due to the inherently low surface reflectance of water. Similarly, surface reflectance values greater than 1.0 can be encountered over bright targets such as snow and playas. These are known computational artifacts in the Landsat surface reflectance products.
Reference

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.