Creating an NDVI of a raster

Introduction

One of the simplest and common operations to be performed on raster objects is computation. This is used to extract specific details about the dataset by using the various aspects of the bands.

In this notebook, I will explore how to perform this computations, which is fairly is easy. I will be using the Normalized Difference Vegetation Index (NDVI). This is a basic index for getting the vegetative cover of anb area by using the surface reflectance of Near Infrared Band and Red Band that is band 5 and 4 respectively.

NDVI

The NDVI can be expressed in a mathematical formula that is :

\[ ndvi = (nir-red)/(nir+red) \]

Libraries

As established the three packages required are terra, tidyverse and sf.

library(terra)
library(tidyverse)

── Attaching core tidyverse packages ─────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.2     ── Conflicts ───────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks terra::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(sf)

Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

sceneDir <- "C:/Users/roywa/Data/LST/LC09_L2SP_168061_20230204_20230311_02_T1/NORMALISED/"
kiambu_stack <- rast(paste0(sceneDir, "KIAMBU_SUBSET_STACK.TIF"))

This is a plot to see how, if at all, the data has been loaded. Though it is unnecessary. I prefer the assurance that is the plot rendering without any skirmishes.

kiambu_stack$`B3` %>%
  plot()

NDVI function

This is a function that I will use to perform the NDVI calculations. I prefere this method since the function can be reused multiple times without error or failure.

calc_ndvi <- function(nir, red){
  ndvi <- (nir-red)/(nir+red)
  return(ndvi)
}

For instance: This here is the performance of the NDVI for the image subset. The function returns a raster with the values defined.

kiambu_ndvi <- lapp(c(kiambu_stack$B5, kiambu_stack$B4), fun=calc_ndvi)

In my case, I have used the data to plot the raster as well as a histogram to see distribution of the various data points.

par(mfrow=c(1,2))
kiambu_ndvi%>%
  plot(main="NDVI OF KIAMBU COUNTY")

kiambu_ndvi%>%
  hist(main = "NDVI values", xlab = "NDVI", ylab= "Frequency",
    col = "wheat", xlim = c(-0.5, 1),  breaks = 30, xaxt = "n"
    )%>%
  axis(side=1, at = seq(-0.6, 1, 0.2), labels = seq(-0.6, 1, 0.2))

This is a threshold cut. I am taking in the values above 0.4 only while the rest are interprated as Null values. This is benefitial if you want to show a specific phenomena that occurs above a certain range of the dataset.

kiambu_ndvi %>%
  clamp(0.4, values=FALSE)%>%
  plot()

Alternatively you can select a spefific range by changing the ranges into a matrix.

m <- c(-Inf, 0.25, NA,  0.25, 0.4, 1,  0.4, Inf, NA)
rcl <- matrix(m, ncol=3, byrow=TRUE)
land <- classify(kiambu_ndvi, rcl)
land %>%
  plot(col="brown")

Finally, the data can also be classified using a vector.

m <- c(-1,0.25, 0.3, 0.4, 0.6, 1)
vegc <- classify(kiambu_ndvi, m)
plot(vegc, col = c("red","brown","orange", "yellow", "green"), main = 'NDVI based thresholding')

Conclusion

This is a very small subset of the indices that are used with remotely sensed images. There are a number of indiced that are way more complicated than the NDVI. However, all the indices that you choose should have a theoretical backing so that you understand eh innner working of these indices.