Imagen LandSat 8, Sibate 2015

Landsat 8 recopilada el 04 de enero de 2015. El subconjunto cubre el area de Sibate, departamento de Cundinamarca, en Colombia

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Seleccionar la carpeta de destino

#C:/Users/David/Desktop/IMAGENES LANDSAT SIBATE
getwd()

[1] "C:/Users/David/Desktop/Landsat_procesamiento"

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Crear objetos RasterLayer para capas individuales Landsat (bandas) y las llamamos enrutandolas como por ejemplo (‘./data/LC08_L1TP_008057_20180317_20180403_01_T1_B2.tif’)

IMPORTANTE: instalar paquetes: raster, sp, rgdal

En este capitulo se describe la forma de acceder y explorar por satelite datos de teledeteccion con R . Tambien mostramos como usarlos para hacer mapas.

Utilizaremos principalmente un subconjunto espacial de una escena Landsat 8 recopilada el 14 de junio de 2017. El subconjunto cubre el area entre Concord y Stockton , en California, EE. UU.

Todas las escenas de imagenes de Landsat tienen una identificacion de producto y metadatos unicos. Puede encontrar la informacion sobre el sensor Landsat, el satelite, la ubicacion en la Tierra (ruta WRS, fila WRS) y la fecha de adquisicion a partir de la ID del producto. Por ejemplo, el identificador de producto de los datos que utilizaremos es ‘LC08_044034_20170614’. Con base en esta guia , puede ver que el Sensor-Satelite es OLI / TIRS combinados Landsat 8, WRS Path 44, WRS Row 34 y recopilados el 14 de junio de 2017. Las escenas de Landsat se entregan mas comunmente como un archivo comprimido, que contiene archivos separados para cada banda

Comenzaremos explorando y visualizando los datos (consulte las instrucciones en el Capitulo 1 para obtener instrucciones de descarga de datos si aun no lo ha hecho).

TITULO: PROPIEDADES DE IMAGEN

Crear objetos RasterLayer para capas individuales Landsat (bandas)

library(raster)
b1 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B1.tif')
# azul
b2 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B2.tif')
# verde
b3 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B3.tif')
# rojo
b4 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B4.tif')
# Infrarojo cercano (NIR)
b5 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B5.tif')
b6 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B6.tif')
b7 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B7.tif')
b8 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B8.tif')
b9 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B9.tif')
b10 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B10.tif')
b11 = raster('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B11.tif')

se inserta un nuevo bloque de codigo ctrl+alt+i: imprimir las variables para verificar, por ejemplo:

b2

class      : RasterLayer 
dimensions : 7741, 7581, 58684521  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 446085, 673515, 362985, 595215  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/David/Desktop/Landsat_procesamiento/data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B2.tif 
names      : LC08_L1TP_008057_20150104_20170415_01_T1_B2 
values     : 0, 65535  (min, max)

Puede ver la resolucion espacial, la extension, el numero de capas, el sistema de referencia de coordenadas y mas.

TITULO: INFORMACION DE IMAGEN Y ESTADISTICAS

A continuacion se muestra como puede acceder a varias propiedades desde un objeto Raster * (esto es lo mismo para cualquier conjunto de datos raster).

# coordinate reference system (CRS)
crs(b2)

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

## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# Number of cells, rows, columns
ncell(b2)

[1] 58684521

## [1] 58684521
dim(b2)

[1] 7741 7581    1

## [1] 7741 7581    1
# spatial resolution
res(b2)

[1] 30 30

## [1] 30 30
# Number of bands
nlayers(b2)

[1] 1

## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)

[1] TRUE

## [1] TRUE

Puede crear un RasterStack (un objeto con varias capas) a partir de los objetos RasterLayer (banda unica) existentes

Las bandas que llamos son las mismas que hemos llamadodo arriba como RasterLayer

s2015 = stack(b5, b4, b3)
# Check the properties of the RasterStack
s2015

class      : RasterStack 
dimensions : 7741, 7581, 58684521, 3  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 446085, 673515, 362985, 595215  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : LC08_L1TP_008057_20150104_20170415_01_T1_B5, LC08_L1TP_008057_20150104_20170415_01_T1_B4, LC08_L1TP_008057_20150104_20170415_01_T1_B3 
min values :                                           0,                                           0,                                           0 
max values :                                       65535,                                       65535,                                       65535 

## class      : RasterStack
## dimensions : 7741, 7581, 58684521, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 446085, 673515, 362985, 595215  (xmin, xmax, ymin, ymax)
## crs        : ++proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_L1TP_008057_20180317_20180403_01_T1_B5, LC08_L1TP_008057_20180317_20180403_01_T1_B4,     LC08_L1TP_008057_20180317_20180403_01_T1_B3
## min values :  0, 0, 0 
## max values :  65535, 65535, 65535

Tambien puede crear el RasterStack usando los nombres de archivo

# first create a list of raster layers to use
filenames = paste0('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B', 1:7, ".tif")
filenames

[1] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B1.tif"
[2] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B2.tif"
[3] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B3.tif"
[4] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B4.tif"
[5] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B5.tif"
[6] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B6.tif"
[7] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B7.tif"

## [1] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B1.tif" 
## [2] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B2.tif" 
## [3] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B3.tif" 
## [4] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B4.tif" 
## [5] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B5.tif" 
## [6] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B6.tif" 
## [7] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B7.tif" 
###### no lee la banda 8 porque tiene diferente extension
## [8] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B8.tif" 
## [9] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B9.tif" 
## [10] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B10.tif"
## [11] "./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B11.tif"
landsat2015= stack(filenames)
landsat2015

class      : RasterStack 
dimensions : 7741, 7581, 58684521, 7  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 446085, 673515, 362985, 595215  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : LC08_L1TP//5_01_T1_B1, LC08_L1TP//5_01_T1_B2, LC08_L1TP//5_01_T1_B3, LC08_L1TP//5_01_T1_B4, LC08_L1TP//5_01_T1_B5, LC08_L1TP//5_01_T1_B6, LC08_L1TP//5_01_T1_B7 
min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
max values :                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535 

## class      : RasterStack
## dimensions : 7741, 7581, 58684521, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 446085, 673515, 362985, 595215  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_L1TP//5_01_T1_B1, LC08_L1TP//5_01_T1_B2, LC08_L1TP//5_01_T1_B3, LC08_L1TP//5_01_T1_B4, LC08_L1TP//5_01_T1_B5, LC08_L1TP//5_01_T1_B6, LC08_L1TP//5_01_T1_B7
## min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values :                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535 

Arriba creamos un RasterStack con 11 capas. Las capas representan la intensidad de la reflexion en las siguientes longitudes de onda: Ultra azul, azul, verde, rojo, infrarrojo cercano (NIR), infrarrojo de onda corta (SWIR) 1, infrarrojo de onda corta (SWIR) 2, pancromatico, cirro, infrarrojo termico (TIRS) 1, Infrarrojo termico (TIRS) 2. No utilizaremos las ultimas cuatro capas y vera como eliminarlas en las siguientes secciones.

TITULO: BANDA UNICA Y MAPAS COMPUESTOS

Puede trazar capas individuales de un RasterStack de una imagen multiespectral.

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

Echa un vistazo a las leyendas de los mapas creados anteriormente. Pueden variar entre 0 y 1. Observe la diferencia en el sombreado y el rango de leyendas entre las diferentes bandas. Esto se debe a que las diferentes caracteristicas de la superficie reflejan la radiacion solar incidente de manera diferente. Cada capa representa la cantidad de radiacion solar incidente que se refleja para un rango de longitud de onda particular. Por ejemplo, la vegetacion refleja mas energia en NIR que otras longitudes de onda y, por lo tanto, parece mas brillante. Por el contrario, el agua absorbe la mayor parte de la energia en la longitud de onda NIR y parece oscura.

No obtenemos tanta informacion de estas parcelas en escala de grises; a menudo se combinan en un “compuesto” para crear tramas mas interesantes. Puede obtener mas informacion sobre los compuestos de color en la teledeteccion aqui y tambien en la seccion a continuacion.

Para hacer una imagen de “color verdadero (o natural)”, es decir, algo que se parece a una fotografia normal (vegetacion en verde, azul agua, etc.), necesitamos bandas en las regiones roja, verde y azul. Para esta imagen Landsat, se pueden usar las bandas 4 (rojo), 3 (verde) y 2 (azul). El - plotRGB - metodo se puede utilizar para combinarlos en un solo compuesto. Tambien puede proporcionar argumentos adicionales - plotRGB - para mejorar la visualizacion (por ejemplo, un estiramiento lineal de los valores, utilizando ).- strecth = “lin” -

landsatRGB2015 = stack(b4, b3, b2)
plotRGB(landsatRGB2015, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

El compuesto de color verdadero revela mucho mas sobre el paisaje que las imagenes grises anteriores. Otro metodo popular de visualizacion de imagenes en la teledeteccion es la imagen conocida de “color falso” en la que se combinan las bandas NIR, rojo y verde. Esta representacion es popular ya que hace que sea facil ver la vegetacion (en rojo).

par(mfrow = c(1,2))
plotRGB(landsatRGB2015, axes=TRUE, stretch="lin", main="Landsat True Color Composite")
landsatFCC2015 = stack(b5, b4, b3)
plotRGB(landsatFCC2015, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

Nota : Compruebe siempre la documentacion del paquete (- help(plotRGB) -) para ver otros argumentos que se pueden agregar (como la escala) para mejorar o modificar la imagen.

Pregunta 1 : Use la funcion plotRGB con RasterStack `` landsat ’’ para crear un compuesto de color verdadero y falso (recuerde la posicion de las bandas en la pila).

TITULO: SUBCONJUNTO Y RENOMBRAR BANDAS

Puede seleccionar capas (bandas) especificas mediante la - subset - funcion o mediante indexacion.

# select first 3 bands only
landsat2015sub01 = subset(landsat2015, 1:3)
# same
landsat2015sub02 = landsat2015[[1:3]]
# Number of bands in the original and new data
nlayers(landsat2015)

[1] 7

## [1] 7
nlayers(landsat2015sub01)

[1] 3

## [1] 3
nlayers(landsat2015sub02)

[1] 3

## [1] 3

No usaremos las ultimas cuatro bandas - landsat -. Puedes eliminar aquellos usando:

landsat2015 = subset(landsat2015, 1:7)
plot(landsat2015)

Para mayor claridad, es util establecer los nombres de las bandas.

names(landsat2015)

[1] "LC08_L1TP_008057_20150104_20170415_01_T1_B1" "LC08_L1TP_008057_20150104_20170415_01_T1_B2"
[3] "LC08_L1TP_008057_20150104_20170415_01_T1_B3" "LC08_L1TP_008057_20150104_20170415_01_T1_B4"
[5] "LC08_L1TP_008057_20150104_20170415_01_T1_B5" "LC08_L1TP_008057_20150104_20170415_01_T1_B6"
[7] "LC08_L1TP_008057_20150104_20170415_01_T1_B7"

## [1] "LC08_L1TP_008057_20180317_20180403_01_T1_B1" "LC08_L1TP_008057_20180317_20180403_01_T1_B2"
## [3] "LC08_L1TP_008057_20180317_20180403_01_T1_B3" "LC08_L1TP_008057_20180317_20180403_01_T1_B4"
## [5] "LC08_L1TP_008057_20180317_20180403_01_T1_B5" "LC08_L1TP_008057_20180317_20180403_01_T1_B6"
## [7] "LC08_L1TP_008057_20180317_20180403_01_T1_B7"
names(landsat2015) = c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat2015)

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

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

TITULO: SUBCONJUNTO ESPACIAL O RECORTE

El subconjunto espacial se puede usar para limitar el analisis a un subconjunto geografico de la imagen. Los subconjuntos espaciales se pueden crear con la - crop - funcion, utilizando un - extent - objeto u otro objeto espacial del que se puede extraer una Extension.

# Using extent
extent(landsat2015)

class      : Extent 
xmin       : 446085 
xmax       : 673515 
ymin       : 362985 
ymax       : 595215 

## class      : Extent
## xmin       : 446085 
## xmax       : 673515
## ymin       : 362985
## ymax       : 595215
sibate2015 = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
e = extent(sibate2015)
# crop landsat by the extent
landsatcrop2015 = crop(landsat2015, e)
plotRGB(landsatcrop2015, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

Pregunta 2 : Tambien es posible la seleccion interactiva de la imagen. Use drawExtent y drawPoly para seleccionar un area de interes

Pregunta 3 : Use el `` Landsatcrop ’’ de RasterStack para crear un compuesto de color verdadero y falso

TITULO: GUARDANDO RESULTADOS EN EL DISCO

En esta etapa, es posible que queramos guardar el raster en el disco usando la funcion - writeRaster -. Se admiten multiples tipos de archivos. Utilizaremos el formato GeoTiff de uso comun. Mientras se conserva el orden de las capas, los nombres de las capas se pierden desafortunadamente en el formato GeoTiff

x = writeRaster(landsatcrop2015, filename="cropped-landsat.tif", overwrite=TRUE)
plot(x)

Alternativamente, puede utilizar el formato ‘raster-grd’.

writeRaster(landsatcrop2015, filename="cropped-landsat.grd", overwrite=TRUE)

Una ventaja de este formato es que guarda los nombres de las capas. La desventaja de este formato es que no muchos otros programas pueden leer los datos, en contraste con el formato GeoTiff.

Nota: Consulte la documentacion del paquete (- help(writeRaster) -) para ver argumentos utiles adicionales que se pueden agregar

TITULO: RELACION ENTRE BANDAS

Una matriz de diagrama de dispersion puede ser util para explorar las relaciones entre capas raster. Esto se puede hacer con la funcion pares () del paquete raster.

Trazado de reflejo en la longitud de onda ultra azul contra el reflejo en la longitud de onda azul.

pairs(landsatcrop2015[[1:2]], main = "Ultra-blue versus Blue")

Trazado de reflejo en la longitud de onda roja contra reflejo en la longitud de onda NIR.

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

La primera grafica revela altas correlaciones entre las regiones de longitud de onda azul. Debido a la alta correlacion, podemos usar una de las bandas azules sin perder mucha informacion.

Esta distribucion de puntos en la segunda grafica (entre NIR y rojo) es unica debido a su forma triangular. La vegetacion se refleja muy bien en el rango NIR que el rojo y crea la esquina superior cerca del eje NIR (y). El agua absorbe energia de todas las bandas y ocupa el lugar cercano al origen. El rincon mas alejado se crea debido a las caracteristicas superficiales altamente reflectantes, como el suelo brillante o el hormigon.

TITULO: EXTRAER VALORES DE PIXELES

A menudo, queremos obtener los valores de las celdas raster para ubicaciones geograficas o areas especificas. La - extract - funcion se utiliza para obtener valores raster en las ubicaciones de otros datos espaciales. Puede usar puntos, lineas, poligonos o un objeto de extension (rectangulo). Tambien puede usar numeros de celda para extraer valores. Al usar puntos, - extract - devuelve los valores de un - Raster* - objeto para las celdas en las que se ubica un conjunto de puntos.

# load the polygons with land use land cover information
sibatecobert = shapefile ('C:/Users/David/Desktop/PERCEPCION REMOTA/TRABAJO_COBERTURAS/cobertu.shp')

Z-dimension discarded

    # samp = readRDS('data/rs/samples.rds')              **en la guia esta asi**
# generate 300 point samples from the polygons
ptsibatecobert2015 = spsample(sibatecobert,300, type='regular')
    # ptsamp = spsample(samp, 300, type='regular')       **en la guia esta asi**
# add the land cover class to the points
ptsibatecobert$Name = over(ptsibatecobert, sibatecobert)$Name
    # ptsamp$class = over(ptsamp, samp)$class            **en la guia esta asi**
# extract values with points
df = extract(landsat2015, ptsibatecobert2015)
   # df = extract(landsat1, ptsamp)                      **en la guia esta asi**
# To see some of the reflectance values
head(df)

     ultra.blue blue green  red   NIR SWIR1 SWIR2
[1,]       8484 7800  7528 6931 14777 10879  8167
[2,]       8395 7782  7532 6774 14395 10388  7847
[3,]       8816 8265  8694 8181 21024 16753 11313
[4,]       8302 7651  7959 6534 25737 11668  7696
[5,]       8419 7791  8206 6847 24761 11801  7962
[6,]       8554 7902  8054 6912 21716 11099  7867

#     ultra.blue blue green  red   NIR SWIR1 SWIR2
#[1,]       8484 7800  7528 6931 14777 10879  8167
#[2,]       8816 8265  8694 8181 21024 16753 11313
#[3,]       8419 7791  8206 6847 24761 11801  7962
#[4,]       8399 7745  8087 6775 24775 12160  8012
#[5,]       8072 7394  6904 6273 11207  8098  6621
#[6,]       8072 7380  6827 6231 10728  8084  6654

TITULO: PERFILES ESPECTRALES

Una grafica del espectro (todas las bandas) para los pixeles que representan ciertas caracteristicas de la superficie terrestre (p. Ej. Agua) se conoce como perfil espectral. Dichos perfiles demuestran las diferencias en las propiedades espectrales de varias caracteristicas de la superficie terrestre y constituyen la base para el analisis de imagenes. Los valores espectrales se pueden extraer de cualquier conjunto de datos multiespectrales utilizando la - extract - funcion. En el ejemplo anterior, extrajimos valores de datos de Landsat para las muestras. Estas muestras incluyen: tierras de cultivo, agua, barbecho, construido y abierto. Primero calculamos los valores medios de reflectancia para cada clase y cada banda.

numero= length(ptsibatecobert)
numero

[1] 297

df_samples = as (ptsibatecobert, "SpatialPointsDataFrame")
df_samples

class       : SpatialPointsDataFrame 
features    : 297 
extent      : 577356.4, 585717, 493599.1, 501108.2  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :         Name 
min values  :         Agua 
max values  : zona desnuda 

df_samples@data=data.frame(ID=1:numero,size=1)
df_samples

class       : SpatialPointsDataFrame 
features    : 297 
extent      : 577356.4, 585717, 493599.1, 501108.2  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 2
names       :  ID, size 
min values  :   1,    1 
max values  : 297,    1 

plot(landsatcrop2015)
plot(sibatecobert, add= TRUE)

plot(df_samples,pch=1, cex=(df_samples$size)/4, add=TRUE)

df_samples$Name =over(df_samples,sibatecobert)$Name
df_samples

class       : SpatialPointsDataFrame 
features    : 297 
extent      : 577356.4, 585717, 493599.1, 501108.2  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       :  ID, size,         Name 
min values  :   1,    1,         Agua 
max values  : 297,    1, zona desnuda 

df1=raster::extract(landsatcrop2015, df_samples)
ms = aggregate(df1, list(ptsibatecobert$Name), mean)
## ms = aggregate(df, list(ptsamp$class), mean)         ** asi esta en la guia**
# 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.1864925 0.1795371 0.17953317 0.1958414 0.25448447 0.24850197
## cropland  0.1129813 0.0909645 0.08596722 0.0550344 0.48335462 0.16142085
## fallow    0.1319198 0.1164869 0.10453764 0.1151243 0.18012962 0.23139228
## open      0.1388014 0.1375235 0.15273163 0.2066425 0.34476670 0.35820877
## water     0.1336242 0.1165728 0.09922726 0.0785947 0.04909201 0.03360047
##               SWIR2
## built    0.20001306
## cropland 0.07314186
## fallow   0.19143030
## open     0.21346343
## water    0.02723398

Ahora trazamos el espectro medio de estas caracteristicas.

# Create a vector of color for the land cover classes for use in plotting
mycolor = c('darkred', 'yellow', 'burlywood', 'cyan', 'blue', 'green','pink')
#transform ms from a data.frame to a matrix
ms = as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,25000), 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")

El perfil espectral muestra (des) similitud en la reflectancia de diferentes caracteristicas en la superficie de la tierra (o por encima de ella). ‘Agua’ muestra una reflexion relativamente baja en todas las longitudes de onda, y ‘construido’, ‘en barbecho’ y ‘abierto’ tienen una reflectancia relativamente alta en las longitudes de onda mas largas.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.