Imagen LandSat 8, California EE.UU _ DAGZ

Landsat 8 recopilada el 14 de junio de 2017. El subconjunto cubre el area entre Concord y Stockton, en California, EE. UU.

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/PERCEPCION REMOTA
getwd()

[1] "C:/Users/David/Desktop/PERCEPCION REMOTA"

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

Se inserta un nuevo bloque de codigo ctrl+alt+i para descargar los datos con los que voy a trabajar

dir.create('data', showWarnings = FALSE)
if (!file.exists('data/rs/samples.rds')) {
    download.file('https://biogeo.ucdavis.edu/data/rspatial/rsdata.zip', dest = './data/rsdata.zip')
    unzip('./data/rsdata.zip', exdir='data')}

Crear objetos RasterLayer para capas individuales Landsat (bandas) y las llamamos enrutandolas como por ejemplo (‘./data/rs/LC08_044034_20170614_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)
# azul
b2 = raster('./data/rs/LC08_044034_20170614_B2.tif')
# verde
b3 = raster('./data/rs/LC08_044034_20170614_B3.tif')
# rojo
b4 = raster('./data/rs/LC08_044034_20170614_B4.tif')
# Infrarojo cercano (NIR)
b5 = raster('./data/rs/LC08_044034_20170614_B5.tif')

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

b2

class      : RasterLayer 
dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/David/Desktop/PERCEPCION REMOTA/data/rs/LC08_044034_20170614_B2.tif 
names      : LC08_044034_20170614_B2 
values     : 0.0748399, 0.7177562  (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=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

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

[1] 1863765

## [1] 1863765
dim(b2)

[1] 1245 1497    1

## [1] 1245 1497    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

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

class      : RasterStack 
dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3 
min values :            0.0008457669,            0.0208406653,            0.0425921641 
max values :               1.0124315,               0.7861769,               0.6924697 

## class      : RasterStack
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3
## min values :            0.0008457669,            0.0208406653,            0.0425921641
## max values :               1.0124315,               0.7861769,    

Tambien puede crear el RasterStack usando los nombres de archivo

# first create a list of raster layers to use
filenames = paste0(' data/rs/LC08_044034_20170614_B', 1:11, ".tif")
filenames

 [1] " data/rs/LC08_044034_20170614_B1.tif"  " data/rs/LC08_044034_20170614_B2.tif" 
 [3] " data/rs/LC08_044034_20170614_B3.tif"  " data/rs/LC08_044034_20170614_B4.tif" 
 [5] " data/rs/LC08_044034_20170614_B5.tif"  " data/rs/LC08_044034_20170614_B6.tif" 
 [7] " data/rs/LC08_044034_20170614_B7.tif"  " data/rs/LC08_044034_20170614_B8.tif" 
 [9] " data/rs/LC08_044034_20170614_B9.tif"  " data/rs/LC08_044034_20170614_B10.tif"
[11] " data/rs/LC08_044034_20170614_B11.tif"

##  [1] "./data/rs/LC08_044034_20170614_B1.tif"
##  [2] "./data/rs/LC08_044034_20170614_B2.tif"
##  [3] "./data/rs/LC08_044034_20170614_B3.tif"
##  [4] "./data/rs/LC08_044034_20170614_B4.tif"
##  [5] "./data/rs/LC08_044034_20170614_B5.tif"
##  [6] "./data/rs/LC08_044034_20170614_B6.tif"
##  [7] "./data/rs/LC08_044034_20170614_B7.tif"
##  [8] "./data/rs/LC08_044034_20170614_B8.tif"
##  [9] "./data/rs/LC08_044034_20170614_B9.tif"
## [10] "./data/rs/LC08_044034_20170614_B10.tif"
## [11] "./data/rs/LC08_044034_20170614_B11.tif"
landsat = stack(filenames)
landsat

class      : RasterStack 
dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029 

## class      : RasterStack
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029

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

There were 30 warnings (use warnings() to see them)

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

landsatRGB = stack(b4, b3, b2)
plotRGB(landsatRGB, 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(landsatRGB, axes=TRUE, stretch="lin", main="Landsat True Color Composite")
landsatFCC = stack(b5, b4, b3)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

Nota : 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
landsatsub1 = subset(landsat, 1:3)
# same
landsatsub2 = landsat[[1:3]]
# Number of bands in the original and new data
nlayers(landsat)

[1] 11

## [1] 11
nlayers(landsatsub1)

[1] 3

## [1] 3
nlayers(landsatsub2)

[1] 3

## [1] 3

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

landsat = subset(landsat, 1:7)

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

names(landsat)

[1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2" "LC08_044034_20170614_B3"
[4] "LC08_044034_20170614_B4" "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
[7] "LC08_044034_20170614_B7"

## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) = c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)

[1] "ultra.blue" "blue"       "green"      "red"        "NIR"        "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(landsat)

class      : Extent 
xmin       : 594090 
xmax       : 639000 
ymin       : 4190190 
ymax       : 4227540 

## class      : Extent
## xmin       : 594090
## xmax       : 639000
## ymin       : 4190190
## ymax       : 4227540
e = extent(624387, 635752, 4200047, 4210939)
# crop landsat by the extent
landsatcrop = crop(landsat, e)

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(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE)

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

writeRaster(landsatcrop, 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(landsatcrop[[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(landsatcrop[[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
samp = readRDS('data/rs/samples.rds')
# generate 300 point samples from the polygons
ptsamp = spsample(samp, 300, type='regular')
# add the land cover class to the points
ptsamp$class = over(ptsamp, samp)$class
# extract values with points
df = extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)

     ultra.blue      blue      green
[1,]  0.1351065 0.1173020 0.10157929
[2,]  0.1326343 0.1151767 0.09889016
[3,]  0.1356270 0.1332415 0.14410640
[4,]  0.1454509 0.1424365 0.16008930
[5,]  0.1598074 0.1680049 0.19600205
[6,]  0.1400944 0.1241549 0.11179360
            red       NIR     SWIR1
[1,] 0.10047328 0.1631471 0.2042646
[2,] 0.09769741 0.1678530 0.2016188
[3,] 0.19365992 0.3442072 0.3302411
[4,] 0.20805971 0.3220870 0.3971004
[5,] 0.26268786 0.4421865 0.4378276
[6,] 0.11508994 0.1862865 0.2442327
         SWIR2
[1,] 0.1674193
[2,] 0.1650121
[3,] 0.1929660
[4,] 0.2689119
[5,] 0.2693673
[6,] 0.2109006

##      ultra.blue      blue      green        red       NIR     SWIR1
## [1,]  0.1367547 0.1197091 0.10429009 0.10507080 0.1670290 0.2161921
## [2,]  0.1343041 0.1163694 0.09889016 0.09752392 0.1686988 0.2066501
## [3,]  0.1383812 0.1375354 0.15377855 0.20988137 0.3602552 0.3594528
## [4,]  0.1293813 0.1254127 0.13582218 0.18546245 0.3094872 0.2950440
## [5,]  0.1481184 0.1531496 0.17986734 0.24896033 0.3882957 0.4010257
## [6,]  0.1342608 0.1158490 0.10029978 0.09932390 0.1649471 0.2108356
##          SWIR2
## [1,] 0.1817324
## [2,] 0.1710843
## [3,] 0.2157801
## [4,] 0.1653591
## [5,] 0.2454254
## [6,] 0.1800408

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.

ms = aggregate(df, list(ptsamp$class), mean)
# instead of the first column, we use row names
rownames(ms) = ms[,1]
ms = ms[,-1]
ms

##          ultra.blue      blue      green       red        NIR      SWIR1
## built     0.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')
#transform ms from a data.frame to a matrix
ms = as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

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