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/Landsat_procesamiento"
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).

## solucion a la pregunta 1
landsatFCC2015_p1 = stack(b5, b6, b4)
plotRGB(landsatFCC2015_p1, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa - Agua / Tierra")

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

se da respuesta a esta pregunta con el ejercicio anterior

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

landsatcrop2015pregunta = crop(stack(b4, b3, b2), e)
plotRGB(landsatcrop2015pregunta, axes = TRUE, stretch = "lin", main = "Landsat Sibate 2015 - Color Verdadero")

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
ptsibatecobert = 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, ptsibatecobert)
   # 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,]       8816 8265  8694 8181 21024 16753 11313
[3,]       8577 7887  8185 7206 22881 12103  8138
[4,]       8399 7745  8087 6775 24775 12160  8012
[5,]       8204 7535  7134 6580 12286  9489  7410
[6,]       9797 9209  8825 8046 13806  9671  7765

#     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] 296

df_samples = as (ptsibatecobert, "SpatialPointsDataFrame")
df_samples

class       : SpatialPointsDataFrame 
features    : 296 
extent      : 577369.2, 585652.4, 493591.4, 501100.5  (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    : 296 
extent      : 577369.2, 585652.4, 493591.4, 501100.5  (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  : 296,    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    : 296 
extent      : 577369.2, 585652.4, 493591.4, 501100.5  (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  : 296,    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.

TITULO: OPERACIONES MATEMATICAS BASICAS

El - raster- paquete admite muchas operaciones matematicas. Las operaciones matematicas generalmente se realizan por pixel (celda de cuadricula). Primero haremos algunas operaciones aritmeticas basicas para combinar bandas. En el primer ejemplo, escribimos una funcion matematica personalizada para calcular el indice de vegetacion de diferencia normalizada (NDVI).

### solo puedo cargar las bandas de la 1 a la 7 porque la 8 no tiene la misma extension de las demas
raslistp3 = paste0('./data2015/LC08_L1TP_008057_20150104_20170415_01_T1_B', 1:7, ".tif")
##raslistp3= practica 3
landsatp3 = stack(raslistp3)
plot(landsatp3)

landsatRGBp3 = landsatp3[[c(4,3,2)]]
## landsatRGBp3 = practica 3
plot(landsatRGBp3)

landsatFCCp3 = landsatp3[[c(5,4,3)]]
## landsatFCCp3 = practica 3
plot(landsatFCCp3)

TITULO: INDICES DE VEGETACION

Definamos una funciOn general para un Indice basado en razOn (vegetaciOn). En la funcion a continuacion, -img- es un objeto Raster * de multiples capas -i- y -k- son los indices de las capas (numeros de capa) utilizados para calcular el indice de vegetacion.

## vi2015= indice de vegetacion para la imagen 2015
vi2015 = function(img, k, i) {
  
  bk = img [[k]]
  bi = img [[i]]
  vi2015=(bk-bi)/(bk+bi)
  return(vi2015)
}

## para landsat NIR = 5, red = 4
## ndvi2015= indice de vegetacion de direfencia normalizada (NDVI) para la imagen landsat 2015
ndvi2015 = vi2015(landsatp3, 5, 4) 
plot(ndvi2015, col=rev(terrain.colors(9)), main= "Landsat 2015 - NVDI")

NDVI1_sibate_2015 = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
eNDVI1_sibate = extent(NDVI1_sibate_2015)
# crop landsat by the extent
NDVI_recorte_1 = crop(ndvi2015, eNDVI1_sibate)
plot(NDVI_recorte_1, col=rev(terrain.colors(10)), main = "Landsat 2015 - NDVI")

Puedes ver la variacion en el verdor desde la trama.

Una forma alternativa de lograr esto es asi

vi2_2015 = function(x, y) {
  (x-y)/(x+y)
}
ndvi2_2015 = overlay(landsatp3[[5]], landsatp3[[4]], fun = vi2_2015)
plot(ndvi2_2015, col=rev(terrain.colors(10)), main = "Landsat 2015 - NDVI v2")

NDVI2_sibate_2015 = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
eNDVI2_sibate = extent(NDVI2_sibate_2015)
# crop landsat by the extent
NDVI_recorte_2 = crop(ndvi2_2015, eNDVI2_sibate)
plot(NDVI_recorte_2, col=rev(terrain.colors(10)), main = "Landsat 2015 - NDVI")

Pregunta 1 : Adapte el codigo que se muestra arriba para calcular los indices para identificar i) agua y ii) acumulados. Sugerencia: Use el diagrama del perfil espectral para encontrar las bandas que tienen reflectancia maxima y minima para estas dos clases.

## p1 es pregunta 1
ndvi2015_p1 = vi2015(landsatp3, 3, 6) 
plot(ndvi2015_p1, col=rev(terrain.colors(10)), main= "Landsat 2015 - NDWI")

vi2_2015_p1 = function(x, y) {
  (x-y)/(x+y)
}
ndvi2_2015_p1 = overlay(landsatp3[[3]], landsatp3[[6]], fun = vi2_2015_p1)
plot(ndvi2_2015_p1, col=rev(topo.colors(10)), main = "Landsat 2015 - NDWI")

### PARA INDICE DE AGUA DE DIFERENCIA NORMALIZADA - NDWI
NDWI_sibate = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
eNDWI = extent(NDWI_sibate)
# crop landsat by the extent
NDWI_recorte = crop(ndvi2_2015_p1, eNDWI)
plot(NDWI_recorte, col=rev(topo.colors(10)), main = "Landsat 2015 - NDWI")

### PARA INDICE INCORPORADO DE DIFERENCIA NORMALIZADA - NDBI
ndbi1_2015_p1 = overlay(landsatp3[[6]], landsatp3[[5]], fun = vi2_2015_p1)
plot(ndbi1_2015_p1, col=rev(heat.colors(10)), main = "Landsat 2015 - NDBI")

## recortado a sibate
NDBI_sibate = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
eNDBI = extent(NDBI_sibate)
# crop landsat by the extent
NDBI_recorte = crop(ndbi1_2015_p1, eNDBI)
plot(NDBI_recorte, col=rev(heat.colors(10)), main = "Landsat 2015 - NDWI")

TITULO: HISTOGRAMA

Podemos explorar la distribucion de valores contenidos en nuestro raster utilizando la funcion -hist ()- que produce un histograma. Los histogramas a menudo son utiles para identificar valores atipicos y valores de datos incorrectos en nuestros datos raster

## ver histograma de datos
hist(NDVI_recorte_1,
     main= "Dsitribucion de valores del NDVI",
     xlab= "NDVI",
     ylab= "Frecuencia",
     col= terrain.colors(50),
     xlim= c (-0.1,0.7),
     breaks = 30,
     xaxt= 'n')
axis(side = 1, at=seq(-0.5,1,0.05), labels = seq(-0.5,1,0.05))

Nos referiremos a este histograma para la siguiente subseccion sobre umbralizacion.

Pregunta 2 : Haga histogramas de los valores que los indices de vegetacion desarrollaron en la pregunta 1

## respuesta a la pregunta 2: ver histograma de datos
hist(NDWI_recorte,
     main= "Dsitribucion de valores del NDWI",
     xlab= "NDWI",
     ylab= "Frecuencia",
     col= topo.colors(50),
     xlim= c (-0.5,0.3),
     breaks = 30,
     xaxt= 'n')
axis(side = 1, at=seq(-0.5,1,0.05), labels = seq(-0.5,1,0.05))

## respuesta a la pregunta 2: ver histograma de datos
hist(NDBI_recorte,
     main= "Dsitribucion de valores del NDBI",
     xlab= "NDBI",
     ylab= "Frecuencia",
     col= heat.colors(50),
     xlim= c (-0.5,0.3),
     breaks = 30,
     xaxt= 'n')
axis(side = 1, at=seq(-0.5,1,0.05), labels = seq(-0.5,1,0.05))

TITULO: UMBRALIZACION

Podemos aplicar reglas basicas para obtener una estimacion de la extension espacial de diferentes caracteristicas de la superficie terrestre. Tenga en cuenta que los valores NDVI estan estandarizados y varian entre -1 y +1. Los valores mas altos indican mas cobertura verde.

Las celdas con valores de NDVI superiores a 0.4 son definitivamente vegetacion. La siguiente operacion enmascara todas las celdas que tal vez no sean vegetacion.

veg = reclassify(NDVI_recorte_1, cbind(-Inf, 0.4, NA))
plot(veg, main = "Vegetacion")

Vamos a mapear el area que corresponde al pico entre 0.25 y 0.3 en el histograma NDVI.

land = reclassify(NDVI_recorte_1, c(-Inf, 0.25, NA, 0.25, 0.3, 1, 0.3, Inf, NA))
plot(land, main = " Que es esto?")

Estas pueden ser las areas abiertas. Puede trazar -land- sobre el -landsatFCC- raster original para obtener mas informacion

plotRGB(landsatcrop2015, r=1, b=2, g=3, axes = TRUE, stretch="lin", main = "Landsat Sibate composicion en falso color")
plot(land, add = TRUE, legend= FALSE)

Tambien puede crear clases para diferentes cantidades de presencia de vegetacion.

vegc = reclassify(NDVI_recorte_1, 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 = 'Umbral basado en NDVI')

Pregunta 3 : Es posible encontrar agua usando el umbral de NDVI o cualquier otro indice

agua = reclassify(NDWI_recorte, cbind(-Inf, 0.1, NA))
plot(agua, main = "Agua", col= topo.colors(10))

aguac = reclassify(NDWI_recorte, c(-Inf,0.025,1, 0.025,0.03,2, 0.03,0.1,3, 0.1,0.2,4, 0.2,Inf, 5))
plot(aguac, col = rev(topo.colors(10)), main = 'Umbral basado en NDwI')

TITULO: ANALISIS DE COMPONENTES PRINCIPALES

Los datos multiespectrales a veces se transforman para ayudar a reducir la dimensionalidad y el ruido en los datos. La transformacion de componentes principales es un metodo generico de reduccion de datos que se puede utilizar para crear algunas bandas no correlacionadas a partir de un conjunto mas grande de bandas correlacionadas.

Puede calcular el mismo numero de componentes principales que el numero de bandas de entrada. El primer componente principal (PC) explica el mayor porcentaje de varianza y otras PC explican la varianza adicional en orden decreciente.

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

Esto se conoce como diagrama de vegetacion y linea de suelo (igual que el diagrama de dispersion en la seccion anterior).

pca = prcomp(sr, scale = TRUE)
pca

Standard deviations (1, .., p=7):
[1] 2.59320781 0.42283706 0.28432083 0.10272416 0.06739963 0.02013456 0.01197457

Rotation (n x k) = (7 x 7):
                                                   PC1         PC2         PC3         PC4         PC5
LC08_L1TP_008057_20150104_20170415_01_T1_B1 -0.3827161  0.13951992 -0.29577337  0.64426612 -0.09248285
LC08_L1TP_008057_20150104_20170415_01_T1_B2 -0.3823702  0.24607769 -0.25399056  0.24551600 -0.03813431
LC08_L1TP_008057_20150104_20170415_01_T1_B3 -0.3830553  0.23811360 -0.16711614 -0.23489380 -0.03847541
LC08_L1TP_008057_20150104_20170415_01_T1_B4 -0.3785436  0.42429526  0.00253034 -0.62149875 -0.06620595
LC08_L1TP_008057_20150104_20170415_01_T1_B5 -0.3624394 -0.75242358 -0.42394861 -0.24244188  0.23460165
LC08_L1TP_008057_20150104_20170415_01_T1_B6 -0.3767577 -0.34216544  0.52654673  0.04094435 -0.67771841
LC08_L1TP_008057_20150104_20170415_01_T1_B7 -0.3794551  0.00608041  0.60256852  0.15084020  0.68541151
                                                     PC6           PC7
LC08_L1TP_008057_20150104_20170415_01_T1_B1  0.070397514 -0.5639086914
LC08_L1TP_008057_20150104_20170415_01_T1_B2  0.262981944  0.7731983872
LC08_L1TP_008057_20150104_20170415_01_T1_B3 -0.842879576  0.0392600097
LC08_L1TP_008057_20150104_20170415_01_T1_B4  0.454485570 -0.2819057041
LC08_L1TP_008057_20150104_20170415_01_T1_B5  0.092049830 -0.0217900795
LC08_L1TP_008057_20150104_20170415_01_T1_B6 -0.007898927  0.0518065166
LC08_L1TP_008057_20150104_20170415_01_T1_B7 -0.018601359  0.0005871392

screeplot(pca)

pci = predict(landsatp3, pca, index = 1:2)

Warning in .Internal(gc(verbose, reset, full)) :
  closing unused connection 3 (C:/Users/David/AppData/Local/Temp/Rtmpw1HOgY/raster/r_tmp_2020-03-29_143139_208_24915.gri)

plot(pci[[1]])

## recortado a sibate
pci_sibate = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
epci = extent(pci_sibate)
# crop landsat by the extent
pci_recorte = crop(pci, epci)
plot(pci_recorte, col=rev(terrain.colors(20)), main = "pci Sibate")

El primer componente principal resalta los limites entre las clases de uso de la tierra o los detalles espaciales, que es la informacion mas comun entre todas las longitudes de onda. Es dificil entender lo que destaca el segundo componente principal. Vamos a intentar umbral de nuevo:

landsatFCC2015 = stack(b5, b4, b3)
# plotRGB(landsatFCC2015, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
## recortado a sibate
pc2 = shapefile('C:/Users/David/Desktop/SIBATE_RURAL/rural.shp')
epc2 = extent(pc2)
# crop landsat by the extent
pc2_rec= crop(landsatFCC2015, epc2)
# plot(pc2_rec, col=rev(terrain.colors(20)), main = "pci Sibate")
pc2 = reclassify(pci_recorte[[2]], c(-Inf,0,1,0,Inf,NA))
par(mfrow = c(1,2))
plotRGB(pc2_rec, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plotRGB(pc2_rec, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")

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

Para obtener mas informacion acerca de la informacion contenida en las parcelas de vegetacion y lineas de suelo, lea este documento de -Gitelson et al- . Una extension de PCA en la teledeteccion se conoce como Transformacion -Tasseled-cap- .

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.