Step 1: Setup and Loading Libraries

Open RStudio and create a new script. Begin by loading the necessary libraries.

# Install and load required libraries

install.packages("raster")

library(raster) # For raster data manipulation

library(rgdal) # For working with vector data

library(sf) # For spatial data operations

library(terra) # For handling large raster datasets efficiently

Step 2: Setting Working Directory and Loading Data

Ensure your working directory is set to the folder containing your Landsat data. Then,
load the Landsat bands into separate raster objects.

# Set working directory

setwd("E:/download")

# Load Landsat bands

f1<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B1.tif"

B1 <- raster(f1)

# Plot Band 1

plot(B1)

f2<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B2.tif"

B2 <- raster(f2)

plot(B2)
f3<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B3.tif"

B3 <- raster(f3)

plot(B3)

f4<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B4.tif"

B4 <- raster(f4)

plot(B4)

f5<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B5.tif"

B5 <- raster(f5)

plot(B5)

f6<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B6.tif"

B6 <- raster(f6)

plot(B6)

f7<-"E:/download/reprojected/reproj_RT_LC09_L1TP_168065_20231002_20231002_02_T1_B7.tif"

B7 <- raster(f7)

plot(B7)

Step 3: Plotting Multiple Bands

Plot multiple bands in a single layout for comparison.

# Set layout for side-by-side plotting

layout(t(1:2))
# Plot Bands 1 and 2

plot(B1, main ="Band 1")

plot(B2, main ="Band 2")

# Open a new plotting device

windows()

# Set layout for side-by-side plotting

layout(t(1:2))

# Plot Bands 3 and 4

plot(B3, main ="Band 3")

plot(B4, main ="Band 4")

Step 4: Combining Bands into a Brick Object

Combine all bands into a single brick object and examine its properties.

# Combine all bands into a brick object

?brick

?stack

landsat_brick <- brick(B1, B2, B3, B4, B5, B6, B7)

# Display information about the brick

names(landsat_brick)

res(landsat_brick)

nlayers(landsat_brick)
xres(landsat_brick)

res(landsat_brick)

yres(landsat_brick)

Step 5: Displaying Color Composite Images

Visualize Landsat images in different color composites.

# Display a true color composite

?plotRGB

plotRGB(landsat_brick)

# Display a true color composite

plotRGB(landsat_brick, r=4, g=3, b=2, stretch = "lin")

Step 5: Saving Raster Images

Save the combined Landsat image and subsets.

#checking the current working directory

getwd()

# Save the bricked Landsat image

?writeRaster

writeRaster(landsat_brick, filename = "brickedlandsat", format = "GTiff", overwrite=TRUE)

Step 6: Subsetting via Defining Extension

Subset the raster data by defining an extent.

# Read the stacked Landsat data

# Display information about the stacked data

stacked

# Plot the second layer of the stacked data

plot(stacked[[2]])

# Define the extent

e = extent(35, 35.8, -7.8, -7.2)

# Plot the extent on top of the second layer

plot(e, add = TRUE)

# Crop the stacked data using the defined extent

s_stacked <- crop(stacked, e)

# Display information about the cropped data

s_stacked

# Plot the cropped RGB image

plotRGB(s_stacked, 4, 3, 2, stretch = 'lin')

Step 7: Subsetting via Shapefile

Subset the raster data using a shapefile as a mask.

# Read the shapefile

AOI <- readOGR("E:/download/AOI.shp")

# Display information about the shapefile

# Plot the stacked data with a true color composite

plotRGB(stacked, r = 4, g = 3, b = 2, stretch = "lin")

# Plot the shapefile boundary

plot(AOI, add = TRUE, border = 'black')

# Crop the stacked data using the shapefile

c_stacked <- crop(stacked, AOI)

# Display information about the cropped data

c_stacked

# Plot the cropped RGB image

plotRGB(c_stacked, 4, 3, 2, stretch = 'lin')

# Plot the shapefile boundary on top of the cropped image

plot(AOI, add = TRUE, border = 'black', lwd = 3)

# Mask the cropped data with the shapefile

c1_stacked <- mask(c_stacked, AOI)

# Display information about the masked data

c1_stacked

# Plot the masked RGB image

plotRGB(c1_stacked, 4, 3, 2, stretch = 'lin')

# Write the masked data to a new raster file

Step 8: Deriving NDVI, Plotting, and Saving

Calculate the Normalized Difference Vegetation Index (NDVI) from the Landsat data,
plot it, and save it as a raster file.

# Read the subset of Landsat data

AOI_Landsat <- stack("E:/download/AOI_Landsat.tif")

# Calculate NDVI

ndvi <- (AOI_Landsat[[5]] - AOI_Landsat[[4]]) / (AOI_Landsat[[5]] + AOI_Landsat[[4]])

# Display NDVI values

ndvi

# Plot NDVI

plot(ndvi)

# Save NDVI as a raster file

writeRaster(ndvi, filename = "ndvi_crop", format = "GTiff", overwrite = TRUE)

Step 9: Land Use and Land Cover Classification

Perform a classification based on the NDVI values to identify different land use and
land cover types.

# Set up plotting layout

par(mfrow = c(3, 2))

# Plot NDVI

plot(ndvi)
# Plot RGB image

plotRGB(AOI_Landsat, 4, 3, 2, stretch = 'lin')

# Define land cover classes based on NDVI thresholds

forest <- ndvi >= 0.3

grassland <- ndvi >= 0.25 & ndvi < 0.3

bushland <- ndvi > 0 & ndvi < 0.25

water <- ndvi <= 0

# Plot land cover classes

plot(forest)

plot(grassland, col = c("white", "yellow3"))

plot(bushland, col = c("white", "azure3"))

plot(water, col = c("white", "blue3"))

# Reset plotting layout

par(mfrow = c(1, 1))

Step 10: Plotting Land Use and Land Cover Classification

Plot the final land use and land cover classification with a legend.

# Combine land cover classes

LULC <- forest + grassland * 2 + bushland * 3 + water * 4

# Plot land use and land cover classification

plot(LULC, col = c("#079408", "#DDDC04", "#D3D5B3", "#0538CF"), legend = FALSE)

legend("topright",

legend = c("Forest", "Grassland", "Bushland", "Waterbody"),

col = c("#079408", "#DDDC04", "#D3D5B3", "#0538CF"),

text.font = 2,

pch = c(17, 19),

bty = "o",

pt.cex = 2,

cex = 1.2,

text.col = "black",

horiz = FALSE,

inset = c(0.02, 0.1))