Imagen faca LandSat 8

Landsat 8, La imagen corresponde al municipio de facatativa, departamento de cundinamarca, colombia, con frecha de 30 de diciembra del 2018

Seleccionar la carpeta de destino

# "C:\percepcionremota\PR\imagenfaca"
getwd()

[1] "C:/percepcionremota/PR/imagenfaca"

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

INTRODUCCION

####Para realizar procesamiento de imagenes satelitales encontramos diversas alternativas en la actualidad, ya sea con software libre o comercial, se afrima que opciones son altamente competitivas y poseen gran capacidad de procesamiento. Para este estudio se utilizaron R, Rstudio, que permiten realizar una clasificación de imagenes no supervisada y supervisada, ArcGIS para la elaboracion de shape, como otros procesos que se desarrollaran durante este estudio, explorando aspectos aspectos técnicos.

####Dentro del proyecto podemos resaltar algunas fases u operaciones como el preprocesado implementado a la imagen satelital donde se han aplicado una serie de transformaciones, la aplicación de técnicas de clasificación supervisada mediante la realización de entrenamiento y testeo, usando como area de interes el municipio de Facatativa.

IMPORTANTE: instalar paquetes: raster, sp, rgdal

Se implementó principalmente un subconjunto espacial de una escena Landsat 8 del municipipo de Facatativa, Cundinamarca, Colombia

TITULO: CAPITULO 1

TITULO: PROPIEDADES DE IMAGEN

Se crearonr objetos RasterLayer para capas individuales Landsat (bandas)

library(raster)
# Costera
b1 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B1.TIF')
# Azul
b2 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B2.TIF')
# Verde
b3 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B3.TIF')
# Rojo
b4 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B4.TIF')
# Infrarojo cercano (NIR)
b5 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B5.TIF')
# Infrarrojo de Onda Corta 1 (SWIR 1)
b6 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B6.tif')
# Infrarrojo de Onda Corta 2 (SWIR 2)
b7 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B7.tif')
# Pancromatica
b8 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B8.tif')
# Cirros (Cirrus)
b9 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B9.tif')
# Sensor Térmico Infrarrojo 1 (TIRS 1)
b10 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B10.tif')
# Sensor Térmico Infrarrojo 2 (TIRS 2)
b11 = raster('./LC08_L1TP_008057_20181230_20190130_01_T1_B11.tif')

A continuación podemos imprimir las variables a verificar, por ejemplo con la banda dos (b2):

b2

class      : RasterLayer 
dimensions : 7741, 7581, 58684521  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 446685, 674115, 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:/percepcionremota/PR/imagenfaca/LC08_L1TP_008057_20181230_20190130_01_T1_B2.TIF 
names      : LC08_L1TP_008057_20181230_20190130_01_T1_B2 
values     : 0, 58849  (min, max)

Podemos observar 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 presentan las diferentes propiedades desde un objeto Raster

# Sistema de coordenadas de referencia (CRS)
crs(b2)

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

# Numero de celdas, filas, columnas
ncell(b2)

[1] 58684521

# Dimension
dim(b2)

[1] 7741 7581    1

# Resolucion espacial
res(b2)

[1] 30 30

# Numero de badas 
nlayers(b2)

[1] 1

# ¿Las bandas tienen la misma extension, numero de filas y columnas, proyeccion, resolucion y origen?
compareRaster(b2,b3)

[1] TRUE

## [1] TRUE

Puede crear un RasterStack con la misma extensión y resolución espacial. Se puede crear un RasterStack a partir de objetos RasterLayer, o de archivos ráster, o ambos.

# Creamos el RasterStack
s2015 = stack(b5, b4, b3)

# Revisamos las propiedades del RasterStack
s2015

class      : RasterStack 
dimensions : 7741, 7581, 58684521, 3  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 446685, 674115, 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_20181230_20190130_01_T1_B5, LC08_L1TP_008057_20181230_20190130_01_T1_B4, LC08_L1TP_008057_20181230_20190130_01_T1_B3 
min values :                                           0,                                           0,                                           0 
max values :                                       65535,                                       62063,                                       58803 

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

# Primero creamos una lista de los raster layers usando:

filenames = paste0('./LC08_L1TP_008057_20181230_20190130_01_T1_B', 1:7, ".tif")
filenames

[1] "./LC08_L1TP_008057_20181230_20190130_01_T1_B1.tif" "./LC08_L1TP_008057_20181230_20190130_01_T1_B2.tif"
[3] "./LC08_L1TP_008057_20181230_20190130_01_T1_B3.tif" "./LC08_L1TP_008057_20181230_20190130_01_T1_B4.tif"
[5] "./LC08_L1TP_008057_20181230_20190130_01_T1_B5.tif" "./LC08_L1TP_008057_20181230_20190130_01_T1_B6.tif"
[7] "./LC08_L1TP_008057_20181230_20190130_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, la banda 9 es cirro, y las bandas 10 y 11 son termales

## [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"

# Creo el RasterStack de las bandas seleccionadas, es decir  de la 1 a la  7 (1:7)
landsatFaca = stack (filenames)
landsatFaca

class      : RasterStack 
dimensions : 7741, 7581, 58684521, 7  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 446685, 674115, 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//0_01_T1_B1, LC08_L1TP//0_01_T1_B2, LC08_L1TP//0_01_T1_B3, LC08_L1TP//0_01_T1_B4, LC08_L1TP//0_01_T1_B5, LC08_L1TP//0_01_T1_B6, LC08_L1TP//0_01_T1_B7 
min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
max values :                 56188,                 58849,                 58803,                 62063,                 65535,                 45351,                 37289 

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

# Trazando capas indivuduales (aZUL, VERDE, ROJO, INFRAROJO CERCANO) a partir del RasterStack
par(mfrow = c(2,2))
plot(b2, main = "Azul", col = gray(0:100 / 100), axes = FALSE)
plot(b3, main = "Verde", col = gray(0:100 / 100), axes = FALSE)

plot(b4, main = "Rojo", col = gray(0:100 / 100), axes = FALSE)
plot(b5, main = "NIR", col = gray(0:100 / 100), axes = FALSE)

En la imagen anterior podemos observar que ls leyendas Pueden variar entre 0 y 1. las diferentes caracteristicas de la superficie reflejan la radiacion solar incidente de manera diferente, esto se observa en la imagen en la diferencial del tono de grises, 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.

#####Cuando las bandas espectrales de un determinado sensor coinciden exactamente con las longitudes de onda de los colores Rojo, Verde y Azul, la composición resultante se conoce como “color verdadero”. Las composiciones RGB se utilizan para analizar y resaltar de forma visual determinada información de la imagen, como ,o presenta la siguiente imagen:

# Ploteamos la combinacion de bandas 4, 3 y 2 que para L8 es color verdadero
landsatRGBFaca <- stack(b4, b3, b2)
plotRGB(landsatRGBFaca, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

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 metodo - plotRGB - 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” -

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 como “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).

# Ploteamos la combinacion de bandas 5, 4 y 3 que para L8 es color falso
par(mfrow = c(1,2))
plotRGB(landsatRGBFaca, 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 : 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, ploteamos la combinacion de bandas 5, 6 y 4 que para L8 es color falso
landsatfaca_p1 = stack(b5, b6, b4)
plotRGB(landsatfaca_p1, axes=TRUE, stretch="lin", main="Composicion de color falsa - Agua / Tierra")

TITULO: SUBCONJUNTO Y RENOMBRAR BANDAS

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

# select first 3 bands only
landsatFacasub1 <- subset(landsatFaca, 1:3)
# same
landsatFacasub2 <- landsatFaca[[1:3]]
# Number of bands in the original and new data
nlayers(landsatFaca)

[1] 7

nlayers(landsatFacasub1)

[1] 3

nlayers(landsatFacasub2)

[1] 3

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

## Eliminando las ultimas cuatro bandas, es decir, seleccionaremos para nuestro ejercicio las bandas de 1 a la 7 (1:7)
landsatFaca <- subset(landsatFaca, 1:7)

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

# Aca vemos que nombres tienen las bandas
names(landsatFaca) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsatFaca)

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

plot (landsatFaca)

## [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"

# Aca las renombramos y las ploteamos
## [1] "ultra.azuk" "azul"       "verde"      "rojo"       "NIR"        "SWIR1"      "SWIR2"     
## [1] "ultra.azul" "azul"       "verde"      "rojo"       "NIR"        "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 funcion - crop -, utilizando un objeto - extent - u otro objeto espacial del que se puede extraer una Extension.

# Usando Extension
extent(landsatFaca)

class      : Extent 
xmin       : 446685 
xmax       : 674115 
ymin       : 362985 
ymax       : 595215 

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

# Cargamos un objeto espacial que se encuentra en un shapefile. Esta extension es el area de faca
Facashape = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
e = extent(landsatFaca)



###crop landsat 

landsatFacacrop = crop (landsatFaca, e)

plotRGB(landsatFacacrop, axes = TRUE, stretch= "lin", main= "Landsat composicion a falso color faca")

NA
NA

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

# Creacion de un compuesto de color verdadero, para landsat 8 la combinacion de color es 4,3,2
landsatcropfacasolucion = crop(stack(b4, b3, b2), e)

# Ploteamos la combinacion
plotRGB(landsatcropfacasolucion, axes = TRUE, stretch = "lin", main = "imagen faca - 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

# Aca guardanmos nuestro raster de faca y lo imprimimos para ver lo que se ha guardado
x = writeRaster(landsatFaca, filename="faca_banda.tif", overwrite=TRUE)
plot(x)

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

# Lo guardamos tambien en formato .grd ya que este formato si guarda los nombres de las capas
writeRaster(landsatFaca, filename="faca_banda.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.

# Comparacion del reflejo en la longitud de onda Ulta azul vs azul
pairs(landsatFacacrop [[1:2]], main = "ultra=blue vs blue" )

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

# Comparacion del reflejo en la longitud de onda Roja vs NIR
pairs(landsatFacacrop[[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 en 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 funcion - extract - 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 objeto - Raster* - para las celdas en las que se ubica un conjunto de puntos.

# Cargar los poligonos con informacion sobre el uso del suelo
facacobert = shapefile ('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
    
# Generar muestras de 300 puntos a partir de los poligonos
ptsfacacobert = spsample(facacobert,300, type='regular')

# Agregar la clase de cobertura del suelo a los puntos
ptsfacacobert$Name = over(ptsfacacobert, facacobert)$Name
    
# Extraer valores con puntos
df = extract(landsatFaca, ptsfacacobert)
   
# Para ver algunos de los valores de reflectancia
head(df)

     ultra.blue  blue green   red   NIR SWIR1 SWIR2
[1,]       8937  8344  8179  7794 13903 12615  9427
[2,]       9169  8481  8421  7857 15754 11541  8527
[3,]      10319 10032  9834  9328 14415 12669 10064
[4,]      10912 10810 10895 11307 15852 15358 12324
[5,]       8739  8084  7965  7107 18447 10927  7865
[6,]       8905  8264  8213  7922 14526 12878  9336

#     ultra.azul azul verde rojo  NIR  SWIR1 SWIR2
#[1,]       7928 7238  6623 6033 12696  7520  5997
#[2,]       7957 7278  6700 6101 13199  7933  6123
#[3,]       7973 7285  6741 6148 13482  8072  6219
#[4,]       7962 7262  6710 6114 13413  8116  6199
#[5,]       7931 7246  6671 6060 13541  7655  6014
#[6,]       7943 7242  6623 6024 12437  7529  6039

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 funcion - extract -. 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(ptsfacacobert)
numero

[1] 306

df_samples = as (ptsfacacobert, "SpatialPointsDataFrame")
df_samples

class       : SpatialPointsDataFrame 
features    : 306 
extent      : 564553.6, 581280.3, 526517.1, 541062.1  (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  :          artificial 
max values  : plantacion_forestal 

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

class       : SpatialPointsDataFrame 
features    : 306 
extent      : 564553.6, 581280.3, 526517.1, 541062.1  (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  : 306,    1 

plot(landsatFacacrop)
plot(facacobert, add= TRUE)

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

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

class       : SpatialPointsDataFrame 
features    : 306 
extent      : 564553.6, 581280.3, 526517.1, 541062.1  (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,          artificial 
max values  : 306,    1, plantacion_forestal 

df1=raster::extract(landsatFacacrop, df_samples)


ms = aggregate(df1, list(ptsfacacobert$Name), mean)

# En lugar de la primera columna, usamos nombres de fila
rownames(ms) = ms[,1]
ms = ms[,-1]
ms

##              ultra.blue       blue      green       red        NIR      SWIR1    SWIR2
# Agua           8251.053     7501.298   6840.658     6066.404  5762.921    5081.289    5050.825
# Bosque           8026.761   7326.134   6773.410     6172.284  11207.410   7962.134    6412.746
# Cultivos       8417.864     7800.318   7603.682     7197.955  18156.818   14063.045   11198.864
# Fresa          8541.667   7870.667     7464.667     6902.000  13157.667   9832.333    7806.000
# Invernadero     11522.083  11035.833  11442.333    11179.083  18368.333   14785.000   11300.083
# pastos           8481.000   7926.727   8267.909     7331.636  21093.273   12449.000   8507.273
# zona desnuda  10975.500    11035.500  12884.750    14428.000  17930.750   20671.750   16907.500

Ahora trazamos el espectro medio de estas caracteristicas.

# Cree un vector de color para las clases de cobertura del suelo para usar en el trazado
mycolor = c('darkred', 'yellow', 'burlywood', 'cyan', 'blue', 'green','magenta', "darkgreen")

# transformar ms de un data.frame a una matriz
ms = as.matrix(ms)

# Primero crea una parcela vacia
plot(0, ylim=c(0,25000), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Reflectancia")

# Agrega las diferentes clases
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

# Titulo
title(main="Perfil espectral de Landsat", font.main = 2)

# Leyenda
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 ‘Pastos’, ‘Invernadero’ y ‘Zona desnuda’ tienen una reflectancia relativamente alta en las longitudes de onda mas largas.

TITULO: OPERACIONES MATEMATICAS BASICAS

El paquete - raster- 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('./LC08_L1TP_008057_20181230_20190130_01_T1_B', 1:7, ".tif")
##raslistp3= practica 3

# Stack de las bandas cargadas
landsatp3 = stack(raslistp3)
plot(landsatp3)

# Combinacion de las bandas 4,3,2
landsatRGBp3 = landsatp3[[c(4,3,2)]]
## landsatRGBp3 = practica 3
plot(landsatRGBp3)

# Combinacion de las bandas en falso color 5,4,3
landsatFCCp3 = landsatp3[[c(5,4,3)]]
## landsatFCCp3 = practica 3
plot(landsatFCCp3)

NA
NA

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.

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

## para landsat NIR = 5, red = 4
## ndvifaca= indice de vegetacion de direfencia normalizada (NDVI) para la imagen landsat 2015

ndvifaca = vifaca(landsatp3, 5, 4) 
plot(ndvifaca, col=rev(terrain.colors(9)), main= "Landsat faca - NVDI")

# Seleccionando el shape del area de faca, cargandolo y recortandolo con respecto del NDVI de la imagen Landsat
NDVI1_faca = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18')
eNDVI1_faca = extent(NDVI1_faca)
NDVI_recorte_1 = crop(NDVI1_faca, eNDVI1_faca)

Puedes ver la variacion en el verdor desde la trama.

Una forma alternativa de lograr esto es asi

vi2_faca = function(x, y) {
  (x-y)/(x+y)
}

ndvi2_faca = overlay(landsatp3[[5]], landsatp3[[4]], fun = vi2_faca)
plot(ndvi2_faca, col=rev(terrain.colors(10)), main = "Landsat faca5 - NDVI v2")

# Seleccionando el shape del area de faca, cargandolo y recortandolo con respecto del NDVI de la imagen Landsat
NDVI2_faca = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
eNDVI2_faca = extent(NDVI2_faca)
NDVI_recorte_2 = crop(NDVI2_faca, eNDVI2_faca)
plot(NDVI_recorte_2, col=rev(terrain.colors(10)), main = "Landsat faca - NDVI")

Pregunta 1 : Se utiliza el codigo anterioemente relacionado y se ajusta a la necesidad del proyecto para de esta forma calcular los indices y 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.

## s1 es solucion 1

ndvifaca_s1 = vifaca(landsatp3, 3, 6) 
plot(ndvifaca_s1, col=rev(terrain.colors(10)), main= "Landsat faca - NDWI")

vi2_faca_s1 = function(x, y) {
  (x-y)/(x+y)
}

ndvi2_faca_s1 = overlay(landsatp3[[3]], landsatp3[[6]], fun = vi2_faca_s1)
plot(ndvi2_faca_s1, col=rev(topo.colors(10)), main = "Landsat faca - NDWI")

### PARA INDICE DE AGUA DE DIFERENCIA NORMALIZADA - NDWI
# Seleccionando el shape del area de faca, cargandolo y recortandolo con respecto del NDwI de la imagen Landsat

NDWI_faca = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
eNDWI = extent(NDWI_faca)
NDWI_recorte = crop(ndvi2_faca_s1, eNDWI)
plot(NDWI_recorte, col=rev(topo.colors(10)), main = "Landsat faca - NDWI")

### PARA INDICE INCORPORADO DE DIFERENCIA NORMALIZADA - NDBI

ndbi1_faca_s1 = overlay(landsatp3[[6]], landsatp3[[5]], fun = vi2_faca_s1)
plot(ndbi1_faca_s1, col=rev(heat.colors(10)), main = "Landsat faca - NDBI")

# Seleccionando el shape del area de faca, cargandolo y recortandolo con respecto del NDBI de la imagen Landsat

NDBI_faca = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
eNDBI = extent(NDBI_faca)
NDBI_recorte = crop(ndbi1_faca_s1, eNDBI)
plot(NDBI_recorte, col=rev(heat.colors(10)), main = "Landsat faca - NDBI")

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(NDBI_recorte,
     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(NDBI_recorte, cbind(-Inf, 0.4, NA))

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

plot(veg, main = "Vegetacion")

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

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

land = reclassify(NDBI_recorte, 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(landsatFaca, r=7, g=6, b=4, axes = TRUE, stretch="lin", main = "faca composicion en falso color")
plot(land, add = TRUE, legend= FALSE)

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

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

vegc = reclassify(NDBI_recorte, 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

# Umbral, en este caso para agua
agua = reclassify(NDBI_recorte, cbind(-Inf, 0.1, NA))
plot(agua, main = "Agua", col= topo.colors(10))

#Reclasifico el Umbral, en este caso para agua
aguac = reclassify(NDBI_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')

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 otros PC explican la varianza adicional en orden decreciente.

set.seed(1)

sr = sampleRandom(landsatp3, 100000)
plot(sr[,c(4,5)], main = "NIR - ROJO")

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.59003689 0.44921866 0.27644178 0.09269922 0.06591655 0.02057198 0.01140578

Rotation (n x k) = (7 x 7):
                                                   PC1        PC2         PC3         PC4         PC5
LC08_L1TP_008057_20181230_20190130_01_T1_B1 -0.3828966  0.2384366 -0.13610614 -0.60007235 -0.32922500
LC08_L1TP_008057_20181230_20190130_01_T1_B2 -0.3817989  0.3228502 -0.06913742 -0.24378800 -0.12020857
LC08_L1TP_008057_20181230_20190130_01_T1_B3 -0.3821098  0.3131713 -0.02221730  0.21450895  0.09840963
LC08_L1TP_008057_20181230_20190130_01_T1_B4 -0.3789077  0.4004274  0.13575699  0.58234331  0.17008071
LC08_L1TP_008057_20181230_20190130_01_T1_B5 -0.3695465 -0.4480919 -0.74943217  0.08558584  0.29581948
LC08_L1TP_008057_20181230_20190130_01_T1_B6 -0.3741646 -0.5129561  0.25823666  0.29187454 -0.66453292
LC08_L1TP_008057_20181230_20190130_01_T1_B7 -0.3761344 -0.3414484  0.57396587 -0.32066662  0.55626948
                                                     PC6         PC7
LC08_L1TP_008057_20181230_20190130_01_T1_B1  0.127061460 -0.54165739
LC08_L1TP_008057_20181230_20190130_01_T1_B2 -0.488236233  0.65799701
LC08_L1TP_008057_20181230_20190130_01_T1_B3  0.781106318  0.29932697
LC08_L1TP_008057_20181230_20190130_01_T1_B4 -0.361200728 -0.42324897
LC08_L1TP_008057_20181230_20190130_01_T1_B5 -0.068502312 -0.03838947
LC08_L1TP_008057_20181230_20190130_01_T1_B6  0.011655828  0.05709746
LC08_L1TP_008057_20181230_20190130_01_T1_B7 -0.007700369 -0.01330462

screeplot(pca)

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

## Recortado a faca
pci_faca = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
epci = extent(pci_faca)
pci_rec = crop(pci, epci)
plot(pci_rec, col=rev(terrain.colors(20)), main = "pci faca")

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:

landsatfcfaca = stack(b5, b4, b3)

## recortado a faca
pc2 = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
epc2 = extent(pc2)
pc2_rec= crop(landsatfcfaca, epc2)

pc2 = reclassify(pci_rec  [[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 = "FFC")
plotRGB(pc2_rec, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "FFC")

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

TITULO: CAPITULO 2

En este capitulo exploramos la clasificacion no supervisada. Existen varios algoritmos de clasificación no supervisados, y la eleccion del algoritmo puede afectar los resultados. Exploraremos un solo algoritmo (k-means) para ilustrar el principio general

faca, Imagen LandSat 8, utilizando 7 bandas

library(raster)
landsat8 = stack(b2,b3,b4,b5,b6,b7)
names(landsat8) = c ('blue', 'green', 'red', 'NIR', 'SWIR1','SWIR2')
plot(landsat8)

Pregunta 1 : Haga una trama compuesta de 3 bandas de falso color de `` landsat8 ’’.

landsat8faca_imagen = stack(b5, b6, b4)
plotRGB(landsat8faca_imagen, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa - Agua / Tierra")

# recortando a faca

facaextend = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
ext = extent(facaextend)
landsatfcc = crop(landsat8faca_imagen, ext)
plotRGB(landsatfcc, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

En la clasificacion no supervisada, utilizamos los datos de reflectancia, pero no proporcionamos ningun dato de respuesta (es decir, no identificamos ningun pixel como perteneciente a una clase en particular). Esto puede parecer extraño, pero puede ser util cuando no tenemos mucho conocimiento previo de un area de estudio. O si tiene un amplio conocimiento de la distribucion de las clases de interes de la cobertura del suelo, pero no tiene datos especificos

El algoritmo agrupa pixeles con caracteristicas espectrales similares en grupos.

Obtenga mas informacion sobre K-means y otros algoritmos supervisados sin supervision .

Realizaremos una clasificacion no supervisada en un subconjunto espacial de la capa -ndvi-. Aqui hay otra forma de calcular -ndvi-. En este caso, no usamos una funcion separada, sino una notacion algebraica directa.

ndviLandSat8=(landsat8[['NIR']]-landsat8[['red']])/(landsat8[['NIR']]+landsat8[['red']])
plot(ndviLandSat8)

Haremos un agrupamiento -kmeans- de los datos -ndvi-. Primero usamos -crop- para hacer un subconjunto espacial de -ndvi-, para permitir un procesamiento mas rapido (puede seleccionar cualquier -extent- usando la función -drawExtent()).

TITULO: CLASIFICACION KMEANS

landsat8faca_imagen = stack(b5, b6, b4)
plotRGB(landsat8faca_imagen, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa - Agua / Tierra")

# recortando a faca

facaextend = shapefile('C:/percepcionremota/ShpFaca/coberturatierrafaca18.shp')
ext = extent(facaextend)
landsatfcc = crop(landsat8faca_imagen, ext)
plotRGB(landsatfcc, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

NA
NA
NA

Tenga en cuenta que -getValues- convirtio el RasterLayer -ndvi- en una matriz (matriz). Ahora realizaremos el agrupamiento -kmeans- en la matriz e inspeccionaremos la salida.

ndviLandSat8=(landsat8[['NIR']]-landsat8[['red']])/(landsat8[['NIR']]+landsat8[['red']])
plot(ndviLandSat8)

-kmeans-devuelve un objeto con 9 elementos. La longitud del elemento -cluster- dentro -kmncluster- es 216567 que es igual a la longitud de -nr- creado desde -ndvi-. Los valores de celda de rango -kmncluster\(cluster- entre 1 a 10 correspondientes al numero de entrada del cluster que proporcionamos en la funcion -kmeans-. -kmncluster\)cluster- indica la etiqueta del cluster para el pixel correspondiente. Necesitamos convertir los valores -kmncluster$cluster- nuevamente a RasterLayer de la misma dimension que -ndvi-.

# Es importante configurar el generador de semillas porque `kmeans` inicia los centros en ubicaciones aleatorias, para eso utilizaremos set.seed
set.seed(99)

# Queremos crear 10 grupos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios utilizando el método "Lloyd"
kmncluster = kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")

Error in na.omit(nr) : objeto 'nr' no encontrado

Podemos ver que -knr- es un RasterLayer, pero no sabemos que grupo (1-10) pertenece a que clase de cobertura terrestre (y si pertenece a una clase que reconoceriamos). Puede averiguarlo trazandolos lado a lado con capas de referencia y usando un color unico para cada grupo.

# Use el objeto ndvi (en este caso se llama "landsat8cropfaca") para establecer los valores del cluster en un nuevo raster
knr = setValues(landsatFacacrop, kmncluster$cluster)

# Tambien puedes hacerlo asi
knr= raster(landsatFacacrop)
values(knr) = kmncluster$cluster
knr

Si bien para otros fines suele ser mejor definir mas clases (y posiblemente fusionar clases mas adelante), una clasificacion simple como esta podria ser util, por ejemplo, fusionar los grupos 4 y 5 para construir una mascara de agua para el año 2015.

Puedes cambiar los colores en mi -mycolor-. Obtenga mas informacion sobre la seleccion de colores en R -aqui- y -aqui- .

Pregunta 2 : Grafique el RGB de 3 bandas de “landsat8” para el subconjunto (extensión “ext”) y el resultado de la agrupacion “kmeans” lado a lado y haga una tabla de cobertura del suelo para el uso del suelo Etiquetas para los grupos. Por ejemplo, los grupos 4 y 5 son agua.

landsat8p2 = stack(b6, b3, b2)
#plotRGB(landsat8p2, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa ")

#crop landsat by the extent
landsatfacap2 = crop(landsat8p2, ext)
plotRGB(landsatfacap2, axes = TRUE, stretch = "lin", main = "Landsat FCC - faca")

# convertir el raster a vector / matriz
nr_p2=getValues(landsatfacap2)
str(nr_p2)

 int [1:300979, 1:3] 8933 9002 10072 10961 11390 12368 12708 13320 13510 13658 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3] "LC08_L1TP_008057_20181230_20190130_01_T1_B6" "LC08_L1TP_008057_20181230_20190130_01_T1_B3" "LC08_L1TP_008057_20181230_20190130_01_T1_B2"

# Es importante configurar el generador de semillas porque `kmeans` inicia los centros en ubicaciones aleatorias, para eso utilizaremos set.seed
set.seed (99)

# Queremos crear 10 grupos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios utilizando el método "Lloyd"
kmncluster_p2 = kmeans(na.omit(nr_p2), centers = 10, iter.max = 500, nstart = 5, algorithm = "Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"
str(kmncluster_p2)

List of 9
 $ cluster     : int [1:300979] 9 9 9 6 6 8 8 8 8 8 ...
 $ centers     : num [1:10, 1:3] 14704 14166 26021 19680 33698 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : chr [1:3] "LC08_L1TP_008057_20181230_20190130_01_T1_B6" "LC08_L1TP_008057_20181230_20190130_01_T1_B3" "LC08_L1TP_008057_20181230_20190130_01_T1_B2"
 $ totss       : num 1.18e+13
 $ withinss    : num [1:10] 3.76e+10 9.01e+10 9.30e+10 9.32e+10 1.04e+11 ...
 $ tot.withinss: num 6e+11
 $ betweenss   : num 1.12e+13
 $ size        : int [1:10] 40372 16656 5230 6859 3016 64330 14028 67531 49945 33012
 $ iter        : int 115
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

# Use el objeto ndvi (en este caso se llama "landsatfacap2") para establecer los valores del cluster en un nuevo raster
knr_p2 = setValues(landsatfacap2, kmncluster_p2$cluster)

# Tambien puedes hacerlo asi
knr_p2= raster(landsatfacap2)
values(knr_p2) = kmncluster_p2$cluster
knr_p2

class      : RasterLayer 
dimensions : 511, 589, 300979  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 564225, 581895, 526425, 541755  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 1, 10  (min, max)

#Defina un vector de color para 10 grupos (mas informacion sobre como configurar el color mas adelante)
#mycolor_p2 = c (topo.colors(10))
plot(knr_p2, main = 'Unsupervised classification', col = rev(topo.colors(10))) 

TITULO: CAPITULO 3

TITULO: CLASIFICACION SUPERVISADA

Aqui exploramos la clasificacion supervisada. Existen varios algoritmos de clasificacion supervisados, y la eleccion del algoritmo puede afectar los resultados. Aqui exploramos dos algoritmos relacionados (CART y RandomForest).

En la clasificacion supervisada, tenemos conocimiento previo sobre algunos de los tipos de cobertura del suelo mediante, por ejemplo, trabajo de campo, datos espaciales de referencia o interpretación de imagenes de alta resolucion (como las disponibles en los mapas de Google). Se identifican sitios especificos en el area de estudio que representan ejemplos homogeneos de estos tipos de cobertura terrestre conocidos. Estas areas se conocen comunmente como sitios de entrenamiento porque las propiedades espectrales de estos sitios se usan para entrenar el algoritmo de clasificacion.

Los siguientes ejemplos utilizan un clasificador de arboles de clasificacion y regresion (CART) (Breiman et al. 1984) ( lecturas adicionales para predecir las clases de cobertura del suelo en el area de estudio).

Realizaremos los siguientes pasos:

Genere sitios de muestra basados en un raster de referencia

Extraer valores de celda de datos de Landsat para los sitios de muestra

Entrenar al clasificador usando muestras de entrenamiento

Clasifique los datos de Landsat usando el modelo entrenado

Evaluar la precision del modelo.

TITULO: DATOS DE REFERENCIA

La National Land Cover Database 2011 (NLCD 2011) es un producto de cobertura terrestre para los EE. UU. NLCD es una base de datos de cobertura del suelo basada en Landsat de 30 m que abarca 4 épocas (1992, 2001, 2006 y 2011). NLCD 2011 se basa principalmente en una clasificación de árbol de decisión de datos Landsat de alrededor de 2011.

Puede encontrar los nombres de clase en NCLD 2011 (aquí) -[ https://www.mrlc.gov/nlcd11_leg.php ]-. Tiene dos pares de valores de clase y nombres que corresponden a los niveles de uso del suelo y el sistema de clasificacion de la cobertura del suelo. Estos niveles generalmente representan el nivel de complejidad, siendo el nivel I el mas simple con amplias categorias de cobertura del suelo. Lea este informe de Anderson et al. Para obtener mas informacion sobre este sistema de clasificacion de uso y cobertura de la tierra.

CAPITULO 3

CLASIFICACION SUPERVISADA

Aqui exploramos la clasificacion supervisada. Existen varios algoritmos de clasificacion supervisados, y la eleccion del algoritmo puede afectar los resultados. Aqui exploramos dos algoritmos relacionados (RandomForest y arbol de decision).

En la clasificacion supervisada, tenemos conocimiento previo sobre algunos de los tipos de cobertura del suelo mediante, por ejemplo, trabajo de campo, datos espaciales de referencia o interpretación de imagenes de alta resolucion (como las disponibles en los mapas de Google). Se identifican sitios especificos en el area de estudio que representan ejemplos homogeneos de estos tipos de cobertura terrestre conocidos. Estas areas se conocen comunmente como sitios de entrenamiento porque las propiedades espectrales de estos sitios se usan para entrenar el algoritmo de clasificacion.

# Directorio de trabajo

setwd("C:/percepcionremota/PR/imagenfaca/")

# Importar las imagenes
B2 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B2.TIF")
B3 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B3.TIF")
B4 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B4.TIF")
B5 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B5.TIF")
B6 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B6.TIF")
B7 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B7.TIF")

# Datos imagen MTL
dia <- 004
SUN_ELEVATION = 51.89134539

RADIANCE_MULT_BAND_2 = 1.3298E-02
RADIANCE_MULT_BAND_3 = 1.2254E-02
RADIANCE_MULT_BAND_4 = 1.0333E-02
RADIANCE_MULT_BAND_5 = 6.3236E-03
RADIANCE_MULT_BAND_6 = 1.5726E-03
RADIANCE_MULT_BAND_7 = 5.3006E-04

RADIANCE_ADD_BAND_2 = -66.49141
RADIANCE_ADD_BAND_3 = -61.27126
RADIANCE_ADD_BAND_4 = -51.66738
RADIANCE_ADD_BAND_5 = -31.61786
RADIANCE_ADD_BAND_6 = -7.86307
RADIANCE_ADD_BAND_7 = -2.65028

REFLECTANCE_MAXIMUM_BAND_2 = 1.210700
REFLECTANCE_MAXIMUM_BAND_3 = 1.210700
REFLECTANCE_MAXIMUM_BAND_4 = 1.210700
REFLECTANCE_MAXIMUM_BAND_5 = 1.210700
REFLECTANCE_MAXIMUM_BAND_6 = 1.210700
REFLECTANCE_MAXIMUM_BAND_7 = 1.210700

RADIANCE_MAXIMUM_BAND_2 = 805.01141
RADIANCE_MAXIMUM_BAND_3 = 741.81116
RADIANCE_MAXIMUM_BAND_4 = 625.53699
RADIANCE_MAXIMUM_BAND_5 = 382.79742
RADIANCE_MAXIMUM_BAND_6 = 95.19823
RADIANCE_MAXIMUM_BAND_7 = 32.08690

# Conversion
SUN_ELEVATION_R <- SUN_ELEVATION*pi/180  # radianes
d <- 1+0.0167*(sin((2*pi*(dia-93.5))/365)) # distancia del sol a la tierra

# Determinacion ESUN
ESUN2 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_2/REFLECTANCE_MAXIMUM_BAND_2
ESUN3 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_3/REFLECTANCE_MAXIMUM_BAND_3
ESUN4 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_4/REFLECTANCE_MAXIMUM_BAND_4
ESUN5 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_5/REFLECTANCE_MAXIMUM_BAND_5
ESUN6 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_6/REFLECTANCE_MAXIMUM_BAND_6
ESUN7 <- pi*(d^2)*RADIANCE_MAXIMUM_BAND_7/REFLECTANCE_MAXIMUM_BAND_7

# Elimar valores nulos
B2[B2==0] <- NA
B3[B3==0] <- NA
B4[B4==0] <- NA
B5[B5==0] <- NA
B6[B6==0] <- NA
B7[B7==0] <- NA

# Determinacion de radiancia sensor
L2 <- RADIANCE_MULT_BAND_2*B2 + RADIANCE_ADD_BAND_2
L3 <- RADIANCE_MULT_BAND_3*B3 + RADIANCE_ADD_BAND_3
L4 <- RADIANCE_MULT_BAND_4*B4 + RADIANCE_ADD_BAND_4
L5 <- RADIANCE_MULT_BAND_5*B5 + RADIANCE_ADD_BAND_5
L6 <- RADIANCE_MULT_BAND_6*B6 + RADIANCE_ADD_BAND_6
L7 <- RADIANCE_MULT_BAND_7*B7 + RADIANCE_ADD_BAND_7

# Determinacion de la radiancia minima
Lmin2 <- RADIANCE_MULT_BAND_2*min(getValues(B2), na.rm = TRUE)+ RADIANCE_ADD_BAND_2
Lmin3 <- RADIANCE_MULT_BAND_3*min(getValues(B3), na.rm = TRUE)+ RADIANCE_ADD_BAND_3
Lmin4 <- RADIANCE_MULT_BAND_4*min(getValues(B4), na.rm = TRUE)+ RADIANCE_ADD_BAND_4
Lmin5 <- RADIANCE_MULT_BAND_5*min(getValues(B5), na.rm = TRUE)+ RADIANCE_ADD_BAND_5
Lmin6 <- RADIANCE_MULT_BAND_6*min(getValues(B6), na.rm = TRUE)+ RADIANCE_ADD_BAND_6
Lmin7 <- RADIANCE_MULT_BAND_7*min(getValues(B7), na.rm = TRUE)+ RADIANCE_ADD_BAND_7

# Determinacion radiancia del objeto oscuro
LDOS2 <- 0.01*ESUN2*sin(SUN_ELEVATION_R)/pi*(d^2)
LDOS3 <- 0.01*ESUN3*sin(SUN_ELEVATION_R)/pi*(d^2)
LDOS4 <- 0.01*ESUN4*sin(SUN_ELEVATION_R)/pi*(d^2)
LDOS5 <- 0.01*ESUN5*sin(SUN_ELEVATION_R)/pi*(d^2)
LDOS6 <- 0.01*ESUN6*sin(SUN_ELEVATION_R)/pi*(d^2)
LDOS7 <- 0.01*ESUN7*sin(SUN_ELEVATION_R)/pi*(d^2)

# Determinar el efecto bruma
LP2 <- Lmin2-LDOS2
LP3 <- Lmin3-LDOS3
LP4 <- Lmin4-LDOS4
LP5 <- Lmin5-LDOS5
LP6 <- Lmin6-LDOS6
LP7 <- Lmin7-LDOS7

# Banda Reflectancia Superficie Dark Object Substrction (DOS)
RS_dos_B2 <- (pi*(L2-LP2)*(d)^2)/(ESUN2*sin(SUN_ELEVATION_R))
RS_dos_B3 <- (pi*(L3-LP3)*(d)^2)/(ESUN3*sin(SUN_ELEVATION_R))
RS_dos_B4 <- (pi*(L4-LP4)*(d)^2)/(ESUN4*sin(SUN_ELEVATION_R))
RS_dos_B5 <- (pi*(L5-LP5)*(d)^2)/(ESUN5*sin(SUN_ELEVATION_R))
RS_dos_B6 <- (pi*(L6-LP6)*(d)^2)/(ESUN6*sin(SUN_ELEVATION_R))
RS_dos_B7 <- (pi*(L7-LP7)*(d)^2)/(ESUN7*sin(SUN_ELEVATION_R))

# Combinacion de bandas
L8_B234567_RS_DOS <- stack(RS_dos_B2,RS_dos_B3,RS_dos_B4,RS_dos_B5,RS_dos_B6,RS_dos_B7)

# Exportacion imagen reflectancia Superfie DOS

writeRaster(L8_B234567_RS_DOS,"L8_B234567_RS_DOS1.tif", drivername="Gtiff", overwrite=TRUE)

## Tema: Temperatura de brillo de Landsat 8 TIRS


# Importar las imagenes
B10 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B10.TIF")
B11 = raster("LC08_L1TP_008057_20181230_20190130_01_T1_B11.TIF")

# K imagen Landsat 8 TIRS
K1_CONSTANT_BAND_10 = 774.8853
K2_CONSTANT_BAND_10 = 1321.0789
K1_CONSTANT_BAND_11 = 480.8883
K2_CONSTANT_BAND_11 = 1201.1442

# Datos de radiancia del sensor imagen Landsat
RADIANCE_MULT_BAND_10 = 3.3420E-04
RADIANCE_MULT_BAND_11 = 3.3420E-04

RADIANCE_ADD_BAND_10 = 0.10000
RADIANCE_ADD_BAND_11 = 0.10000

# Elimar valores nulos
B10[B10==0] <- NA
B11[B11==0] <- NA

# Determinancion de la Radianca
L10 <- RADIANCE_MULT_BAND_10*(na.omit(B10)) + RADIANCE_ADD_BAND_10
L11 <- RADIANCE_MULT_BAND_11*(na.omit(B11)) + RADIANCE_ADD_BAND_11


# Determinacion Temperatura de brillo Celsius
TB_B10 <- K2_CONSTANT_BAND_10/log(K1_CONSTANT_BAND_10/L10+1) - 273.15
TB_B11 <- K2_CONSTANT_BAND_11/log(K1_CONSTANT_BAND_11/L11+1) - 273.15

# Exportacion imagen Temperatura de brillo Celsius
writeRaster(TB_B10,"L8_TB_B10.tif", drivername="Gtiff",overwrite=TRUE)
writeRaster(TB_B11,"L8_TB_B11.tif", drivername="Gtiff",overwrite=TRUE)

L8_B234567_RS_DOS_py = projectRaster(L8_B234567_RS_DOS, crs = "+init=epsg:32618")

##Clasificacion supervisada con Random Forest - Landsat 8 OLI

# Agregar raster multibandas de landsat 8
L8_faca <- stack("clip_L8_B234567_RS_DOS_py.tif")

Error in .local(.Object, ...) : 

Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) : 
  Cannot create a RasterLayer object from this file. (file does not exist)

## Clasificacion supervisada con Decision_Tree - Landsat 8 OLI

# Agregar raster multibandas de landsat 8
L8_2015 <- stack("clip_L8_B234567_RS_DOS_py.tif")

Error in .local(.Object, ...) : 

Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) : 
  Cannot create a RasterLayer object from this file. (file does not exist)