Chapter 2 Data exploration
2.1 Get remote sensing data
2.2 Data visualization
## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/saif/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/saif/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')
# Single band and composite maps
# You can plot individual layers of a RasterStack 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))
To make a true (natural) color image that looks like a normal photograph (vegetation in green, water in blue…). The true-color reveals much about landscape than gray images.
landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
Anoher popular image visualization method in remote sensing is known as “False color” image. This representation makes it easy to see the vegetaion (in red).
par(mfrow = c(1,2))
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
landsatFCC = stack(b5,b4,b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
2.3 Subset and rename bands
2.3.1 load all the layers
We can create a RasterStack withh the 11 layers. These layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panachromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
## [1] "data/rs/LC08_044034_20170614_B1.tif"
## [2] "data/rs/LC08_044034_20170614_B2.tif"
## [3] "data/rs/LC08_044034_20170614_B3.tif"
## [4] "data/rs/LC08_044034_20170614_B4.tif"
## [5] "data/rs/LC08_044034_20170614_B5.tif"
## [6] "data/rs/LC08_044034_20170614_B6.tif"
## [7] "data/rs/LC08_044034_20170614_B7.tif"
## [8] "data/rs/LC08_044034_20170614_B8.tif"
## [9] "data/rs/LC08_044034_20170614_B9.tif"
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"## class : RasterStack
## dimensions : 1245, 1497, 1863765, 11 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values : 9.641791e-02, 7.483990e-02, 4.259216e-02, 2.084067e-02, 8.457669e-04, -7.872183e-03, -5.052945e-03, 3.931751e-02, -4.337332e-04, 2.897978e+02, 2.885000e+02
## max values : 0.73462820, 0.71775615, 0.69246972, 0.78617686, 1.01243150, 1.04320455, 1.11793602, 0.82673049, 0.03547901, 322.43139648, 317.99530029We can select specific layers (bands)
## [1] 11## [1] 7For clarity, we can set the names of the bands
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"2.3.2 Spatial subset
We can use spatial subsetting in order to limit the analysis to a geographical region. We can specify the extent coordinates similar to the example above and we can use interactive selection from the image with “drawExtent” and “drawPoly” functions.
## class : Extent
## xmin : 594090
## xmax : 639000
## ymin : 4190190
## ymax : 4227540e <- extent(610000, 4200000, 630000,4220000)
# crop landsat by the extent
landsatcrop <- crop(landsat, e)
# plot
par(mfrow = c(1,2))
plot(landsat[[5]])
rect(e[1],e[2],e[3],e[4], border='red', lwd=2)
plot(landsatcrop[[5]])
The subbset image can be saved with “writeRaster” function
2.4 Spectral profiles
The spectral profile is a plot of the spectrum for pixels representing a certain earth surface features. It demonstrates the difference in spectral properties of various earth surface features.
To do that, we need to join information about the lad use and land cover with pixel values of the raster data.
# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')
# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
# add the land cover class to the points
ptsamp$class <- over(ptsamp, samp)$class
# extract values with points
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)## ultra.blue blue green red NIR SWIR1 SWIR2
## [1,] 0.1348463 0.1171285 0.10090701 0.10077689 0.1620628 0.2171030 0.1806697
## [2,] 0.1325042 0.1150899 0.09971425 0.09843475 0.1831637 0.2129608 0.1747276
## [3,] 0.1334800 0.1288825 0.13744865 0.18316367 0.3254918 0.3292869 0.1989731
## [4,] 0.1443883 0.1468389 0.16507718 0.22690523 0.3690165 0.3911583 0.2422809
## [5,] 0.1321572 0.1143743 0.09774078 0.09529021 0.1667254 0.2016405 0.1645567
## [6,] 0.1372752 0.1200344 0.10405154 0.10411660 0.1676145 0.2153030 0.1779589# plot
plot(landsat[[5]])
plot(samp, add = TRUE, border = 'red', lwd = 2)
plot(ptsamp, pch = 21, bg = "black", col="black", cex = 0.45, add = TRUE)
Now, we can compute he 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## ultra.blue blue green red NIR SWIR1
## built 0.1770095 0.16882171 0.16839923 0.18151631 0.23660067 0.23342241
## cropland 0.1122196 0.08998683 0.08391113 0.05331576 0.46224299 0.15512451
## fallow 0.1325273 0.11680835 0.10406427 0.11311789 0.17682463 0.23101104
## open 0.1392278 0.13822205 0.15368719 0.20788541 0.34339042 0.35489589
## water 0.1340949 0.11714744 0.09997358 0.07941412 0.04923834 0.03382958
## SWIR2
## built 0.19172661
## cropland 0.07002481
## fallow 0.19389564
## open 0.21224384
## water 0.02742395Now we plot the mean spectra of these features.
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), 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 differences in the reflectance of different eatures of the earth’s surface. Water shows relatively low reflection in all wavelengths and built surfaces have relatively high reflectance in the longer wavelength.