Practicas 1 a 3 - Percepción Remota

library(knitr)

Contenido

1. Practica 1: Descarga de datos básicos para análisis
2. Practica 2: Exploración de imágenes
   2.1. Cargando imágenes a R
   2.2. Información contenida en las variables
   2.3. Otra alternativa para cargar imágenes
   2.4. Banda única y mapas compuestos
   2.5. Composición de imágenes en color verdadero y falso
   2.6. Subconjuntos de archivos
   2.7. Recortes espaciales de imágenes
   2.8. Composición de color verdadero y falso para imágenes recortadas
   2.9. Almacenamiento de datos en disco
   2.10. Relación entre bandas
   2.11. Valores de pixeles
   2.12. Perfiles espectrales
   
3. Practica 3: Operaciones matemáticas básicas
   3.1. Índice de vegetación NDVI
   3.1.1. Histograma de NDVI
   3.2. Índice de agua de diferencia normalizada NDWI
   3.2.1. Mc Feeters (2006)
   3.2.2. Histograma de NDWI – Mc Feeters (2006)
   3.2.3. M-NDWI Xu (2006)
   3.2.4. Histograma de MNDWI – Xu (2006)
   3.3. Índice de diferencia normalizada de área construida NDBI
   3.3.1. Histograma de NDBI
   3.4. Umbralización
   3.5. Análisis de componentes principales

1. Practica 1: Descarga de datos básicos para análisis

Los datos pueden ser descargados, y guardados directamente en la dirección de preferencia para su posterior análisis, mediante la creación de un directorio (dir.create) y posteriormente descomprimidos (unzip).

dir.create('data', showWarnings = FALSE)
if (!file.exists('data/rs/samples.rds')) {
    download.file('https://biogeo.ucdavis.edu/data/rspatial/rsdata.zip', dest = 'data/rsdata.zip')
    unzip('data/rsdata.zip', exdir='data')
}

2. Practica 2: Exploración de imágenes

2.1. Cargando imágenes a R

Las imágenes previamente descargadas y descomprimidas en la carpeta /data/rs/, se cargan empleando múltiples opciones, mediante las librerías raster y sp, serán subidas las primeras 7 bandas de forma individual, las cuales se encuentran identificadas a continuación.

Banda Landsat 8	Codificación
Azul	b2
Verde	b3
Roja	b4
NIR	b5
SWIR1	b6
SWIR 2	b7

library(sp)
library (raster)

b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')
b6 <- raster('data/rs/LC08_044034_20170614_B6.tif')
b7 <- raster('data/rs/LC08_044034_20170614_B7.tif')

2.2. Información contenida en las variables

Debido a que las bandas han sido almacenadas como objetos (ya que R emplea un sistema de programación enfocado en objetos), se pueden visualizar la información contenida dentro de estas variables recientemente creadas.

Banda 2

b2

class      : RasterLayer 
dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
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 
source     : C:/Users/Elkin/Documents/R/data/rs/LC08_044034_20170614_B2.tif 
names      : LC08_044034_20170614_B2 
values     : 0.0748399, 0.7177562  (min, max)

crs(b2)

CRS arguments:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

ncell(b2)

[1] 1863765

dim(b2)

[1] 1245 1497    1

res(b2)

[1] 30 30

nlayers(b2)

[1] 1

Banda 3

b3

class      : RasterLayer 
dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
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 
source     : C:/Users/Elkin/Documents/R/data/rs/LC08_044034_20170614_B3.tif 
names      : LC08_044034_20170614_B3 
values     : 0.04259216, 0.6924697  (min, max)

crs(b3)

CRS arguments:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

ncell(b3)

[1] 1863765

dim(b3)

[1] 1245 1497    1

res(b3)

[1] 30 30

nlayers(b3)

[1] 1

Banda 4

b4

class      : RasterLayer 
dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
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 
source     : C:/Users/Elkin/Documents/R/data/rs/LC08_044034_20170614_B4.tif 
names      : LC08_044034_20170614_B4 
values     : 0.02084067, 0.7861769  (min, max)

crs(b4)

CRS arguments:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

ncell(b4)

[1] 1863765

dim(b4)

[1] 1245 1497    1

res(b4)

[1] 30 30

nlayers(b4)

[1] 1

Banda 5

b5

class      : RasterLayer 
dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
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 
source     : C:/Users/Elkin/Documents/R/data/rs/LC08_044034_20170614_B5.tif 
names      : LC08_044034_20170614_B5 
values     : 0.0008457669, 1.012432  (min, max)

crs(b5)

CRS arguments:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

ncell(b5)

[1] 1863765

dim(b5)

[1] 1245 1497    1

res(b5)

[1] 30 30

nlayers(b5)

[1] 1

Comparación entre las diferentes bandas, un valor de TRUE permitirá la interacción correcta entre ellas posteriormente.

compareRaster(b2,b3,b4,b5)

[1] TRUE

2.3. Otra alternativa para cargar las imágenes

Al crear un objeto que contenga el listado de nombres que serán cargados, puede emplearse un (stack) donde estén almacenadas las bandas requeridas.

filenames <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- stack(filenames)
landsat

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.99530029 

2.4. Banda única y mapas compuestos

Mediante una variación del código, donde la función (stretch) que da franjas mínimas y máximas en forma de cuantiles y estira en este intervalo de salida los valores originales, se logra una mejor resolución de la imagen.

par(mfrow = c(1,2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(stretch(b2, minq=0.5, maxq=0.9), main = "Blue", col = gray(0:100 / 100))

plot(stretch(b3, minq=0.5, maxq=0.9), main = "Green", col = gray(0:100 / 100))

plot(stretch(b4, minq=0.5, maxq=0.9), main = "Red", col = gray(0:100 / 100))

plot(stretch(b5, minq=0.5, maxq=0.9), main = "NIR", col = gray(0:100 / 100))

2.5. Composición de imágenes en color verdadero Y falso

Debido a que las imágenes proyectadas en blanco y negro de las bandas independientes no arrojan mucha información acerca de la cobertura del suelo, es posible emplear una combinación de bandas que sea procesada por el programa, con el fin de generar una imagen RGB (Red-Green-Blue), lo que permitirá en principio emplear las bandas 4, 3 y 2 respectivamente para elaborar una imagen en color verdadero de la superficie del terreno.

landsatRGB <- landsat[[c(4,3,2)]]
plotRGB(landsatRGB, axes =TRUE, stretch = "lin", main = "Landsat True Color Composite")

Sin embargo, reemplazando las bandas anteriormente mencionadas, por otras diferentes, y generando combinaciones adecuadas, es posible identificar información que en color verdadero no son fácilmente clasificables, es allí donde las imágenes en color falso, donde las bandas no corresponden con el sistema RGB, presentan una gran utilidad.
Debido a que la banda NIR presenta una gran reflectancia de la vegetación, esta es representada por coloraciones rojas intensas.

landsatFCC_veg <- landsat[[c(5,4,3)]]
plotRGB(landsatFCC_veg, axes =TRUE, stretch = "lin", main = "Landsat False Color Composite-Vegetation")

Una opción adicional es clasificar las zonas que tienen reflectancias infrarrojas de onda corta mayores, que por lo general son áreas urbanas, con coberturas de concreto.

landsatFCC_urb <- landsat[[c(7,6,4)]]
plotRGB(landsatFCC_urb, axes =TRUE, stretch = "lin", main = "Landsat False Color Composite-Urban Zones")

De igual forma, la vegetación y las zonas en preparación y explotación agrícola pueden resultar identificadas usando las bandas SWIR1, NIR y azul.

landsatFCC_agr <- landsat[[c(6,5,2)]]
plotRGB(landsatFCC_agr, axes =TRUE, stretch = "lin", main = "Landsat False Color Composite-Agriculture")

landsatFCC_veghealt <- landsat[[c(5,6,2)]]
plotRGB(landsatFCC_veghealt, axes =TRUE, stretch = "lin", main = "Landsat False Color Composite-
Healthy Veg.")

2.6. Subconjuntos de archivos

La función (subset) y su variación landsat[[c(x1,x2,x3)]] o landsat[[x1:x3]] ya fue empleada anteriormente para graficar los mapas landsatRGB, landsatFCC y se verá posteriormente en los respectivos recortes (landsatcrop), cómoda y practica ya que el orden de enlace puede ser ascendente o descendente teniendo almacenado un solo bloque con todas las variables y no estas independientes y por aparte.

landsatsub1 <- subset(landsat, 1:3)

landsatsub2 <- landsat[[1:3]]

nlayers(landsat)

[1] 11

nlayers(landsatsub2)

[1] 3

Tambien puede ser empleada para sobreescribir un objeto y reemplazar las bandas almacenadas en un bloque, eliminando aquellas que no nos interesan.

landsat <- subset(landsat, 1:7)

names(landsat)

[1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2" "LC08_044034_20170614_B3" "LC08_044034_20170614_B4" "LC08_044034_20170614_B5"
[6] "LC08_044034_20170614_B6" "LC08_044034_20170614_B7"

Renombrandolas, permitirá una facil identificación de estas.

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

[1] "ultra.blue" "blue"       "green"      "red"        "NIR"        "SWIR1"      "SWIR2"     

2.7. Recortes espaciales de imágenes

El caso más común es obtener una imagen satelital de gran tamaño, que abarca el área de interés y mucho más, por lo que se requiere delimitar únicamente la zona donde se desean conocer las características, evitando así ralentizar los procesos.

(extent) requiere dimensionar previamente el área de estudio, y empleando coordenadas cartesianas puede ser delimitado un rectángulo (Ymin, Ymax, Xmin y Xmax), posteriormente la función (crop) permitirá almacenar en un nuevo objeto el recorte realizado entre el raster y el recuadro recientemente creado.

extent(landsat)

class      : Extent 
xmin       : 594090 
xmax       : 639000 
ymin       : 4190190 
ymax       : 4227540 

e <-extent(624384, 635752, 4200047, 4210939)
landsatcrop <- crop(landsat,e)
extent(landsatcrop)

class      : Extent 
xmin       : 624390 
xmax       : 635760 
ymin       : 4200060 
ymax       : 4210950 

2.8. Composición de Colores Verdadero y Falso para las Imágenes Recortadas

A continuación se presentan procesos de ploteo con el detalle de las imágenes recortadas.

landsatRGBcrop <- landsatcrop[[c(4,3,2)]]
plotRGB(landsatRGBcrop, axes =TRUE, stretch = "lin", main = "Landsat True Color Composite-croped")

landsatFCCcrop_veg <- landsatcrop[[c(5,4,3)]]
plotRGB(landsatFCCcrop_veg, axes =TRUE, stretch = "lin", main = "Landsat False Color Composite-croped-vegetation")

2.9. Almacenamiento de Datos en el Disco

Los diferentes archivos en los que se puede almacenar la información procesada, tienen características que deben ser claras.

Archivos “.tif”	archivos “.grd”
Los archivos GeoTiff son de uso común, conservan el orden de almacenamiento de las capas del raster, pero sus nombres se pierden en este tipo de formato	Los archivos de esta característica, permiten almacenar el nombre de las capas, sin embargo, no muchos programas logran leerlos.

x1 <- writeRaster(landsatcrop, filename="data/rs/cropped-landsat.tif", overwrite=TRUE)


x2 <- writeRaster(landsatcrop, filename="data/rs/cropped-landsat.grd", overwrite=TRUE)

2.10. Relación entre bandas

Encontrar la proporcionalidad de la relación entre algunas bandas permite identificar cuales pueden ser usadas sin perder información, como sucede con la banda 1 y 2 de Landsat 8, mediante la función (pairs), la comparación es llevada a cabo, así mismo, puede darse el caso de las bandas 4 y 5, que generan una nube de puntos cercana al eje Y, donde la reflectancia vegetal es muy marcada, mostrando una baja correlación.

pairs(landsatcrop[[1:2]], main = "Ultra-Blue Vs Blue")

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

2.11. Valores de los Pixeles

Teniendo un archivo de polígonos que identifiquen fácilmente las diferentes áreas de explotación agrícola, forestal, urbana entre otras, se pueden extraer coordenadas o puntos mediante la función (spsample), escogiendo la cantidad de puntos y la forma de extracción de estos.

#Esta primera linea, carga el archivo RDS que contiene los polígonos de clasificación de suelo.
lc_samples <- readRDS('data/rs/samples.rds')
#spsample, genera para este caso, 300 puntos de forma regular sobre los polígonos almacenados en el objeto lc_samples
ptsamp <- spsample(lc_samples, 300, type = 'regular')
#Se agreaga la clase de cobertura de suelo a los puntos extraidos. 
ptsamp$class <- over(ptsamp, lc_samples)$class
#Extraer los valores de radiancia en el raster "landsat"de los puntos
df <- extract(landsat, ptsamp)
head(df)

     ultra.blue      blue     green       red       NIR     SWIR1     SWIR2
[1,]  0.1412655 0.1346294 0.1439546 0.1778288 0.3226725 0.3632045 0.2504567
[2,]  0.1400077 0.1235477 0.1126394 0.1163911 0.1883034 0.2332810 0.1958719
[3,]  0.1360607 0.1181911 0.1034443 0.1039431 0.1727325 0.2174066 0.1815155
[4,]  0.1382511 0.1212055 0.1052660 0.1049841 0.1787830 0.2135247 0.1745975
[5,]  0.1370583 0.1364945 0.1499184 0.2007297 0.3393711 0.3296990 0.1946141
[6,]  0.1383378 0.1363426 0.1537352 0.2059128 0.3326266 0.3665659 0.2251053

2.12. Perfiles espectrales

Una de las herramientas mas importantes es el perfil espectral de las bandas contenidas en un objeto, elaborada con los pixeles extraídos anteriormente, donde están representadas algunas características de la superficie.

ms <- aggregate(df, list(ptsamp$class), mean)

rownames(ms) <- ms[,1]
#La primera columna, tendrá los nombres de las clases
ms <- ms[,-1]
ms

mycolor <- c('darkred','yellow','burlywood','cyan','blue')
ms <- as.matrix(ms)

plot(0, ylim=c(0,0.6), xlim=c(1,7),type='n',xlab="Bands",ylab="Reflectance")
for (i in 1:nrow(ms)){
  lines(ms[i,],type="l", lwd=3,lty=1,col=mycolor[i])
}

title(main="Spectral Profile from Landsat", font.main=2)
legend("topleft",rownames(ms),
       cex=0.8, col=mycolor,lty=1,lwd=3,bty="n")

Practica 3: Operaciones matemáticas básicas

Índices de vegetación

vi <- function(img,k,i){
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk-bi)/(bk+bi)
  return(vi)
}
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main="Landsat-NDVI")

Otra alternativa:

vi2 <- function(x,y){
  (x-y)/(x+y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")

Histograma de NDVI

hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "green",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

Índice de Agua de Diferencia Normalizada

wi <- function(img,k,i){
  bk <- img[[k]]
  bi <- img[[i]]
  wi <- (bk-bi)/(bk+bi)
  return(vi)
}
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(topo.colors(10)), main="NDWI-Mc Feeters (2006)")

Histograma de NDWI - Mc Feeters (2006)

hist(ndwi,
     main = "Distribution of NDWI values",
     xlab = "NDWI",
     ylab= "Frequency",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

ndwifunc <- function(img, k, i){
  bcrk <- img[[k]]
  bcri <- img[[i]]
  ndwifunc <- (bcrk - bcri) / (bcrk + bcri)
  return(ndwifunc)
}

mndwi <- ndwifunc(landsat, 3, 6)
plot(mndwi, col = rev(topo.colors(10)), main = "M-NDWI - Xu (2006)")

Histograma de MNDWI - Xu (2006)

hist(mndwi,
     main = "Distribution of M-NDWI values",
     xlab = "MNDWI",
     ylab= "Frequency",
     col = "yellow",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

Índice de Diferencias Normalizadas de Área Construida

(Empleando la imagen recortada)

ndbifunc <- function(img, k, i){
  bcrk <- img[[k]]
  bcri <- img[[i]]
  ndbifunc <- (bcrk - bcri) / (bcri + bcrk)
  return(ndbifunc)
}

ndbi <- ndbifunc(landsatcrop, 6, 5)
plot(ndbi, col = rev(topo.colors(10)), main = "NDBI")

Histograma de NDBI

hist(ndbi,
     main = "Distribution of NDBI values",
     xlab = "NDBI",
     ylab= "Frequency",
     col = "blue",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

Umbralización

veg <- reclassify(ndvi, cbind(-Inf,0.4,NA))
plot(veg, main='Vegetation')

land <- reclassify(ndvi, c(-Inf, 0.25, NA,  0.25, 0.3, 1,  0.3, Inf, NA))
plot(land, main = 'Open Areas')

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
plot(land, add=TRUE, legend=FALSE)

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

La respuesta del perfil espectral, muestra que para las 7 bandas, la reflectancia del agua es la unica que se encuentra entre 0 y 0.2 como valores máximos aproximados, así que acotamos el subconjunto de reclasificación en esta franja.

water <- reclassify(ndwi, cbind(-Inf, 0.1, NA))
#Solo es necesario acotar el limite inferior, ya que todo sobre este valor (0.1) se muestra con clasificación de interés
colors <- colorRampPalette(c("wheat","cyan","blue"))
plot(water,col = colors(10), main = 'water')

Y para el caso de áreas construidas, el nivel de reflectancia muestra el mismo comportamiento sobre los valores de 0.1 o 0.15 hasta 0.2

urbanreclass <- reclassify(ndbi, cbind(-Inf, 0.0, NA))
colors <- colorRampPalette(c("wheat","grey","black"))
plot(urbanreclass,col = colors(50), main = 'Urban zones')

Bajando el límite superior hasta 0, se logran evidenciar incluso, lo que parecen ser vias que rodean los cuerpos de agua.

Análisis de componentes Principales

set.seed(1)
sr <- sampleRandom(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")

pca <- prcomp(sr, scale = TRUE)
pca

Standard deviations (1, .., p=7):
[1] 2.20737373 1.17428246 0.73633573 0.42388639 0.12500714 0.09389890 0.04741454

Rotation (n x k) = (7 x 7):
                 PC1           PC2         PC3         PC4         PC5        PC6         PC7
ultra.blue 0.3662344 -0.4524940727 -0.15415991  0.52595682 -0.21276947  0.1337062  0.54550658
blue       0.4125282 -0.3336052918 -0.15782033  0.09509849 -0.21537031  0.0836529 -0.79447772
green      0.4372284 -0.1056685710 -0.26457625 -0.20443752  0.49716419 -0.6553920  0.09570391
red        0.4333279 -0.0003693985 -0.05074204 -0.67055674  0.01651199  0.5619375  0.20966468
NIR        0.1762472  0.6825385571 -0.58744574  0.32192859  0.09371265  0.2090330 -0.04325163
SWIR1      0.3815340  0.4119308383  0.29523433 -0.08811728 -0.65630314 -0.3878540  0.09301774
SWIR2      0.3743258  0.1929889704  0.66820379  0.33388919  0.47051465  0.1889099 -0.08709876

screeplot(pca)

pci <- predict(landsat, pca, index = 1:2)
plot(pci[[1]])

par(mfrow=c(1,2))
pc2 <- reclassify(pci[[2]], c(-Inf,0,1,0,Inf,NA))
plotRGB(landsatFCC_veg, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plotRGB(landsatFCC_veg, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")

plot(pc2, legend = FALSE, add = TRUE)