Series temporales de imagenes raster. Parte 2

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. Explorar librerías y funciones en R para el tratamiento y análisis de imágenes raster multitemporales, en particular para series temporales de imágenes raster.

Temario

A- Cargar librerias y datos
B-Lectura y visualiación de wRF usando funcion “raster”
C-Lectura y visualiación de wRF usando funcion “brick”
D-Libreria RasterVis
E- Tratar como serie temporal de imagenes
F- Visualización 3D y animación
G- Visualizar series temporales
H-Análisis de Componentes Principales

A-Cargar librerías y datos

Datos para esta práctica Aqui

# library(rgdal)# lectura de archivos vectoriales
library(sf)# tipo de datos vectoriales para R

## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3

library(raster)# lectura y tratamiento de datos raster

## Loading required package: sp

library(rasterVis)#visualización raster multitemporales

## Loading required package: lattice

## Loading required package: latticeExtra

## Loading required package: RColorBrewer

setwd("C:/R_ecohidrologia/seriestemporalesraster") #definir directorio de trabajo 

Lectura de archivo shp de zona de estudio

ecuador<-st_read("provincias_disuelto.shp") #archivo sin shp, fuente de datos: editado INEC, 

## Reading layer `provincias_disuelto' from data source `C:\R_ecohidrologia\seriestemporalesraster\provincias_disuelto.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 490690.4 ymin: 9445216 xmax: 1147852 ymax: 10160820
## epsg (SRID):    32717
## proj4string:    +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs

ecuador<-st_transform(ecuador, 4326)
plot(ecuador["DPA_PROVIN"])

B-Lectura y visualiación de wRF usando funcion “raster”

WRF <- raster("pr_day_Ecuador_CSIRO-Mk3-6-0_2000.nc")#cargar imagen

## Loading required namespace: ncdf4

class(WRF)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

summary(WRF)

##         Precipitation
## Min.        0.0000000
## 1st Qu.     0.9080394
## Median      4.1287465
## 3rd Qu.     8.9990506
## Max.       35.8813210
## NA's        0.0000000

res(WRF) #resolución del raster

## [1] 0.09009783 0.08998303

extent(WRF)

## class       : Extent 
## xmin        : -82.62985 
## xmax        : -74.25075 
## ymin        : -5.74901 
## ymax        : 2.079515

#Sistemas de referencia
proj4string(WRF) 

## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

projection(WRF) 

## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Visualizacion

image(WRF)
plot(as(ecuador,"Spatial"), add = TRUE)

WRF_masked <- mask(WRF, as(ecuador,"Spatial"), snap="out")
image(WRF_masked)
plot(as(ecuador,"Spatial"), add = TRUE)

C- Lectura y visualizacion de WRF con función “brick”

Lectura de raster multitemporal

pr_diaria <- brick("pr_day_Ecuador_CSIRO-Mk3-6-0_2000.nc")#cargar imagen
names(pr_diaria)

##   [1] "X2000.01.01" "X2000.01.02" "X2000.01.03" "X2000.01.04" "X2000.01.05"
##   [6] "X2000.01.06" "X2000.01.07" "X2000.01.08" "X2000.01.09" "X2000.01.10"
##  [11] "X2000.01.11" "X2000.01.12" "X2000.01.13" "X2000.01.14" "X2000.01.15"
##  [16] "X2000.01.16" "X2000.01.17" "X2000.01.18" "X2000.01.19" "X2000.01.20"
##  [21] "X2000.01.21" "X2000.01.22" "X2000.01.23" "X2000.01.24" "X2000.01.25"
##  [26] "X2000.01.26" "X2000.01.27" "X2000.01.28" "X2000.01.29" "X2000.01.30"
##  [31] "X2000.01.31" "X2000.02.01" "X2000.02.02" "X2000.02.03" "X2000.02.04"
##  [36] "X2000.02.05" "X2000.02.06" "X2000.02.07" "X2000.02.08" "X2000.02.09"
##  [41] "X2000.02.10" "X2000.02.11" "X2000.02.12" "X2000.02.13" "X2000.02.14"
##  [46] "X2000.02.15" "X2000.02.16" "X2000.02.17" "X2000.02.18" "X2000.02.19"
##  [51] "X2000.02.20" "X2000.02.21" "X2000.02.22" "X2000.02.23" "X2000.02.24"
##  [56] "X2000.02.25" "X2000.02.26" "X2000.02.27" "X2000.02.28" "X2000.02.29"
##  [61] "X2000.03.01" "X2000.03.02" "X2000.03.03" "X2000.03.04" "X2000.03.05"
##  [66] "X2000.03.06" "X2000.03.07" "X2000.03.08" "X2000.03.09" "X2000.03.10"
##  [71] "X2000.03.11" "X2000.03.12" "X2000.03.13" "X2000.03.14" "X2000.03.15"
##  [76] "X2000.03.16" "X2000.03.17" "X2000.03.18" "X2000.03.19" "X2000.03.20"
##  [81] "X2000.03.21" "X2000.03.22" "X2000.03.23" "X2000.03.24" "X2000.03.25"
##  [86] "X2000.03.26" "X2000.03.27" "X2000.03.28" "X2000.03.29" "X2000.03.30"
##  [91] "X2000.03.31" "X2000.04.01" "X2000.04.02" "X2000.04.03" "X2000.04.04"
##  [96] "X2000.04.05" "X2000.04.06" "X2000.04.07" "X2000.04.08" "X2000.04.09"
## [101] "X2000.04.10" "X2000.04.11" "X2000.04.12" "X2000.04.13" "X2000.04.14"
## [106] "X2000.04.15" "X2000.04.16" "X2000.04.17" "X2000.04.18" "X2000.04.19"
## [111] "X2000.04.20" "X2000.04.21" "X2000.04.22" "X2000.04.23" "X2000.04.24"
## [116] "X2000.04.25" "X2000.04.26" "X2000.04.27" "X2000.04.28" "X2000.04.29"
## [121] "X2000.04.30" "X2000.05.01" "X2000.05.02" "X2000.05.03" "X2000.05.04"
## [126] "X2000.05.05" "X2000.05.06" "X2000.05.07" "X2000.05.08" "X2000.05.09"
## [131] "X2000.05.10" "X2000.05.11" "X2000.05.12" "X2000.05.13" "X2000.05.14"
## [136] "X2000.05.15" "X2000.05.16" "X2000.05.17" "X2000.05.18" "X2000.05.19"
## [141] "X2000.05.20" "X2000.05.21" "X2000.05.22" "X2000.05.23" "X2000.05.24"
## [146] "X2000.05.25" "X2000.05.26" "X2000.05.27" "X2000.05.28" "X2000.05.29"
## [151] "X2000.05.30" "X2000.05.31" "X2000.06.01" "X2000.06.02" "X2000.06.03"
## [156] "X2000.06.04" "X2000.06.05" "X2000.06.06" "X2000.06.07" "X2000.06.08"
## [161] "X2000.06.09" "X2000.06.10" "X2000.06.11" "X2000.06.12" "X2000.06.13"
## [166] "X2000.06.14" "X2000.06.15" "X2000.06.16" "X2000.06.17" "X2000.06.18"
## [171] "X2000.06.19" "X2000.06.20" "X2000.06.21" "X2000.06.22" "X2000.06.23"
## [176] "X2000.06.24" "X2000.06.25" "X2000.06.26" "X2000.06.27" "X2000.06.28"
## [181] "X2000.06.29" "X2000.06.30" "X2000.07.01" "X2000.07.02" "X2000.07.03"
## [186] "X2000.07.04" "X2000.07.05" "X2000.07.06" "X2000.07.07" "X2000.07.08"
## [191] "X2000.07.09" "X2000.07.10" "X2000.07.11" "X2000.07.12" "X2000.07.13"
## [196] "X2000.07.14" "X2000.07.15" "X2000.07.16" "X2000.07.17" "X2000.07.18"
## [201] "X2000.07.19" "X2000.07.20" "X2000.07.21" "X2000.07.22" "X2000.07.23"
## [206] "X2000.07.24" "X2000.07.25" "X2000.07.26" "X2000.07.27" "X2000.07.28"
## [211] "X2000.07.29" "X2000.07.30" "X2000.07.31" "X2000.08.01" "X2000.08.02"
## [216] "X2000.08.03" "X2000.08.04" "X2000.08.05" "X2000.08.06" "X2000.08.07"
## [221] "X2000.08.08" "X2000.08.09" "X2000.08.10" "X2000.08.11" "X2000.08.12"
## [226] "X2000.08.13" "X2000.08.14" "X2000.08.15" "X2000.08.16" "X2000.08.17"
## [231] "X2000.08.18" "X2000.08.19" "X2000.08.20" "X2000.08.21" "X2000.08.22"
## [236] "X2000.08.23" "X2000.08.24" "X2000.08.25" "X2000.08.26" "X2000.08.27"
## [241] "X2000.08.28" "X2000.08.29" "X2000.08.30" "X2000.08.31" "X2000.09.01"
## [246] "X2000.09.02" "X2000.09.03" "X2000.09.04" "X2000.09.05" "X2000.09.06"
## [251] "X2000.09.07" "X2000.09.08" "X2000.09.09" "X2000.09.10" "X2000.09.11"
## [256] "X2000.09.12" "X2000.09.13" "X2000.09.14" "X2000.09.15" "X2000.09.16"
## [261] "X2000.09.17" "X2000.09.18" "X2000.09.19" "X2000.09.20" "X2000.09.21"
## [266] "X2000.09.22" "X2000.09.23" "X2000.09.24" "X2000.09.25" "X2000.09.26"
## [271] "X2000.09.27" "X2000.09.28" "X2000.09.29" "X2000.09.30" "X2000.10.01"
## [276] "X2000.10.02" "X2000.10.03" "X2000.10.04" "X2000.10.05" "X2000.10.06"
## [281] "X2000.10.07" "X2000.10.08" "X2000.10.09" "X2000.10.10" "X2000.10.11"
## [286] "X2000.10.12" "X2000.10.13" "X2000.10.14" "X2000.10.15" "X2000.10.16"
## [291] "X2000.10.17" "X2000.10.18" "X2000.10.19" "X2000.10.20" "X2000.10.21"
## [296] "X2000.10.22" "X2000.10.23" "X2000.10.24" "X2000.10.25" "X2000.10.26"
## [301] "X2000.10.27" "X2000.10.28" "X2000.10.29" "X2000.10.30" "X2000.10.31"
## [306] "X2000.11.01" "X2000.11.02" "X2000.11.03" "X2000.11.04" "X2000.11.05"
## [311] "X2000.11.06" "X2000.11.07" "X2000.11.08" "X2000.11.09" "X2000.11.10"
## [316] "X2000.11.11" "X2000.11.12" "X2000.11.13" "X2000.11.14" "X2000.11.15"
## [321] "X2000.11.16" "X2000.11.17" "X2000.11.18" "X2000.11.19" "X2000.11.20"
## [326] "X2000.11.21" "X2000.11.22" "X2000.11.23" "X2000.11.24" "X2000.11.25"
## [331] "X2000.11.26" "X2000.11.27" "X2000.11.28" "X2000.11.29" "X2000.11.30"
## [336] "X2000.12.01" "X2000.12.02" "X2000.12.03" "X2000.12.04" "X2000.12.05"
## [341] "X2000.12.06" "X2000.12.07" "X2000.12.08" "X2000.12.09" "X2000.12.10"
## [346] "X2000.12.11" "X2000.12.12" "X2000.12.13" "X2000.12.14" "X2000.12.15"
## [351] "X2000.12.16" "X2000.12.17" "X2000.12.18" "X2000.12.19" "X2000.12.20"
## [356] "X2000.12.21" "X2000.12.22" "X2000.12.23" "X2000.12.24" "X2000.12.25"
## [361] "X2000.12.26" "X2000.12.27" "X2000.12.28" "X2000.12.29" "X2000.12.30"

class(pr_diaria)

## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"

Visualizacion

image(pr_diaria[[1]])

plot(pr_diaria[[1]])

plot(pr_diaria[[1:4]])

D-Libreria RasterVis

levelplot(pr_diaria[[1]]) #mean # funcion de rasterVis

Color adecuado a precipitación

mytheme <- rasterTheme(region = brewer.pal(9, "Blues"))
levelplot(pr_diaria[[1]], par.settings = mytheme)

levelplot(pr_diaria[[1]], par.settings = mytheme, margin=F)#eliminar los margenes

levelplot(pr_diaria[[1]], par.settings = mytheme, 
          margin = list(FUN = 'sd'))#plotar en margenes otra funcion

levelplot(pr_diaria[[1]], par.settings = mytheme, 
          margin = list(FUN = 'mean'),
         contour = TRUE) # Isolineas y promedio en margenes

levelplot(pr_diaria[[1]], par.settings = mytheme, 
          margin = list(FUN = 'mean'),
          contour = TRUE)+ layer(sp.polygons(as(ecuador,"Spatial"), lwd=2, col='black')) # Plotar con poligono de ecuador, + layer es funcion de libreria lattice

Plotear imagenes multitemporales

levelplot(pr_diaria[[1:4]])

levelplot(pr_diaria[[1:3]], par.settings = mytheme,  
                colorkey = list(space = "bottom"),                
                layout = c(3, 1, 1)) + layer(sp.polygons(as(ecuador,"Spatial")))

levelplot(pr_diaria[[1:9]], par.settings = mytheme,  
          colorkey = list(space = "bottom"),                
          layout = c(3, 3, 1)) + layer(sp.polygons(as(ecuador,"Spatial")))

E- Tratar como serie temporal de imagenes

library(rts) # raster time serie, requiere xts y zoo

## Warning: package 'rts' was built under R version 3.4.4

## Loading required package: xts

## Warning: package 'xts' was built under R version 3.4.4

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 3.4.4

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: RCurl

## Warning: package 'RCurl' was built under R version 3.4.4

## Loading required package: bitops

## rts 1.0-45 (2018-03-17)

d<- seq(as.Date("2000-01-01"), as.Date("2000-12-30"), by="day")
d <- as.Date(d)
rt <- rts(pr_diaria,d)
rt[[1]]

## class       : RasterLayer 
## band        : 1  (of  365  bands)
## dimensions  : 87, 93, 8091  (nrow, ncol, ncell)
## resolution  : 0.09009783, 0.08998303  (x, y)
## extent      : -82.62985, -74.25075, -5.74901, 2.079515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\R_ecohidrologia\seriestemporalesraster\pr_day_Ecuador_CSIRO-Mk3-6-0_2000.nc 
## names       : X2000.01.01 
## z-value     : 2000-01-01 
## zvar        : precip

summary(rt)

##        Length         Class          Mode 
##             1 RasterBrickTS            S4

plot(rt[[1:4]])

write.rts(rt,"rtfile", overwrite = TRUE) # writing into the working directory
# rt <- read.rts("rtfile") # reading from the working directory

#Generar medias mensuales de los rasters
mes <- endpoints(rt,'months') # extract the end index on each year period
mes

##  [1]   0  31  60  91 121 152 182 213 244 274 305 335 365

rt.m <- period.apply(rt,mes,mean) 

plot(rt.m)

levelplot(rt.m@raster)

levelplot(rt.m@raster, par.settings = mytheme, margin = FALSE)+ layer(sp.polygons(as(ecuador,"Spatial")))

names(rt.m@raster) <- c(1:12)
levelplot(rt.m@raster, par.settings = mytheme, margin = FALSE)+ layer(sp.polygons(as(ecuador,"Spatial")))

Mascara para ecuador

rt.m.c <- mask(rt.m@raster, as(ecuador,"Spatial"))
levelplot(rt.m.c[[1]], par.settings = mytheme, margin = FALSE)+ layer(sp.polygons(as(ecuador,"Spatial")))

Exploracion valores del raster

histogram(rt.m.c)

histogram(rt.m.c[[1]])

densityplot(rt.m.c)

bwplot(rt.m.c, violin=FALSE)

bwplot(rt.m.c, violin=TRUE)  # boxplot and a kernel density plot

Hovmoller. A Hovmöller diagram is a commonly used way of plotting meteorological data to highlight the role of waves. The axes of a Hovmöller diagram are typically longitude or latitude (abscissa or x-axis) and time (ordinate or y-axis) with the value of some field represented through color or shading.

rt.m.c <- setZ(rt.m.c, 1:12, name='time')
getZ(rt.m.c)

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12

hovmoller(rt.m.c)

hovmoller(rt.m.c, xlab='Latitude', FUN=mean)#Latitude as default

hovmoller(rt.m.c, xlab='Latitude', FUN=sd)#Latitude as default

F- Visualización 3D y animación

library(rgl)

## Warning: package 'rgl' was built under R version 3.4.4

plot3D(rt.m.c[[1]], col=heat.colors, rev=TRUE, useLegend=TRUE, adjust=TRUE)

Animación. No funciona en markdown, pero si en script de R.

# animate(rt.m@raster, pause=0.25, main="Animación", n=1)

G- Visualizar series temporales

Agregado mensual

rt.s <- period.apply(rt,mes,sum) 
rt.s <- mask(rt.s@raster, as(ecuador,"Spatial"))
plot(rt.s)

sample <- sampleRandom(rt.s, size=10, cells=TRUE, sp=TRUE)
sample

## class       : SpatialPointsDataFrame 
## features    : 10 
## extent      : -80.51255, -75.64727, -2.554612, 0.1448793  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : cell,      layer.1.1.1,      layer.2.1.1,      layer.1.2.1,      layer.2.2.1,      layer.1.1.2,       layer.2.1.2,      layer.1.2.2,      layer.2.2.2,        layer.1.1,       layer.2.1,           layer.1.2,        layer.2.2 
## min values  : 2027, 32.5315439067781, 75.7004984617233, 179.399980321527, 167.942185103893, 16.2649090441409, 0.404609836929012,                0,                0,                0,               0, 0.00761787407100201, 6.86941778846085 
## max values  : 4799, 172.808274716139, 312.974317312241, 519.445694819093, 583.775461673737,   216.7229834795,  218.253309011459, 27.9461622685194, 58.9670548886061, 139.645025104284, 286.64189812867,    210.965171545744, 215.367280811071

names(sample)

##  [1] "cell"        "layer.1.1.1" "layer.2.1.1" "layer.1.2.1" "layer.2.2.1"
##  [6] "layer.1.1.2" "layer.2.1.2" "layer.1.2.2" "layer.2.2.2" "layer.1.1"  
## [11] "layer.2.1"   "layer.1.2"   "layer.2.2"

names(sample) <- c("cell", 1:12)

## Warning in checkNames(value): attempt to set invalid names: this may lead
## to problems later on. See ?make.names

names(sample)

##  [1] "cell" "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"   
## [11] "10"   "11"   "12"

head(sample@data)

##   cell         1        2        3        4         5           6
## 1 4488  32.53154 213.9163 295.6176 220.8481  16.26491   0.4046098
## 2 3572 149.74809 312.9743 519.4457 529.5451 123.12730  47.2477319
## 3 2566 124.30773 150.5755 213.2002 229.0522 139.53071 191.8188863
## 4 4799 172.80827 213.7065 411.6247 417.9084 118.48879 135.0662365
## 5 2027 163.74062 193.6031 261.5027 219.1627 146.91337 218.2533090
## 6 3147 144.15452 203.9573 242.5269 234.3999 155.57775 152.9813079
##           7         8         9        10           11         12
## 1  0.000000  0.000000   0.00000   0.00000 7.617874e-03   6.869418
## 2  1.807479  5.930349  15.42963  11.62776 3.018348e+01 134.681461
## 3 16.881745 36.630706  73.85829  98.45021 1.438681e+02 192.848706
## 4 27.946162 17.115015 139.64503 141.94953 1.577578e+02 193.108377
## 5 23.624356 58.967055  67.89457 225.77705 2.109652e+02 177.928950
## 6 19.411887 21.699736  99.82526 286.64190 1.908095e+02 154.137128

plot(t(sample@data[2,2:13]), type="l", col="red")
lines(t(sample@data[3,2:13]), type="o", col="blue")
lines(t(sample@data[4,2:13]), type="o", col="green")

pts.s <- list("sp.points", sample, col="red", pch=16)
spplot(as(ecuador,"Spatial"), "DPA_PROVIN",sp.layout = list(pts.s))

Libreria dygraphs

library(dygraphs)

## Warning: package 'dygraphs' was built under R version 3.4.4

dysample <- data.frame( c(1:12), t(sample@data[,2:13]))

dygraph(dysample)

dygraph(dysample[1:4])

dygraph(dysample) %>% dyRangeSelector()

Libreria highchart

library(highcharter)

## Warning: package 'highcharter' was built under R version 3.4.4

## Highcharts (www.highcharts.com) is a Highsoft software product which is

## not free for commercial and Governmental use

#convertir a serie temporal
myts <- ts(dysample[,2:11] ,start=c(2000, 1), end=c(2000, 12), frequency=12)

highchart() %>% 
  hc_add_series(myts[,"X1"], name = "X1") %>%
  hc_add_series(myts[,"X2"], name = "X2")%>%
  hc_add_series(myts[,"X3"], name = "X3")

H-Análisis de Componentes Principales

library(RStoolbox)

## Warning: package 'RStoolbox' was built under R version 3.4.4

rpc<- rasterPCA(rt.m@raster, nComp = nlayers(rt.m@raster), spca = TRUE, maskCheck = FALSE)
rpc

## $call
## rasterPCA(img = rt.m@raster, nComp = nlayers(rt.m@raster), spca = TRUE, 
##     maskCheck = FALSE)
## 
## $model
## Call:
## princomp(cor = spca, covmat = covMat[[1]])
## 
## Standard deviations:
##     Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6 
## 2.79534395 1.40281899 0.94331808 0.69696889 0.51772111 0.45531770 
##     Comp.7     Comp.8     Comp.9    Comp.10    Comp.11    Comp.12 
## 0.40902077 0.28539647 0.24174019 0.17892048 0.13425973 0.09980677 
## 
##  12  variables and  8091 observations.
## 
## $map
## class       : RasterBrick 
## dimensions  : 87, 93, 8091, 12  (nrow, ncol, ncell, nlayers)
## resolution  : 0.09009783, 0.08998303  (x, y)
## extent      : -82.62985, -74.25075, -5.74901, 2.079515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       :        PC1,        PC2,        PC3,        PC4,        PC5,        PC6,        PC7,        PC8,        PC9,       PC10,       PC11,       PC12 
## min values  : -4.1648510, -3.0809120, -1.6985005, -3.7922670, -1.7089155, -1.6678382, -1.0543438, -2.1382520, -1.0378210, -0.5795459, -0.4802221, -0.5983044 
## max values  :  7.4514435,  5.3617540,  3.5118317,  2.3967042,  1.2708613,  1.2286601,  1.8132104,  0.9542458,  1.0663749,  0.5912115,  0.6119950,  0.4232415 
## 
## 
## attr(,"class")
## [1] "rasterPCA" "RStoolbox"

summary(rpc$model)

## Importance of components:
##                           Comp.1    Comp.2     Comp.3     Comp.4
## Standard deviation     2.7953440 1.4028190 0.94331808 0.69696889
## Proportion of Variance 0.6511623 0.1639918 0.07415408 0.04048047
## Cumulative Proportion  0.6511623 0.8151541 0.88930816 0.92978863
##                            Comp.5     Comp.6    Comp.7      Comp.8
## Standard deviation     0.51772111 0.45531770 0.4090208 0.285396475
## Proportion of Variance 0.02233626 0.01727618 0.0139415 0.006787596
## Cumulative Proportion  0.95212489 0.96940108 0.9833426 0.990130173
##                            Comp.9     Comp.10    Comp.11      Comp.12
## Standard deviation     0.24174019 0.178920480 0.13425973 0.0998067659
## Proportion of Variance 0.00486986 0.002667712 0.00150214 0.0008301159
## Cumulative Proportion  0.99500003 0.997667745 0.99916988 1.0000000000

ggRGB(rpc$map,1,2,3, stretch="lin", q=0)

if(require(gridExtra)){
plots <- lapply(1:3, function(x) ggR(rpc$map, x, geom_raster = TRUE))
grid.arrange(plots[[1]],plots[[2]], plots[[3]], ncol=2)
}

## Loading required package: gridExtra