Raster_NDVI

Raster data manipulation with terra

terra is a recently developed package that’s becoming very popular. It can handle both vector and raster data. The package is especially powerful for raster data manipulation.

SpatRaster

We discussed in class that photographs are rasters. They are single layered rasters, in the sense that one pixel only contains information about the color represented in the pixel.

Raster data used in urban planning are more complex. Take Satellite data for example. Consider a satellite that can capture data from Earth at a 5 m resolution. It means that every 5 m X 5 m unit in the Earth surface is one pixel of that image. Satellites capture different kinds of information for the same unit, for example satellite sensor readings on infrared waves, microwaves, visual spectrum (red, green, blue) etc. Collectively, each different kind of information is a different layer; leading to a multilayered raster data.

SpatRaster format in terra represents multi-layered raster data.

To understand raster data in depth, let’s create some ourselves.

library(terra)
terra 1.7.29
library(repr)
r <- rast(ncol=10, nrow=10, xmin=-150, xmax=-80, ymin=20, ymax=60)
r
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84

  • rast(): This is the function used to create the raster object. It initializes a new raster with the specified properties.

  • ncol=10, nrow=10: These are the arguments to the rast() function, specifying the number of columns and rows of the raster. In this case, the raster will have 10 columns and 10 rows.

  • xmin=-150, xmax=-80, ymin=20, ymax=60: These are the arguments to the rast() function, specifying the extent or spatial coverage of the raster. The xmin and xmax values represent the minimum and maximum x-coordinate (longitude) values, while the ymin and ymax values represent the minimum and maximum y-coordinate (latitude) values. In this case, the raster will cover the geographic extent from -150 to -80 in longitude and from 20 to 60 in latitude.

By executing the code, a new raster object will be created with the specified dimensions (10x10) and spatial extent (-150 to -80 in longitude, 20 to 60 in latitude). The variable r will hold this newly created raster object, which can be further manipulated and analyzed using appropriate raster functions and operations in R.

The object r is a raster geometry without any value.

plot(r)

values(r) <- runif(ncell(r))
r
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source(s)   : memory
name        :       lyr.1 
min value   : 0.004707775 
max value   : 0.987785708

By executing the code, random values will be generated using the runif() function, and these values will be assigned to the cells of the raster object r. Each cell in the raster will have a random value drawn from a uniform distribution. This is a common technique for assigning random values to raster cells for simulation or testing purposes.

Finally, we have some values, albeit random, to plot.

plot(r)

We could also assign cell numbers (in this case overwriting the previous values).

values(r) <- 1:ncell(r)
r
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source(s)   : memory
name        : lyr.1 
min value   :     1 
max value   :   100 
plot(r)

This creates this spectrum because the value is equal to the cell number (counted like you would read a comic book: from left to right, then top to bottom).

Exercise
  • Create a 5x5 raster with normally distributed values between 1 and 20.
  • Hint: Look up rnorm().

Creating multiple layers

r2 <- r * r
r3  <- sqrt(r)
s <- c(r, r2, r3)
s
class       : SpatRaster 
dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source(s)   : memory
names       : lyr.1, lyr.1, lyr.1 
min values  :     1,     1,     1 
max values  :   100, 10000,    10 
plot(s)

We see that s is a multilayered object.

We can access individual layer with [[]] notation and their numbering in the layers.

plot(s[[3]])

Exercise
  • To s, add another layer that is the sum of values in r2 and r3
  • Plot that layer only by indexing from s
  • Look up naming the layers in terra, implement it to add names to these layers, and access a layer of s using its name

Indexing vegetation in satellite images

We can create an index of presence of vegetation in a given area from satellite images using NDVI (Normalized Difference Vegetation Index). NDVI has significant applications in various fields, including urban planning.

NDVI is a vegetation index that quantifies the presence and vigor of vegetation based on the reflectance of different wavelengths of light. It is calculated using the following formula:

NDVI = (NIR - Red) / (NIR + Red)

where NIR represents near-infrared reflectance and Red represents red reflectance. NDVI values range from -1 to +1, with higher values indicating denser and healthier vegetation.

Vegetation cover is related to air quality, well being, surface temperature, and other issues that urban planners are interested in.

In this example, we calculate NDVI using Landsat 8 satellite data. Landsat 8 is a satellite mission operated by the United States Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA). It is part of the Landsat program, which has been capturing Earth observation data since 1972. Landsat 8 was launched in 2013 and continues to provide valuable imagery for various applications, including environmental monitoring, land cover classification, and natural resource management.

Landsat 8 captures data in the visible, near-infrared, and thermal infrared portions of the electromagnetic spectrum. It has a set of different spectral bands, each capturing light at specific wavelengths.

Importing data

First we import tif files for Red and NearInfrared bands. TIF (Tagged Image File Format) is used to store raster data.

redFile <- './Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B4_clip.tif'
nirFile <- './Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B5_clip.tif'
redBand <- rast(redFile)
nirBand <- rast(nirFile)
plot(redBand)

plot(nirBand)

red <- redBand[[1]] # we do this to extract the first band of the image. Both images have only one band, so we could alternatively just use redBand and nirBand
nir <- nirBand[[1]]

Function to estimate NDVI

ndviCal <- function(red, nir) {

    ndviArray <- (nir - red)/(nir + red)
    return(ndviArray)
}

In the following chunk, we calculate NDVI.

Notice the use of Sys.time() before and after the calculation. It’s a programming practice called benchmarking. This helps us estimate how long some processes take. We can also use benchmarking to find out alternatives that process faster.

iniTime <- Sys.time()
#ndviBand <- ndviCal(redBand,nirBand)
ndviBand <- ndviCal(red,nir)
endTime <- Sys.time()
print(endTime - iniTime)
Time difference of 0.142606 secs

We will plot the NDVI highlighting areas with high vegetation.

options(repr.plot.width=20, repr.plot.height=12)
plot(ndviBand, main="Landsat-NDVI")

We can also save the raster as tif files.

writeRaster(ndviBand, './Landsat8/NDVIByR.tif', overwrite=TRUE)

Distribution of values

hist(ndviBand)
Warning: [hist] a sample of35% of the cells was used

Categorization of raster

# Define the cutoff value
cutoff <- 0.2

# Categorize the raster into "green" or "not green"
green_raster <- ifel(ndviBand >= cutoff, as.logical(1), as.logical(0))

# Visualize the categorized raster
plot(green_raster, main = "Categorized Raster")

The cutoff value is somewhat arbitrary. Why this cutoff? The NDVI capitalizes on the contrast between characteristics of two types of radiation: the chlorophyll pigment absorption in the Red region, and the high reflectivity of plant material in the Near Infrared (NIR) region. When the NDVI is above 0.2 or 0.3, it’s a good indicator of green vegetation since it suggests a higher degree of reflectivity in the NIR region and absorption in the Red region – both characteristics of healthy, green vegetation.

Bear in mind that these are just general guidelines and can vary depending on the specific type of vegetation and the environmental conditions. For more accurate results, it may be useful to ground truth the data or use additional remote sensing indices and/or data sources.

Reclassification of raster

m <- c(-1,0.25, 0.3, 0.4, 0.5, 1)
vegc <- classify(ndviBand, m)
plot(vegc, col = rev(terrain.colors(4)), main = 'NDVI based thresholding')

Exercise
  • Categorize the raster into 3 groups, 0-0.2, 0.2-0.4, and 0.4 to 1 and visualize it using plot.
  • Benchmark the time taken

References