# Processing landsat 8 to make NDVI in R

I want to do the NDVI with images landsat 8 in R, radiometric correction, radiance and reflectance, I make it from the models of the USGS (model here)

But taking the image (ND) for the calculation of the radiance and this for the calculation of the reflectance, it gives me errors that I then show pixel values within the image in a cloudless zone.

Evidently the values of the radiance are wrong the calculation of the NDVI gives wrong. The research has led me to know that it is something related to the 16 bits in which the image is. These are the pinxel values of the original image.

As I can solve this problem in R and that the radiance d values give me good to be able to do the NDVI and me between -1 and +1...

Something similar was treated in another post but was solved given that ENVI internally corrects these values without doing intermediate processes (Ref:Processing landsat 8 in ENVI)

values ND, Radiance, Reflectance

• I had the same problem with invalid reflectance values. Check your sin() function description, does it require degrees or radians? – Mr. Che Apr 7 '17 at 8:55

In the post you link, ENVI is not doing any internal stuff, and R reliably deals with conversions between the 16-bits in the image and the 64-bit internally in R.

The most likely reason for your problems are based around the sin(Ose). In R, it calculates in radians, while the formula requires degrees. As such, you need to make your formula go sin(Ose*Pi/180).

• The equations I use are these: radiancia = ((MLQcal)+Al) reflectancia <- ((MpQcal)+Ap)/sin(sun*(pi/180)) The error as I read this in the form of save the cut of the study area, where it is stored different to give consistent values, as I understand the values of ND should be 0-255. – Esteven Muriillo Apr 10 '17 at 15:10
• @EstevenMuriillo - Landsat 8 is a 16-bit sensor, not 8-bits. As such the ND / DN / digital number should be between 0 and 65535. – Mikkel Lydholm Rasmussen Apr 10 '17 at 15:27
• You have all the reason, although I still do not understand is because then the reflectance gives me negative values of -0.132897958670363 if I use the formula well even using it so, reflectance <- ((0.000020000 * Qcal) - (0.100000)) / (without (sun * (Pi / 180))) .. equals the same value .. In the post that you attach you mention the same error but to me it gives me in the whole image not in the limit values of it, as seen in the image of The air pixel values.... – Esteven Muriillo Apr 10 '17 at 16:14
• @EstevenMuriillo - reflectance <- (0.00002*DN)-0.1 – Mikkel Lydholm Rasmussen Apr 10 '17 at 16:52
• If, exactly that I do, above, place the code in case you want to go – Esteven Muriillo Apr 10 '17 at 18:10

You are over complicating the equation a bit. For the OLI sensor, two of the expected coefficients (multiplicative and additive rescaling) are constants and the entire equation can be simplified to (skipping the radiance step):

at.sensor.reflectance = (x * multiplicative.rescaling + additive.rescaling) / sin(sun.elev * (pi /180))

Where;
x = pixel value
multiplicative.rescaling = 0.00002
sun.elev = scene specific (in metadata)


As pointed out previously, the bit depth of the OLI sensor is 16-bit and not 8-bit (0-255). Given the constants, the expected range of pixel values is 0-65535.

Also, please do look at the math behind NDVI. It is a ratio so, the bit-depth of the NIR and VIS bands do not matter (16-bit or floating point). Where the index gets dicey is if the reflectance values range into the negative. If this does happen, it is over a very small range of pixels and it is valid to just bound these outliers to a [0-1] range.

Commonly, you correct to at sensor reflectance, before deriving NDVI, to remove atmospheric attenuation. Honestly, if you are not comparing NDVI through time or across disjunct scenes, this is not necessary and the measure of photosynthetically active radiation (PAR) should not be too biased.

If, exactly that I do, I calculate the radiance with the DNs then I calculate the reflectance with the radiance and still give me negative values of the reflectance.

metadato <- list.files(pattern = "txt")

REFLECTANCE_MULT_BAND = m[c(169:174),]
valor_REFLECTANCE_MULT_BAND <- substr(REFLECTANCE_MULT_BAND,31,41)

SUN_ELEVATION = m[68,]
valor_SUN_ELEVATION <- substr(SUN_ELEVATION,21,41)
#valor_SUN_ELEVATION

imagenes <- list.files(pattern = "Radiancia")
imagenes
b = lapply(imagenes[c(3:8)], FUN=raster)

#raciancia=(Mp*Qcal+Ap)/sun --Mp=radiance mult band -- Qcal=bandas -- Ap=radiance add band

for (i in 1:6) {
nom = names(b[[i]])
print(paste("la imagen usada sera: ",nom))
print(paste("Mp: ",valor_REFLECTANCE_MULT_BAND[i]," y Ap: ",valor_REFLECTANCE_ADD_BAND[i]))

Mp <-as.numeric(as.character(valor_REFLECTANCE_MULT_BAND[i]))
DN <- b[[i]]
sun <- as.numeric(as.character(valor_SUN_ELEVATION))
reflectancia <- ((0.000020000*DN)-0.100000)/(sin(sun*(pi/180)))

outname <- paste(substr(nom,1,24),sep='')
writeRaster(reflectancia, paste(outname,"_Reflectancia.tif"))