Part I

#to avoid R not recognising some files
load("E:/Documents/R Studio is here/R Studio files projects/gis_coding/.RData")

Doing GIS or Remote sensing work without touching satellite data will be a deservice to the term “remote sensing”. Remote sensing, is, in a nutshell, the perception and collection of data without coming into physical contact with it. For example, your eyes are a perfect example of remote sensors. They can perceive color, size or sometimes texture without directly contacting the object. Why would they even be used to contact it anyway? ;)

For this purpose, we downloaded the following Sentinel-2 images which will be used in this exercise. The images were downloaded from Earth Explorer, which also hosts Sentinel 2 data. For a small tutorial on how to download images from Earth Explorer, see my tutorial here: - https://www.researchgate.net/publication/341343109_BASIC_RASTER_AND_VECTOR_OPERATIONS_USING_QGIS_A_TUTORIAL

Below are the dates for the Sentinel-2 images downloaded from Earth Explorer for this tutorial. Use the address term, “Nakuru”.

Sentinel 2A Level 1c - 12th January 2019

Sentinel 2A Level 1C - 8th January 2016

Here is the shapefile that was used to clip the above satellite imagery, as you will see later on in the tutorial.

aoi shapefile - https://www.dropbox.com/sh/0ooii4f2tnojggr/AAAfNCalu8kpN8IKc1QJ9LLSa?dl=0

When analysing satellite images for two or more different time periods, make as much effort to use images from the same season, so that the conditions are similar, as sometimes even the angle of the sun and prevailing seasons affect reflectance values.

Extract both imagery data. We will begin analysing the 2019 image.

#load the packages that assist in raster analysis
library(raster)
## Warning: package 'raster' was built under R version 4.1.3
## Loading required package: sp
library(sp)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/User/Documents/R/win-library/4.1/rgdal/proj
library(rgeos)
## rgeos version: 0.5-9, (SVN revision 684)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-6 
##  Polygon checking: TRUE

Sentinel-2 tile data consists 13 bands, that is, bands 1 - 12 and the last is TCI. We will call each individual band into R. Instead of copying the text address from our directory for all the bands, we will only copy for the first, and copy the calling code for band1 for all the other bands. Of course, the last suffix before the file format B.. will be changed to match the band number of interest.

#call individual bands into R

band1 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B01.jp2")
band2 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B02.jp2")
band3 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B03.jp2")
band4 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B04.jp2")
band5 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B05.jp2")
band6 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B06.jp2")
band7 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B07.jp2")
band8 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B08.jp2")
band9 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B09.jp2")
band10 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B10.jp2")
band11 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B11.jp2")
band12 <- raster("./sentinel-2-data/L1C_T36MZE_A018579_20190112T075614/S2A_MSIL1C_20190112T075301_N0207_R135_T36MZE_20190112T092302.SAFE/GRANULE/L1C_T36MZE_A018579_20190112T075614/IMG_DATA/T36MZE_20190112T075301_B12.jp2")

The bands should now be loaded into your R environment with their respective names.

To view our image tile, we will load just three bands. Only three are needed to visualize an image. Our area of interest is around Lake Nakuru. Worth a visit.

#to create true color composite, the band combination is 4, 3, 2
trucol_2019 <- stack(band4, band3, band2)
trucol_2019
## class      : RasterStack 
## dimensions : 10980, 10980, 120560400, 3  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 799980, 909780, 9890200, 1e+07  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## names      : T36MZE_20190112T075301_B04, T36MZE_20190112T075301_B03, T36MZE_20190112T075301_B02 
## min values :                          0,                          0,                          0 
## max values :                      65535,                      65535,                      65535
#save trucol_2019 to a file directory
writeRaster(trucol_2019, filename = "truecolor_2019.tif", format = "GTiff", overwrite = T)

After stacking the band combination for our image, we will view a true color image of our tile in R.

#plot true color image
plotRGB(trucol_2019, stretch = 'lin', axes = T, main = "Sentinel-2 2019")

Now, that we have seen the full extent of our tile, its time we cropped it to a particular region of interest, around Lake Nakuru. For the purpose of this tutorial, we had created an ROI using the function drawExtent(). Type this function on the help page to see what it does and follow the instructions to create your own ROI.

#save our ROI values into our environment after using drawExtent()
roi <- roi
roi
## class      : Extent 
## xmin       : 835473.8 
## xmax       : 866160 
## ymin       : 9934047 
## ymax       : 9965870
#view our ROI before cropping
plotRGB(trucol_2019, stretch = 'lin', axes = T, ext = c(835473.8, 866160, 9934047, 9965870),  main = "Sentinel-2 2019") #the coordinates of our extent have been put into ext argument

#crop image to our extent
nakuru2019 <- crop(x = trucol_2019, y = c(835473.8, 866160, 9934047, 9965870))
#draw cropped image
plotRGB(nakuru2019, stretch = 'lin', axes = T, main = "Nakuru 2019")

Let’s save our newly cropped image of Nakuru into our directory.

#save cropped image to our directory
writeRaster(nakuru2019, filename = "nakuru2019.tif", format = "GTiff", overwrite = T)

You may have noticed that our satellite image appears hazy, and is far from the clear images we see in other places, say from Google Earth. This is because it has not yet been atmospherically corrected. This is quite a huge topic, but there are lots of sites on the web that summarize it in a very simple way. Atmospheric correction of Sentinel 2 images can be done with the open source software Qgis. For a step by step procedure of conducting atmospheric correction, see pages 63 - 76 of my tutorial. If you need a quick job, see this video - https://www.youtube.com/watch?v=5_T85oJ_gWc

The atmospherically corrected Nakuru 2019 image, and the shapefile used for clipping to our region of interest with the lake at the center are both downloadable from here:

nak2019_atmos_clip - https://www.dropbox.com/sh/tryfml2j6kpo6cj/AADBaRUwRHd773Pae0srSxsWa?dl=0

aoi shapefile - https://www.dropbox.com/sh/0ooii4f2tnojggr/AAAfNCalu8kpN8IKc1QJ9LLSa?dl=0

We used a shapefile just in case the next satellite image may be in a different CRS, which makes the use of uniform coordinates lead to extraction of misplaced ROIs. Why not use the coordinates from drawExtent()? These were used, with the evidence shown below, however the resulting images was a mind-blogging several kilometres below the image shown above.

previously used coordinates

Things didn’t work out as we had expected. Anyway, the shapefile made things get back into shape!

Let’s load the atmospherically corrected image of the same area as that above and see if all this preprocessing clears the mist.

#load the requisite raster packages
library(raster)
library(rgdal)
library(rgeos)
library(sp)
library(sf)
#load the atmospherically corrected raster 
nak2019_atmos <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B0stack_raster.tif")
#call all the bands for the atmospherically corrected image into R
atmb1 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B01.tif")
atmb2 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B02.tif")
atmb3 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B03.tif")
atmb4 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B04.tif")
atmb5 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B05.tif")
atmb6 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B06.tif")
atmb7 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B07.tif")
atmb8 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B08.tif")
atmb8A <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B8A.tif")
atmb9 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B09.tif")
atmb10 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B10.tif")
atmb11 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B11.tif")
atmb12 <- raster(".//sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B12.tif")
#stack the required bands for true color composite image - 4, 3, 2
atmos19_trucol <- stack(atmb4, atmb3, atmb2)
#plot the true color composite image (4, 3, 2)
plotRGB(atmos19_trucol, stretch = 'lin', axes = T, main = "Atmospherically corrected 2019 image")

Viewer, viewer off the screen, which of the two is the clearest of them all? You bet this one.

Part II We have already seen how an unprocessed images looks like. That it is fuzzy and hazzy. We have to be crazy to deal with it for serious analysis at first sight. Therefore, in remote sensing, unless there is a good reason, preprocessing, which includes atmospheric or geometric processing (or both) is a necessary component of satellite image analysis.

Therefore, without further much ado, we will load into R a preprocessed satellite image of the same area but for 2016.

nak2016_atmos_clip - https://www.dropbox.com/sh/j4orvqsnhx56bhz/AACSFnqZ84s4bUnF2R8trbvea?dl=0

#call the raster and requisite libraries
library(raster)
library(rgdal)
library(rgeos)
library(sp)
library(sf)
#load the bands for 2016 Sentinel-2 image
b1_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B01.tif")
b2_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B02.tif")
b3_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B03.tif")
b4_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B04.tif")
b5_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B05.tif")
b6_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B06.tif")
b7_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B07.tif")
b8_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B08.tif")
b8a_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B8A.tif")
b9_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B09.tif")
b10_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B10.tif")
b11_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B11.tif")
b12_2016 <- raster("E:/Documents/R Studio is here/R Studio files projects/gis_coding/sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B12.tif")
#stack the bands for truecolor (4, 3, 2)
trucol_2016 <- stack(b4_2016, b3_2016, b2_2016)
trucol_2016
## class      : RasterStack 
## dimensions : 2223, 2526, 5615298, 3  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 831490, 856750, 9948280, 9970510  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## names      : clip_RT_T36MZE_20160108T075312_B04, clip_RT_T36MZE_20160108T075312_B03, clip_RT_T36MZE_20160108T075312_B02 
## min values :                             0.0149,                                  ?,                                  ? 
## max values :                              0.464,                                  ?,                                  ?

Good thing is that the clipped preprocessed images for 2016 (above) and for 2019 (below) are all in the same datum.

atmos19_trucol
## class      : RasterStack 
## dimensions : 2223, 2526, 5615298, 3  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 831490, 856750, 9948280, 9970510  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## names      : clip_RT_T36MZE_20190112T075301_B04, clip_RT_T36MZE_20190112T075301_B03, clip_RT_T36MZE_20190112T075301_B02 
## min values :                             0.0174,                                  ?,                                  ? 
## max values :                             0.3951,                                  ?,                                  ?
#plot truecolor image of Sentinel-2 2016 image
plotRGB(trucol_2016, stretch = 'lin', axes = T, main = "Atmospherically corrected 2016 image")

#save the 2016 truecolor image as a raster file in our directory
writeRaster(trucol_2016, filename = "trucolor_2016.tif", format = "GTiff", overwrite = T)

Take a look at the 2019 and 2016 folders you downloaded from the dropbox links. You will notice that there is a .tif file ending with B0stack_raster.tif. This raster file is a stack of all the 12 bands of the initial satellite tile. We will try to see if we can plot a RGB straight from this composite rasterstack without calling the individual bands into R. Talk of killing all birds with a single shot.

#call the clipped and preprocessed 12-band composite into R
stack2019 <- brick("./sentinel-2-data/nak2019_atmos_clip/clip_RT_T36MZE_20190112T075301_B0stack_raster.tif")
stack2019
## class      : RasterBrick 
## dimensions : 2223, 2526, 5615298, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 831490, 856750, 9948280, 9970510  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : clip_RT_T36MZE_20190112T075301_B0stack_raster.tif 
## names      : clip_RT_T36MZE_20190112T075301_B0stack_raster.1, clip_RT_T36MZE_20190112T075301_B0stack_raster.2, clip_RT_T36MZE_20190112T075301_B0stack_raster.3, clip_RT_T36MZE_20190112T075301_B0stack_raster.4, clip_RT_T36MZE_20190112T075301_B0stack_raster.5, clip_RT_T36MZE_20190112T075301_B0stack_raster.6, clip_RT_T36MZE_20190112T075301_B0stack_raster.7, clip_RT_T36MZE_20190112T075301_B0stack_raster.8, clip_RT_T36MZE_20190112T075301_B0stack_raster.9, clip_RT_T36MZE_20190112T075301_B0stack_raster.10, clip_RT_T36MZE_20190112T075301_B0stack_raster.11, clip_RT_T36MZE_20190112T075301_B0stack_raster.12, clip_RT_T36MZE_20190112T075301_B0stack_raster.13

The band names are darn too long. Let’s shorten them. But I’m curious, is B0stack_raster13 band 8A or Band 12? We are not sure about this.

#give the bands for 2019 image shorter names
names(stack2019) <- c('B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B10', 'B11', 'B12')
#view the names of the 2019 image
names(stack2019)
##  [1] "B1"  "B2"  "B3"  "B4"  "B5"  "B6"  "B7"  "B8"  "B8A" "B9"  "B10" "B11"
## [13] "B12"

Let’s do the same for the clipped and preprocessed 2016 image.

#call the clipped and preprocessed 12-band 2016 image composite
stack2016 <- brick("./sentinel-2-data/nak2016_atmos_clip/clip_RT_T36MZE_20160108T075312_B0stack_raster.tif")
stack2016
## class      : RasterBrick 
## dimensions : 2223, 2526, 5615298, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 831490, 856750, 9948280, 9970510  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : clip_RT_T36MZE_20160108T075312_B0stack_raster.tif 
## names      : clip_RT_T36MZE_20160108T075312_B0stack_raster.1, clip_RT_T36MZE_20160108T075312_B0stack_raster.2, clip_RT_T36MZE_20160108T075312_B0stack_raster.3, clip_RT_T36MZE_20160108T075312_B0stack_raster.4, clip_RT_T36MZE_20160108T075312_B0stack_raster.5, clip_RT_T36MZE_20160108T075312_B0stack_raster.6, clip_RT_T36MZE_20160108T075312_B0stack_raster.7, clip_RT_T36MZE_20160108T075312_B0stack_raster.8, clip_RT_T36MZE_20160108T075312_B0stack_raster.9, clip_RT_T36MZE_20160108T075312_B0stack_raster.10, clip_RT_T36MZE_20160108T075312_B0stack_raster.11, clip_RT_T36MZE_20160108T075312_B0stack_raster.12, clip_RT_T36MZE_20160108T075312_B0stack_raster.13
#give shorter names to the bands of the stack2016 image
names(stack2016) <- c('b1_2016', 'b2_2016', 'b3_2016', 'b4_2016', 'b5_2016', 'b6_2016', 'b7_2016', 'b8_2016', 'b8a_2016', 'b9_2016', 'b10_2016', 'b11_2016', 'b12_2016')
#view new band names of stack2016
names(stack2016)
##  [1] "b1_2016"  "b2_2016"  "b3_2016"  "b4_2016"  "b5_2016"  "b6_2016" 
##  [7] "b7_2016"  "b8_2016"  "b8a_2016" "b9_2016"  "b10_2016" "b11_2016"
## [13] "b12_2016"
#visualize stack 2019 as true color
plotRGB(stack2019, r = 'B4', g = 'B3', b = 'B2', stretch = 'lin', axes = T, main = "Nakuru 2019 image")

#visualize stack 2016 as true color
plotRGB(stack2016, r = 'b4_2016', g = 'b3_2016', b = 'b2_2016', stretch = 'lin', axes = T, main = "Nakuru 2016 image")

Calling bands from a raster stack is way simpler than calling them individually. However, whenever there are bands with same number with a slight deviation, such as B8 and B8A, you’d better call the bands individually before their use.

Confirmed. In the rasterstacks we have dealt with above, when the preprocessing tools came across a band with the same digit as the previous, but with an additional suffix such as the A in band 8A the tool will simply assign this band with the next numerical as the suffix ie. band8A will become band 9, and band10 will become band11 until we end with band 13.

In short, Band8A is assigned as _B0stack_raster.9 in our rasterstacks.