Exploration of Landsat satellite imagery with R

1. Introduction

In this notebook I illustrate how to find, download and explore satellite remote sensing data with R. I also show how to create color composites and explore spectral profiles.

A lot of the code that follows is based on the notebook written by Aniruddha Ghosh and Robert J. Hijmans which you may find here. However, I expand here both code and explanations to help my students get started on remote sensing image analysis. I use here the recently released terra and luna packages. Thanks a lot, Robert and Ani, for such a valuable work!!!

I will primarily use a spatial subset of a Landsat 8 scene collected in 2013. The subset covers the area known as Montes de Maria, in the Colombian Caribbean region. I chose 2013 Landsat imagery because the most recent official land cover classification in Colombia corresponds to the 2010-2012 time period. By this, I will be able to visually compare land cover data and satellite imagery.

2. Finding and downloading Landsat Image

Landsat images are organized according to the Worldwide Reference System (WRS). There is a Landsat scene that covers every given location on Earth’s landmasses between ~82° N and S. Landsat scenes cover a region of approximately 182 km x 185 km (113 x 115 miles). Each Landsat scene is identified by a Path and Row number.

So if you want to find out which Landsat scene(s) cover a given location you can use the Landsat Acquisition Tool. I did it, and found that the Landsat image scenes covering Montes de Maria have the WRS Path 9 and the WRS Row 53). This information may be useful in different scenarios.

Landsat Acquisition Tool The luna R library provides a function to find and download Landsat data for a specific product, area, and time period. However, I tried it several times with no success.

For such a reason, in this notebook, I will use the *U.S. Global Visualization Viewer (GloVis) to find a Landsat 8 imagery over Montes de Maria:

Global Visualization Viewer

I downloaded the Landsat 8 image with ID LC08_L1TP_009053_20131223_20170427_01_T1. Based on this Landsat naming convention, you can see that the Sensor-Satellite is OLI/TIRS combined Landsat 8, WRS Path 009, WRS Row 053 and it was acquired on 23 December 2013.

You can also download a Landsat scene which covers the region you are interested in. I would suggest to select a region you already know. Furthermore, I would advice to select a time period for which land cover information is available.

After downloading your Landsat scene, take your time to answer the following questions:

What is the processing date of this scene?
What is the processing correction level of this scene?
What does mean that level?

Note that the downloaded Landsat scene was delivered as a compressed file: LC08_L1TP_009053_20131223_20170427_01_T1.tar.gz. The double extension is because it is a tar file compressed with gzip.

Once you decompress the file, note that it contains separate files for each band: Landsat individual band files

Do you want to understand what is the data included in the quality assessment band (denoted as BQA)? See this link

Note that there is a file with the prefix MTL.txt. It contains the metadata. See the following figure which shows a subset of the Landsat scene metadata: Landsat metadata

3. Exploring and visualizing the data

3.1 Reading individual bands

Let’s use the terra library to create SpatRaster objects for single Landsat layers (bands)

library(terra)

## terra version 1.0.10

# Blue
b2 <- rast('./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B2.TIF')
# Green
b3 <- rast('./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF')
# Red
b4 <- rast('./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B4.TIF')
# Near Infrared (NIR)
b5 <- rast('./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B5.tif')

Let’s check what is every variable:

b5

## class       : SpatRaster 
## dimensions  : 7741, 7581, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : LC08_L1TP_009053_20131223_20170427_01_T1_B5.TIF 
## name        : LC08_L1TP_009053_20131223_20170427_01_T1_B5

You can see the spatial resolution, extent, number of layers, coordinate reference system and more.

3.2 Image information and statistics

The below shows how you can access various properties from a Raster* object (this is the same for any raster data set).

# coordinate reference system (CRS)
crs(b5)

## [1] "PROJCRS[\"WGS 84 / UTM zone 18N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World - N hemisphere - 78°W to 72°W - by country\"],\n        BBOX[0,-78,84,-72]],\n    ID[\"EPSG\",32618]]"

# Number of cells, rows, columns
ncell(b5)

## [1] 58684521

dim(b5)

## [1] 7741 7581    1

# spatial resolution
res(b5)

## [1] 30 30

# Number of bands
nlyr(b5)

## [1] 1

## histogram
hist(b5, maxcell=5000000, plot=TRUE)

## mean
(b5mean <- global(b5, fun="mean"))

## range
(b5range <- global(b5, fun="range"))

Note that the digital values at each band represent a scaled spectral radiance. If we want to convert them to aparent reflectance, we need to review the metadata:

Reflectance conversion factors

## standard deviation
(b5sd <- global(b5, fun="sd"))

# Do the bands have the same extent, number of rows and columns, projection, 
# resolution, and origin
compareGeom(b2,b3)

## [1] TRUE

You can create a SpatRaster with multiple layers from the existing SpatRaster (single layer) objects.

s <- c(b5, b4, b3)

# Check the properties of the multi-band image
s

## class       : SpatRaster 
## dimensions  : 7741, 7581, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## sources     : LC08_L1TP_009053_20131223_20170427_01_T1_B5.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B4.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF  
## names       : LC08_L1TP_~7_01_T1_B5, LC08_L1TP_~7_01_T1_B4, LC08_L1TP_~7_01_T1_B3

You can also create the multi-band image using the filenames.

# first create a list of raster layers to use
filenames <- paste0('./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B', 1:7, ".TIF")

filenames

## [1] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B1.TIF"
## [2] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B2.TIF"
## [3] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF"
## [4] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B4.TIF"
## [5] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B5.TIF"
## [6] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B6.TIF"
## [7] "./landsat/LC08_09_53_2013/LC08_L1TP_009053_20131223_20170427_01_T1_B7.TIF"

(landsat <- rast(filenames))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## sources     : LC08_L1TP_009053_20131223_20170427_01_T1_B1.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B2.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF  
##               ... and 4 more source(s)
## names       : LC08_~T1_B1, LC08_~T1_B2, LC08_~T1_B3, LC08_~T1_B4, LC08_~T1_B5, LC08_~T1_B6, ...

Above we created a SpatRaster with 7 layers. The layers contains digital values associated to reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, and Shortwave Infrared (SWIR) 2. We leave out four bands: Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS).

3.2 Single band and composite visualization

You can plot individual layers of a multi-spectral image.

par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))

Have a look at the legends of the maps created above. They can range between 0 and 1. Notice the difference in shading and range of legends between the different bands. This is because different surface features reflect the incident solar radiation differently. Each layer represent how much incident solar radiation is reflected for a particular wavelength range. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.

Note also that it is possible to improve image contrast using the stretch parameter. It can take one of the following values: “lin” (for linear stretch) or “hist” (for histogram equalization):

We do not gain that much information from these grey-scale plots; they are often combined into a “composite” to create more interesting plots. You can learn more about color composites in remote sensing here and also in the section below.

To make a “true (or natural) color” image, that is, something that looks like a normal photograph (vegetation in green, water blue etc), we need bands in the red, green and blue regions. For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used. The plotRGB method can be used to combine them into a single composite. You can also supply additional arguments to plotRGB to improve the visualization (e.g. a linear stretch of the values, using stretch = “lin”).

Before anything, let’s convert the raw digital numbers into apparent reflectance (in percent):

(nb4 <- round(b4*2E-03))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        : LC08_L1TP_009053_20131223_20170427_01_T1_B4 
## min value   :                                           0 
## max value   :                                          84

(nb3 <- round(b3*2E-03))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        : LC08_L1TP_009053_20131223_20170427_01_T1_B3 
## min value   :                                           0 
## max value   :                                          77

(nb2 <- round(b2*2E-03))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        : LC08_L1TP_009053_20131223_20170427_01_T1_B2 
## min value   :                                           0 
## max value   :                                          74

(landsatRGB <- c(nb4, nb3, nb2))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## sources     : memory  
##               memory  
##               memory  
## names       : LC08_L1TP_~7_01_T1_B4, LC08_L1TP_~7_01_T1_B3, LC08_L1TP_~7_01_T1_B2 
## min values  :                     0,                     0,                     0 
## max values  :                    84,                    77,                    74

# for better visualization we apply linear contrast stretch
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

The true-color composite reveals much more about the landscape than the earlier gray images. Another popular image visualization method in remote sensing is known “false color” image in which NIR, red, and green bands are combined. This representation is popular as it makes it easy to see the vegetation (in red).

landsatFCC <- c(b5, b4, b3)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

Note: Always check for terra package documentation (help(plotRGB)) for other arguments that can be added (like scale) to improve or modify the image.

Use the plotRGB function with multi-band (7 bands) landsat image to create several color composites.

3.4 Subset and rename bands

You can select specific layers (bands) using subset function, or via indexing.

# select first 3 bands only
(landsatsub1 <- subset(landsat, 1:3))

## class       : SpatRaster 
## dimensions  : 7741, 7581, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## sources     : LC08_L1TP_009053_20131223_20170427_01_T1_B1.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B2.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF  
## names       : LC08_L1TP_~7_01_T1_B1, LC08_L1TP_~7_01_T1_B2, LC08_L1TP_~7_01_T1_B3

# same
(landsatsub2 <- landsat[[1:3]])

## class       : SpatRaster 
## dimensions  : 7741, 7581, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 409485, 636915, 1002285, 1234515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## sources     : LC08_L1TP_009053_20131223_20170427_01_T1_B1.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B2.TIF  
##               LC08_L1TP_009053_20131223_20170427_01_T1_B3.TIF  
## names       : LC08_L1TP_~7_01_T1_B1, LC08_L1TP_~7_01_T1_B2, LC08_L1TP_~7_01_T1_B3

# Number of bands in the original and new data
nlyr(landsat)

## [1] 7

nlyr(landsatsub1)

## [1] 3

nlyr(landsatsub2)

## [1] 3

If you will not use all bands in landsat, you can remove those using:

# uncommnent if needed
#landsat <- subset(landsat, 1:5)

For clarity, it is useful to set the names of the bands. (source)[https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites?qt-news_science_products=0#qt-news_science_products]

names(landsat)

## [1] "LC08_L1TP_009053_20131223_20170427_01_T1_B1"
## [2] "LC08_L1TP_009053_20131223_20170427_01_T1_B2"
## [3] "LC08_L1TP_009053_20131223_20170427_01_T1_B3"
## [4] "LC08_L1TP_009053_20131223_20170427_01_T1_B4"
## [5] "LC08_L1TP_009053_20131223_20170427_01_T1_B5"
## [6] "LC08_L1TP_009053_20131223_20170427_01_T1_B6"
## [7] "LC08_L1TP_009053_20131223_20170427_01_T1_B7"

names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

names(landsat)

## [1] "ultra-blue" "blue"       "green"      "red"        "NIR"       
## [6] "SWIR1"      "SWIR2"

3.5 Spatial subset or crop

Spatial subsetting can be used to limit analysis to a geographic subset of the image. Spatial subsets can be created with the crop function, using an extent object, or another spatial object from which an Extent can be extracted.

# Using extent
ext(landsat)

## SpatExtent : 409485, 636915, 1002285, 1234515 (xmin, xmax, ymin, ymax)

# enter the values you need
#e <- ext(500000, 600000, 1150000, 1200000)
# crop landsat by the extent
#landsatcrop <- crop(landsat, e)

I will use an spatial object to subset the Landsat scene:

# the spatial object must have the same CRS as the Landsat image
(montes <- vect("./shapes/MontesdeMaria_UTM18N.shp"))

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 15, 11  (geometries, attributes)
##  extent      : 422753.2, 525002.6, 1019785, 1121539  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       : DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR     MPIO_CRSLC MPIO_NAREA
##  type        :      <chr>      <chr>           <chr>          <chr>      <num>
##  values      :         13      13212         CÓRDOBA           1909      597.3
##                        13      13244 EL CARMEN DE B~ Ley 13 de 1857      946.3
##                        13      13248        EL GUAMO           1853      383.2
##  MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area      layer            path
##      <num>      <chr>      <num>      <num>      <chr>           <chr>
##       2017    BOLÍVAR      1.206    0.04918 bol_montes /Users/ivanliz~
##       2017    BOLÍVAR      2.069    0.07793 bol_montes /Users/ivanliz~
##       2017    BOLÍVAR     0.8978    0.03159 bol_montes /Users/ivanliz~

(e <- ext(montes))

## SpatExtent : 422753.23555044, 525002.59027505, 1019785.4193537, 1121539.31032031 (xmin, xmax, ymin, ymax)

# crop landsat by the extent
(landsatcrop <- crop(landsat, e))

## class       : SpatRaster 
## dimensions  : 3392, 3409, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 422745, 525015, 1019775, 1121535  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## names       : ultra-blue,  blue, green,   red,   NIR, SWIR1, ... 
## min values  :          0,     0,     0,     0,     0,     0, ... 
## max values  :      20788, 21783, 23294, 27555, 43270, 65535, ...

plotRGB(landsatcrop, r=5, g=6, b=4, axes=TRUE, stretch="lin", main="Landsat subset RGB564 color composite")

3.6 Saving results to disk

At this stage we may want to save the raster to disk using the function writeRaster. Multiple file types are supported. We will use the commonly used GeoTiff format. While the layer order is preserved, layer names are unfortunately lost in the GeoTiff format.

writeRaster(landsatcrop, filename="./landsat/montes-landsat.TIF", 
            datatype="INT1U", overwrite=TRUE)

Note: Check for package documentation (help(writeRaster)) for additional helpful arguments that can be added.

What are the possible datatype values?
What does mean every possible datatype values? See here

3.7 Relation between bands

A scatterplot matrix can be helpful in exploring relationships between raster layers. This can be done with the pairs() function of the terra package.

Let’s make a plot of digital values in the blue wavelength against digital values in the green wavelength.

pairs(landsatcrop[[2:3]], main = "Blue versus Green")

The first plot reveals high correlations between the blue and the green wavelength regions. Because of the high correlation, we could use one of the two bands without losing much information.

Now, let’s make a plot of digital values in the red wavelength against digital values in the near infrared wavelength.

pairs(landsatcrop[[4:5]], main = "Red versus NIR")

This distribution of points in second plot (between NIR and red) is unique due to its triangular shape. Vegetation reflects very highly in the NIR range than red and creates the upper corner close to NIR (y) axis. Water absorbs energy from all the bands and occupies the location close to origin. The furthest corner is created due to highly reflecting surface features like bright soil or concrete. Read more about soil lines (here)[http://www.ipgp.fr/~jacquemoud/publications/baret1993a.pdf].

3.8 Extract pixel values

Often we want to get the values of raster cells for specific geographic locations or area. The extract function is used to get raster values at the locations of other spatial data. You can use points, lines, polygons or an Extent (rectangle) object. You can also use cell numbers to extract values. When using points, extract returns the values of a Raster* object for the cells in which a set of points fall.

I will use a 2012 landcover data set provided by the IDEAM, the Colombian environmental studies agency:

## (landcover <- vect("./shapes/MontesLC_UTM18N.shp"))
library(raster)

## Loading required package: sp

(landcover <- shapefile("./shapes/MontesLC_UTM18N.shp"))

## class       : SpatialPolygonsDataFrame 
## features    : 3426 
## extent      : 422788.8, 524926.8, 1019785, 1121539  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 4
## names       : OBJECTID, CODIGO,                                LEYENDA3N,  LC 
## min values  :   157436,    111,            1.1.1. Tejido urbano continuo, 111 
## max values  :   174297,    523, 5.2.3. Estanques para acuicultura marina, 523

## S4 method for signature 'SpatExtent'
#ptsamp <- spatSample(landcover, 100, method="random", strata="LC")
library(sp)
ptsamp <- spsample(landcover, 100, type = "random")

## Warning in proj4string(obj): CRS object has comment, which is lost in output

# add the land cover class to the points
ptsamp$class <- over(ptsamp, landcover)$LC

ptsamp <- vect(ptsamp)
# We use the x-y coordinates to extract the spectral values for the locations
xy <- as.matrix(geom(ptsamp)[,c('x','y')])
df <- extract(landsat, xy)
# To see some of the digital values
head(df)

3.9 Spectral profiles

A plot of the spectrum (all bands) for pixels representing a certain earth surface features (e.g. water) is known as a spectral profile. Such profiles demonstrate the differences in spectral properties of various earth surface features and constitute the basis for image analysis. Spectral values can be extracted from any multispectral data set using extract function. In the above example, we extracted values of Landsat data for the samples. These samples include: forests, shrubs, grassland, cropland, water, fallow, built and open. First we compute the mean reflectance values for each class and each band.

ms <- aggregate(df, list(ptsamp$class), mean)
# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms <- ms * 2E-05 
ms

Now we plot the mean spectra of these features.

# Create a vector of color for the land cover classes for use in plotting
mycolor <- c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3" ,"#2166AC", "#9ea832", "#32a873", "#3292a8", "#5532a8", "#9e32a8",
             "#a8325d", "#0ce2ed", "#32a8a0", "#3267a8", "#324ea8")
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.8), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

The spectral profile shows (dis)similarity in the reflectance of different features on the earth’s surface (or above it). ‘Water’ shows relatively low reflection in all wavelengths, and ‘built’, ‘fallow’ and ‘open’ have relatively high reflectance in the longer wavelengts.

NOTE: This notebook illustrated terra functionalities for remote sensing analysis. It should not be taken as an example of a technical report.

4. Reproducibility

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] raster_3.4-5 sp_1.4-5     terra_1.0-10
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        codetools_0.2-18  lattice_0.20-41   digest_0.6.27    
##  [5] grid_4.0.4        jsonlite_1.7.2    magrittr_2.0.1    evaluate_0.14    
##  [9] highr_0.8         rlang_0.4.10      stringi_1.5.3     vctrs_0.3.6      
## [13] rmarkdown_2.7     rgdal_1.5-23      tools_4.0.4       stringr_1.4.0    
## [17] xfun_0.21         yaml_2.2.1        compiler_4.0.4    htmltools_0.5.1.1
## [21] knitr_1.31