Introduccion

Las imagenes nos han pemitido cada vez mas monitorear y mejorar la gestion de los recursos naturales del planeta, ademas de ser valiosos aliados en investigacion y cuantificacion sobre los mismos. El presente texto nos permite identificar algunas formas basicas de exploracion de las mismas.

library(raster)
## Loading required package: sp
library(sp)

#blue
b2<-raster("C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B1.tif")# Se crea el layer por cada banda
#green
b3<-raster("C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B2.tif")#creacion layer para banda verde
#red
b4<-raster("C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B4.tif")
#infrarrojo cercano
b5<-raster("C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B5.tif")
b2 #Permite visualizar la infomacion de la capa
## 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:/Files_R/rsdata/rs/LC08_044034_20170614_B1.tif 
## names      : LC08_044034_20170614_B1 
## values     : 0.09641791, 0.7346282  (min, max)
crs(b2)#Permite identificar caractristicas del sistema de referencia
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
ncell(b2)#Permite identificar el numero de celdas de la imagen 
## [1] 1863765
dim(b2)#dimensiones del pixel #resolucion espacial
## [1] 1245 1497    1
res(b2)#permite identificar numero de bandas
## [1] 30 30
nlayers(b2)
## [1] 1
compareRaster(b2,b3)#las bandas tienen la misma resolucion, origen
## [1] TRUE
s<- stack(b5,b4,b3)#crear un rasterstack(UN OBJETO CON MULTIPLES LAYERS)desde una capa singular 
s# se observan las propiedades de rasterstack
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (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_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B2 
## min values :            0.0008457669,            0.0208406653,            0.0748399049 
## max values :               1.0124315,               0.7861769,               0.7177562
filenames<-paste0("C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B",1:11,".tif")
filenames
##  [1] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B1.tif" 
##  [2] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B2.tif" 
##  [3] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B3.tif" 
##  [4] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B4.tif" 
##  [5] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B5.tif" 
##  [6] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B6.tif" 
##  [7] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B7.tif" 
##  [8] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B8.tif" 
##  [9] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B9.tif" 
## [10] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B10.tif"
## [11] "C:\\Files_R\\rsdata\\rs\\LC08_044034_20170614_B11.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

las capas representan la intensidad de reflecion de las siguientes longitudes de onda,Ultra blue, blue, green, red, near infrared (TIRS) 2, no podriamos usar

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))

La leyenda de los mapas puede tomar un rango entre o - 1, la diferencia de las sombras y el rango de leyenda entre las diferentes bandas es por que la superficie cada capa representa un valor de rango de longitud de onda particular por ejemplo la vegetacion refleja mas energia en NIR que en otra longitud de onda

landsatRGB <- stack (b4,b5,b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "landsat True Color Composite")

La opcion true color o color natural , esto es que algunas de las vistas de de las fotografias necesitan bandas en el rojo, para una imagen landsat la banda 4 (rojo), 3 verde (green), y 2 azul,plotRGB metodo que puede usar la combinacion de ellas en una composicion singular. Se puede usar el argumento adicional plotRGB improvisar la visualizacion estira la linea de valores usnado strecth = “lin” # El color verdadero revela mucho mas acerca del paisaje que la imgen en tonos grises # otra visualizacion es el falso color esta combina las bandas NIR, rojo, y verde es popular por que visualiza la vegetacion mejor.

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" )

Nota:siempre verificar a documentacion del paquete
help(plotRGB)

landsatsub1<- subset(landsat, 1:3)#seleccionar las primeras tres bandas
landsatsub2 <- landsat[[1:3]]# consruye un subconjunto de las bandas 1 hasta la 3
nlayers(landsat) #numero de capas en el stack 
## [1] 11
help(plotRGB)
## starting httpd help server ... done
nlayers(landsatsub1)# Numero de capas den el subconjunto landsatsub1
## [1] 3

Mediante la funcion crop podemos crear subconjuntos para limitar el analisis

landsat <- subset(landsat, 1:7)
names(landsat)
## [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"
names(landsat)<- c("ultra-blue", "blue","Green", "red","NIR", "SWIR1", "SWIR2" )
names(landsat)
## [1] "ultra.blue" "blue"       "Green"      "red"        "NIR"       
## [6] "SWIR1"      "SWIR2"

Subconjuntos de datos y renombrar bandas

extent(landsat) # Podemos  con la funcion crop usando un objeto extent, extraer un limite o un espacio geografico
## class      : Extent 
## xmin       : 594090 
## xmax       : 639000 
## ymin       : 4190190 
## ymax       : 4227540
e<- extent(624387,635752,4200047,4210939)
landsatcrop<-crop(landsat,e) # Generando zona de interés

Guaradar Resultados en el disco

x<- writeRaster (landsatcrop, filename = "cropped-landsat.tif", overwrite=TRUE)#usando writerRaster para guardar

writerRaster(landsatcrop,filename=“cropped-landsat.grd”,overwrite=TRUE)# usando el formato raster-grd ventaja guarda los nombres de las capas, desventaja no es de facil lectura como elformato Geotiff por ejemplo.

pairs(landsatcrop[[1:2]],main =“ultra-blue versus blue”)# mediante una matriz de dispersion se pueden explorar las relaciones entre lasdiferentes capas del raster, esto se realiza con la funcion pairs la cual pertenece al paquete raster.

Relacion entre bandas

Una matriz de dispersion puede ayudar a explorar la relacion entre las capas esta se puede realizar con la funcion pairs del paquete raster.

pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")# generando grafica comparativa

Ahora ploteamos la refleccion de la longitud de onda de la banda roja contra la longitud de onda de la banda NIR

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

Extrayendo valores del pixel

samp <- readRDS("C:\\Users\\yarvis\\Music\\Files\\rsdata\\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.1350848 0.1182778 0.1069575 0.1081720 0.1722988 0.2221125 0.1884552
## [2,]  0.1351933 0.1185164 0.1045937 0.1040732 0.1973249 0.2271872 0.1883467
## [3,]  0.1319186 0.1280151 0.1386414 0.1845299 0.3165787 0.3147353 0.1835107
## [4,]  0.1422847 0.1410703 0.1541255 0.2053706 0.3348170 0.3712935 0.2343654
## [5,]  0.1524557 0.1567930 0.1787180 0.2424761 0.3958426 0.4278518 0.2684998
## [6,]  0.1370149 0.1405281 0.1642965 0.2321750 0.3678454 0.3849343 0.2171030
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.1768226 0.16683270 0.16598042 0.17885783 0.24130394 0.22307757
## cropland  0.1119888 0.08976627 0.08398548 0.05419306 0.46686841 0.15577980
## fallow    0.1318800 0.11590648 0.10285030 0.11147207 0.17197724 0.22916015
## open      0.1398763 0.13925354 0.15532887 0.21085210 0.35093515 0.36149541
## water     0.1335521 0.11657636 0.09921908 0.07862744 0.04873286 0.03336133
##               SWIR2
## built    0.18181153
## cropland 0.06986625
## fallow   0.19461506
## open     0.21561460
## water    0.02700576
# 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 (main ="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
      cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")