Masking Landsat 9 images using a shapefile of an ROI

Introduction

This is a brief notebook to use a stacked Landsat 9 image to extract an ROI. In my case, I will be using the county of Kiambu in Kenya. This is because the county is perfectly inside one scene, this offers ease in extracting the ROI. I will come up with another notebook to detail how to extract the ROI of an area where the scenes overlap.

Importing the files

I will be using the previously stacked image to get my area extract LC09_L2SP_STACK_1_7.TIF. I will utilise three libraries:

• tidyverse: This is primarily for various data manipulation since the raster stack is just another data-frame.

• terra: This is obviously for the raster manipulation.

• sf: This is for vector data manipulation. In our case, a shapefile.

library(tidyverse)
library(plotly)
library(terra)
library(sf)

imgDir = "C:/Users/roywa/Data/LST/LC09_L2SP_168061_20230119_20230313_02_T1/"
imgName = "LC09_L2SP_STACK_1_7.TIF"
imgPath = paste0(imgDir, imgName)

kiambuFilePath = "C:/Users/roywa/Data/SHP/Kiambu/"

imgData <- rast(imgPath)
kiambu <- st_read(kiambuFilePath, drivers="ESRI Shapefile")

options:        ESRI Shapefile 
Reading layer `Kiambu' from data source `C:\Users\roywa\Data\SHP\Kiambu' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 36.49081 ymin: -1.312168 xmax: 37.3626 ymax: -0.7566282
Geodetic CRS:  WGS 84

print("IMAGE SUMMARY")

[1] "IMAGE SUMMARY"

print(imgData)

class       : SpatRaster 
dimensions  : 7731, 7591, 7  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 175185, 402915, -276315, -44385  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 37N (EPSG:32637) 
source      : LC09_L2SP_STACK_1_7.TIF 
names       : BAND 1, BAND 2, BAND 3, BAND 4, BAND 5, BAND 6, ... 
min values  :      1,      5,    130,     74,   6152,   7554, ... 
max values  :  42745,  42834,  41997,  42694,  44407,  52291, ... 

Observing the data

Looking at the data we notice a few things:

1. The Image has bands so that means that it is a stacked image. The one we created in the previous notebook. Additionally we can check this with the names() function

imgData %>%
  names()

[1] "BAND 1" "BAND 2" "BAND 3" "BAND 4" "BAND 5" "BAND 6" "BAND 7"

2. The two datasets are of different CRS; this means that we cannot perform any operation that involves the two together

crs(kiambu, proj=TRUE)

[1] "+proj=longlat +datum=WGS84 +no_defs"

crs(imgData, proj=TRUE)

[1] "+proj=utm +zone=37 +datum=WGS84 +units=m +no_defs"

Re-projection

The natural next step is to re-project these datasets so that they conform to each other and therefore operations can be performed on them. In this case of the data, I would rather re-project the vector dataset than the image. The image has 7 layers that would require warping and the image takes a whooping 700MB on disk.

kiambu_utm <- kiambu %>%
  st_transform(crs(imgData))

par(mfrow=c(1,2))
vect(kiambu)%>%
  plot(main="Kiambu latlong", col="red")
vect(kiambu_utm) %>%
  plot(main="Kiambu UTM", col="blue")

In the plot above we can see the red plot being the WGS84 plot with latlong graduation while the blue one is the one that we aim to use since its in the same projection

crs(kiambu_utm) == crs(imgData)

[1] TRUE

For the image that I want to mask:

plotRGB(c(
  imgData$`BAND 4`,
  imgData$`BAND 3`,
  imgData$`BAND 2`
))

Masking

So now I have an imgData that for the sake of understanding I will use the 4th band for visualizing what is happening.

plot(imgData$`BAND 4`)
lines(vect(kiambu_utm))

This is proof that these two objects are in fact in the same CRS and that operations can be performed on both of them.

kiambu_subset_name = "KIAMBU_LC09_L2SP_STACK_1_7.TIF"
kiambu_crop <- terra::crop(imgData, kiambu_utm)
kiambu_mask <- terra::mask(kiambu_crop, kiambu_utm, 
                           inverse=FALSE,
                           overwrite=TRUE, 
                           filename=paste0(imgDir,kiambu_subset_name))

kiambu_mask %>%
  plot()

plotRGB(c(
  kiambu_mask$`BAND 4`,
  kiambu_mask$`BAND 3`,
  kiambu_mask$`BAND 2`
))

And that is how you get a mask of a stacked image. The same method can be used for an individual image, however, I don’t fancy that method due to the loops that you have to implement to get the job done.