Normalized Differential Vegetation Index in Mexico - Remote Sensing Analysis

This project shows the analysis of Landsat 8 satellite images of Mexico to calculate the NDVI. It includes preliminary raster manipulation from raster files I created using Google Earth Enginle.

Carmen Galaz-García true
05-28-2021

About the Data

I generated the raster file mexico_median.tif using Google Earth Engine. It is available in my GitHub repository for this project. This raster shows the median of the images in the USGS Landsat 8 Collection 1 (Tier 1 TOA Reflectance) from May 1st 2013 to May 1st 2014. I will probably write another blog post about how I used Google Earth Engine! Meanwhile, here’s a screenshot of the Java script.

Google Earth Engine script to create Mexico RGB raster file

Raster exploration

The Landsat 8 raster I produced using Google Earth Engine is a coarse approximation of Mexico’s shape. It is calculated using all the Landsat 8 scenes that intersect a coarse polygon enclosing Mexico. We can see here the red, green and blue layers (RGB).

hide 
# -----  EXPLORING RGB LAYERS -----
# read all bands as raster stack
mex_rgb <- raster::stack(here('_posts/2021-05-28-mexico-rasters/mexico_median_rgb.tif')) 
#mex_rgb
# looking at stack information we get that:
#   B4 band = red
#   B3 band = green
#   B2 band = blue

# ---- Define Plot Function ----
color_graph <- function(band,color){
  return (plot(band, 
               col = hcl.colors(n = 100, palette = color),
               axes = FALSE,
               legend = FALSE)
          )
}
# --- Colors
RGB = c('Reds 2', 'Greens 2', 'Blues 2')

# --- plot RGB bands
for (i in 1:3){
  color_graph(mex_rgb[[i]],RGB[i])
}



Masking

To do the masking and select only the region corresponding to Mexico I used a polygon shapefile made by Hijmans, Robert J. from the Museum of Vertebrate Zoology at UC Berkeley. To use the polygon shapefile for masking we need to convert it into a raster. Then we can use it mask each layer of the raster stack.

hide 
# --- read shapefile
mex_sf <- read_sf(here('_posts/2021-05-28-mexico-rasters/stanford-zc863pb5331-shapefile/zc863pb5331.shp')) %>%
  st_transform(crs(mex_rgb))   # match CRS

# --- transform shapefile into raster using red band as guide
mex_rast <- fasterize::fasterize(mex_sf, mex_rgb[[1]])
#writeRaster(mex_rast, 'data/county.tif') #export mask raster
plot(mex_rast,axes=FALSE, box=FALSE)

hide 
# --- mask each layer
for (i in 1:3){
  mex_rgb[[i]] <- mask(mex_rgb[[i]],mex_rast)
}
rm(mex_rast)

# --- plot masked RGB bands
for (i in 1:3){
  color_graph(mex_rgb[[i]],RGB[i])
}



True Color

Now that we have the masked red, green and blue bands, we can use these to make a true color composite image of Mexico. It is a close replica to what we can see with our eyes. Notice too the storm around Tabasco. This is slightly inconvenient because it doesn’t allow to have a clear view of the ground.

hide 
plotRGB(mex_rgb, scale=0.63, colNA='skyblue1')



Normalized Differential Vegetation Index (NDVI)

Near-Infrared Band

To calculate the NDVI we need the near-infrared band too. Since the previous raster only contained RGB bands, here I added a raster I generated with Google Earth Engine that shows the median of the near-infrared band from May 1, 2013 to May 1, 2014.

hide 
# --- read nir raster
nir <- raster(here('_posts/2021-05-28-mexico-rasters/mexico_median_nir.tif')) 

# --- masking (note: had to update mask to match extent of nir raster)
mex_rast2 <- fasterize::fasterize(mex_sf, nir)
nir <- mask(nir, mex_rast2)
rm(mex_rast2)

# --- plot
df_nir <- nir %>% 
  raster::rasterToPoints() %>%
   as.data.frame()

ggplot() +
  geom_tile(data = df_nir, aes(x = x, y = y, fill = mexico_median_nir)) +
  coord_sf(expand = 0) +
  scale_fill_gradient(low = 'goldenrod1', high = 'red') +
  labs(fill='NIR')+
  theme_void()
hide 
rm(df_nir)

NDVI calculation

According to the Earth Observing System

The results of the NDVI calculation range from -1 to 1. Negative values correspond to areas with water surfaces, manmade structures, rocks, clouds, snow; bare soil usually falls within 0.1- 0.2 range; and plants will always have positive values between 0.2 and 1. Healthy, dense vegetation canopy should be above 0.5, and sparse vegetation will most likely fall within 0.2 to 0.5.

The NDVI is calculated using the near-infrared and red bands of the satellite image. The formula is \[NDVI = \frac{NIR - Red}{NIR + Red}.\] In this section we create a map that outlines zones in Mexico with healthy and heavy vegetation canopy (\(0.5\leq NDVI\)), with code based on this repository by C. O’Hara.

hide 
# ---- NDVI calculation ----
ndvi <- (nir - mex_rgb[[1]]) / (nir + mex_rgb[[1]]) 


# --- Dense Vegetation Function ---
Fdense_veg <- function(x) {
  return(ifelse(0.5<=x, 1, NA))
}
# --- dense vegetation data frame
df_dense <- calc(ndvi, fun = Fdense_veg) %>% 
   raster::rasterToPoints() %>%
   as.data.frame()
# --- ndvi data frame
df_ndvi <- ndvi %>% 
  raster::rasterToPoints() %>%
  as.data.frame()

# --- create plot 
ggplot() +
  geom_tile(data = df_ndvi, aes(x = x, y = y, fill = layer)) +
  geom_tile(data = df_dense, aes(x = x, y = y),fill = 'green4') +
  coord_sf(expand = 0) +
  scale_fill_gradient(low = 'black', high = 'white') +
  labs(fill='NDVI')+
  theme_void()

Fig 1. Zones in Mexico with median healthy and heavy vegetation canopy (\(0.5\leq NDVI\)) in color green, from 05/01/2013 to 05/01/2014 mapped over an NDVI map of the country.

Vegetation density coverage

We can also use the raster data to compute the percentage of the land covered by vegetation of a certain density.

hide 
# --- Count # of pixels within each vegetation range
veg_thresh <- c(0.2, 0.45, 0.7)
c_sparse <- sum(veg_thresh[1]<=ndvi[] & ndvi[] <= veg_thresh[2], na.rm=TRUE)
c_mod <- sum(veg_thresh[2]<=ndvi[] & ndvi[] <= veg_thresh[3], na.rm=TRUE)
c_dense <- sum(veg_thresh[3]<=ndvi[] & ndvi[] <= 1, na.rm=TRUE)

c_total <- sum(is.na(ndvi[]) != TRUE)
c_nonveg <- c_total - (c_dense+c_mod+c_sparse)


# --- Create data set
veg_pixels <- data.frame(
  category=c("Sparse", "Moderate", "Dense", "No vegetation"),
  count=c(c_sparse, c_mod, c_dense, c_nonveg)
)
# -- add info for ring chart
# Compute percentages
veg_pixels$fraction = veg_pixels$count / sum(veg_pixels$count)
# Compute the cumulative percentages (top of each rectangle)
veg_pixels$ymax = cumsum(veg_pixels$fraction)
# Compute the bottom of each rectangle
veg_pixels$ymin = c(0, head(veg_pixels$ymax, n=-1))
# Compute label position
veg_pixels$labelPosition <- (veg_pixels$ymax + veg_pixels$ymin) / 2
# Compute a good label
veg_pixels$label <- paste0(veg_pixels$category, "\n value: ", round(veg_pixels$fraction*100,2), "%")

 
# ---- Make  plot ----
ggplot(veg_pixels, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=category)) +
     geom_rect() +
     coord_polar(theta="y") +
    geom_label( x=3.5, aes(y=labelPosition, label=label), size=3) +
    scale_fill_manual(values = c("darkgreen", "forestgreen", "gray70", "lightgreen"))+
    xlim(c(2, 4)) +
    theme_void() +
  theme(legend.position = "none")

Figure 2. Percentage of land in Mexico covered by dense vegetation (\(0.7\leq NDVI\)), moderately dense vegetation (\(0.45\leq NDVI < 0.7\)) and sparse vegetation (\(0.2 \leq NDVI <0.45\)), based on data of median satellite images (05/01/2013 to 05/01/2014).

References

Boundary, Mexico, 2015. Hijmans, Robert J. and University of California, Berkeley. Museum of Vertebrate Zoology.
https://geodata.lib.berkeley.edu/catalog/stanford-zc863pb5331

USGS Landsat 8 Collection 1 Tier 1 TOA Reflectance. Dataset Provider: USGS/Google. https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_TOA

GitHub repository: esm244_w2021_lab6_rasters C. O’Hara, 2021. https://github.com/oharac/esm244_w2021_lab6_rasters