Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/

Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ



Objetivo de la práctica: Introducir los conceptos básicos de sensores remotos y el análisis de imágenes satelitales, utilizando la librería RSToolBox.Esta práctica se basa en el material del curso “Spatial Data Analysis and Modeling for Agricultural Development” (https://gfc.ucdavis.edu/events/arusha2016/), en particular en el material del día 4.

Slides (https://gfc.ucdavis.edu/events/arusha2016/_static/talks/day4_1_remotesensing.pdf)
Handouts (https://gfc.ucdavis.edu/events/arusha2016/_static/labs/day4/day4_lab1_remote-sensing.pdf)
Datos (https://gfc.ucdavis.edu/events/arusha2016/_static/labs/day4/day4_remotesensing_data.zip).

Una presentación teórica sobre teledetección y clasificación puede accederse aqui [https://docs.google.com/presentation/d/1d5iH5vPFYqNrV0TnI-o8_m1YtAD8XInxls6P9BWWKRM/edit#slide=id.p14]

1 Basic image manipulation and visualization

1.0.0.1 Read image and display basic image properties

library(raster)
## Loading required package: sp
r <- brick("landsat8-2016march.tif") #load Landsat data
s <- brick('sentinel.tif') #load Sentinel data

r #6 bandas, 30m resolution
## class       : RasterBrick 
## dimensions  : 642, 593, 380706, 6  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 251940, 269730, -382140, -362880  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\R_ecohidrologia\remotesensing\landsat8-2016march.tif 
## names       : landsat8.2016march.1, landsat8.2016march.2, landsat8.2016march.3, landsat8.2016march.4, landsat8.2016march.5, landsat8.2016march.6
s #10 bandas, 10m resolution
## class       : RasterBrick 
## dimensions  : 1934, 1786, 3454124, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 10, 10  (x, y)
## extent      : 251900.9, 269760.9, -382182.3, -362842.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\R_ecohidrologia\remotesensing\sentinel.tif 
## names       : sentinel.1, sentinel.2, sentinel.3, sentinel.4, sentinel.5, sentinel.6, sentinel.7, sentinel.8, sentinel.9, sentinel.10 
## min values  :   575.5850,   401.0961,   235.0271,   307.0000,   464.0000,   489.0000,   470.8530,   371.0000,   149.0000,     89.0000 
## max values  :   6397.991,   6278.674,   6372.010,   5482.000,   6586.000,   7242.000,  11990.447,   7095.000,   7738.000,    5612.000

1.0.0.2 Image information and statistics

# Summarize a Raster* object. A sample is used for very large files.
# summary(lsat1, maxsamp = 5000)
# number of bands
nlayers(r)
## [1] 6
nlayers(s)
## [1] 10
# coordinate reference system (CRS)
crs(r)
## CRS arguments:
##  +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
crs(s)
## CRS arguments:
##  +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# get spatial resolution
xres(r)
## [1] 30
yres(r)
## [1] 30
res(r)
## [1] 30 30
res(s)
## [1] 10 10
# Number of rows, columns, or cells
ncell(r)
## [1] 380706
dim(s)
## [1] 1934 1786   10

For Landsat and Sentinel, these are different spectral bands. The term ‘band’ is frequently used in remote sensing to refer to a variable (layer) in a multi-variable dataset as these variables typically represent reflection in different bandwidths in the electromagnetic spectrum.

1.0.0.3 Visualize single and multi-band imagery

To display 3-band color image, we use plotRGB. We have select the index of bands we want to render in the red, green and blue regions. For this Landsat image, r = 3 (red), g = 2(green), b = 1(blue) will plot the true color composite (vegetation in green, water blue etc). Selecting r = 4 (NIR), g = 3 (red), b = 2(green) will plot the false color composite (very popular in remote sensing with vegetation as red).

nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(r, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
plotRGB(r, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")

nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(s, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "Senitnel True Color Composite")
plotRGB(s, r = 7, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Sentinel False Color Composite")

You can see the effect of differences in spatial resolution. Sentinel image has 10 m pixel size compared to 30 m pixels of Landsat.

1.0.0.4 Subset and rename spectral bands

# select first 3 bands only
rsub <- subset(r, 1:3)
# Number of bands in orginal and new data
nlayers(r)
## [1] 6
nlayers(rsub)
## [1] 3
#Set the names of the bands using the following:
# For LANDSAT
names(r)
## [1] "landsat8.2016march.1" "landsat8.2016march.2" "landsat8.2016march.3"
## [4] "landsat8.2016march.4" "landsat8.2016march.5" "landsat8.2016march.6"
names(r) <- c('blue','green','red','NIR','SWIR1','SWIR2')

# For SENTINEL
names(s)
##  [1] "sentinel.1"  "sentinel.2"  "sentinel.3"  "sentinel.4"  "sentinel.5" 
##  [6] "sentinel.6"  "sentinel.7"  "sentinel.8"  "sentinel.9"  "sentinel.10"
names(s) <- c('blue','green','red','rededge1','rededge2','rededge3','NIR','rededge4','SWIR1','SWIR2')

1.0.0.5 Spatial subset/crop

Spatial subsetting can be used to limit analysis to a geographic subset of the image. Spatial subsets can be selected using the following methods: using extent object, using another spatial object from which an Extent can be extracted/created or entering row and column numbers.

# Using extent
extent(r)
## class       : Extent 
## xmin        : 251940 
## xmax        : 269730 
## ymin        : -382140 
## ymax        : -362880
e <- extent(255000, 260000, -375000, -370000)

# crop LANDSAT
rr <- crop(r, e)

# crop SENTINEL
ss <- crop(s, e)

# Now plot the subsets side-by-side to compare the resolution
nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(rr, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat Subset")
plotRGB(ss, r = 7, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Sentinel Subset")

Activity: Explore the functions drawExtent and drawPoly to select area of your interests.

1.0.0.6 Change spatial resolution or pixel size

aggregate creates a lower resolution (larger cells) image from finer pixels. Aggregation groups rectangular areas to create larger cells. Reverse process is known as disaggreagte which creates higher resolution (smaller cells) from larger cells. These methods are particularly useful to compare (or analyze) data from different sources; for example, MODIS data 250 m with Landsat 30 m

ssa <- aggregate(ss, fact=3)

nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(rr, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat 30m")
plotRGB(ssa, r = 7, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Sentinel 30m")

1.0.0.7 Extract raster values

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

# extract values with points
samp <- readRDS('samples.rds')
df <- extract(s,samp)

# To see the reflectance values
head(df)
##          blue    green      red rededge1 rededge2 rededge3      NIR
## [1,] 1531.706 1466.710 1399.800 1428.046 2236.057 2477.101 2626.426
## [2,] 2735.996 2643.311 2737.977 2908.017 3597.897 3957.831 3979.245
## [3,] 2487.741 2424.093 2505.525 2699.732 3583.802 3971.218 3929.050
## [4,] 2262.449 2195.070 2176.647 2424.051 3416.316 3811.576 3733.057
## [5,] 2299.131 2219.972 2194.162 2381.674 3268.081 3622.779 3632.741
## [6,] 1962.716 2012.093 2078.173 2367.000 3553.000 4057.000 3810.069
##      rededge4    SWIR1    SWIR2
## [1,] 2698.097 1831.076 1083.050
## [2,] 4163.326 3667.964 2560.216
## [3,] 4230.438 3371.761 2221.030
## [4,] 4067.465 3254.127 2198.960
## [5,] 3860.959 2988.901 1980.681
## [6,] 4398.000 3109.000 1950.000

1.0.0.8 Plot spectra

Plot of the spectrum (all bands) for the pixel/area is also known as spectral profiles. These profiles demonstrate the differences in spectral propitiates of various Earth surface features and constitute basis for numerous techniques for remote sensing 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: cloud, forest, crop, fallow, built-up, open-soil, water and grassland. In the next example, we will plot the mean spectra of these features.

The spectral profile shows (dis)similarity between different features. Clouds are generally bright and highly reflective in all wavelengths. Crop and forest show similar spectral feature (also see built-up and open areas). Water shows relatively low reflection.

df <- round(df)
# create an empty data.frame to store the mean spectra
ms <- matrix(NA, nrow = length(unique(samp$id)), ncol = nlayers(ss))

for (i in unique(samp$id)){
 x <- df[samp$id==i,]
 ms[i,] <- colMeans(x)
}

# Specify the row- and column names
rownames(ms) <- unique(samp$class)
colnames(ms) <- names(ss)

# We will create a vector of color for the land use land cover classes and will resuse it for other pl
mycolor <- c('cyan', 'darkgreen', 'yellow', 'burlywood', 'darkred', 'darkgray', 'blue', 'lightgreen')

# First create an empty plot
plot(1, ylim=c(300, 4000), xlim = c(1,10), xlab = "Bands", ylab = "Reflectance", xaxt='n')

# Custom X-axis
axis(1, at=1:10, lab=colnames(ms))
# add the other spectra
for (i in 1:nrow(ms)){
 lines(ms[i,], type = "o", lwd = 3, lty = 1, col = mycolor[i])
 }
title(main="Spectral Profile from Sentinel", font.main = 2)# Title
legend("topleft", rownames(ms), cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")# Finally the legend

2 Basic mathematical operations

2.0.0.1 Compute vegetation indices

Normalized Difference Vegetation Index (NDVI). Let’s define a general function for ratio based vegetation index. i and k are the index of bands to be used for the indices computation.

vi <- function(img, i, k){
bi <- img[[i]]
bk <- img[[k]]
vi <- (bk-bi)/(bk+bi)
return(vi) }

# For Sentinel NIR = 7, red = 3.

ndvi <- vi(ss, 3,7)
plot(ndvi, col = rev(terrain.colors(30)), main = 'NDVI from Sentinel')

2.0.0.2 Thresholding

We can apply basic rules to get an estimate of spatial extent of different Earth surface features. Note that NDVI values are standardized and ranges between -1 to +1. Higher values indicate more green cover. Pixels having NDVI values greater than 0.4 are definitely vegetation. Following operation masks all nonvegetation pixels.

veg <- calc(ndvi, function(x){x[x < 0.4] <- NA; return(x)})
plot(veg, main = 'Veg cover')

You can also create classes for different density of vegetation, 0 being no-vegetation.

vegc <- reclassify(veg, c(-Inf,0.2,0, 0.2,0.3,1, 0.3,0.4,2, 0.4,0.5,3, 0.5, Inf, 4))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI based thresholding')

3 Image classification

We will explore two classification methods: unsupervised and supervised. Various unsupervised and supervised classification algorithms may be used. Different choice of classifier may produce different results. We will explore two k-means (unsupervised) and decision tree (supervised) algorithms.

3.0.0.1 Unsupervised classification

In unsupervised classification, we don’t supply any training data. This particularly is useful when we don’t have prior knowledge of the study area. The algorithm groups pixels with similar spectral characteristics into unique clusters/classes/groups following some statistically determined conditions. You have to re-label and combines these spectral clusters into information classes (for e.g. land-use land-cover). We will perform unsupervised classification of the ndvi layer generated previously. To perform unsupervised classification, a good practice is to start with large number of centers (more clusters) and merge/group/recode similar clusters by inspecting the original imagery. Unsupervised algorithms are often referred as clustering.

nr <- getValues(ndvi)
nr.km <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 3, algorithm="Lloyd")
knr <- ndvi
knr[] <- nr.km$cluster
plot(knr, main = 'Unsupervised classification of Sentinel data')

4 Trabajo Calificado 6p

1- Descargue dos imagenes Landsat para su zona de estudio y de dos fechas diferentes y con poca cobertura de nubes. Utilice la página web https://earthexplorer.usgs.gov/ y el producto Landsat Surface Reflectance Level-2 Science Products (Landat collection 1-level2 on demand), ya que presenta valores de reflectancia. 2p

Si llega a requerir (productos con valores digitales y no reflectancia) debe aplicar la correccion radiometrica para convertir los valores digitales a valores de reflectancia, es decir se transforman los valores de numeros digitales a valores de reflectancia mediante el metodo apreciativo y temperatura de brillo. Para ello necesita descargar los archivos de metadatos y de las imagenes y guardarlos en el directorio *C:-library.4.

Guiese por el siguiente codigo. -Import meta-data and bands based on MTL file. Important: copy all landsat images file within *C:-library.4

Ayuda funcion radCor: http://bleutner.github.io/RStoolbox/rstbx-docu/radCor.html.

mtlFile  <- system.file("LC08_L1TP_010062_20161120_20170318_01_T1_MTL.txt", package="RStoolbox", mustWork = T)
metaData <- readMeta(mtlFile)
lsat     <- stackMeta(mtlFile)
lsat_ref <- radCor(lsat, metaData = metaData, method = "apref")
lsat_ref
plot(lsat_ref)

2- Calcule el indice NDVI de las dos imagenes. 2p

3- Cuantifique el aumento y disminucion NDVI de la segunda imagen respecto de la primera. 2p

Has llegado al final de este material!