Practica 3 - Exploración
Ejercicio desarrollado en Marzo 11 de 2020, en el curso de Sensores remotos, y orientado por la practica de la página https://rspatial.org/raster/rs/2-exploration.html
Llamamos los datos del ejercicio
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')
}Se crean los elementos tipo Raster para las bandas b2, b3, b4 y b5
library(raster)## Loading required package: sp# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')“————————————” PROPIEDADES DE LAS IMAGENES “————————————”
Se puede comprobar los atributos de elemento Raster
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/usuario/R/data/rs/LC08_044034_20170614_B2.tif
## names : LC08_044034_20170614_B2
## values : 0.0748399, 0.7177562 (min, max)## 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:/github/rspatial/rspatial-web/source/rs/_R/data/rs/LC08_044034_20170614_B2.tif
## names : LC08_044034_20170614_B2
## values : 0.0748399, 0.7177562 (min, max)“————————————” INFORMACIÓN DE IMAGÉN Y ESTADISTICAS “————————————”
Se recupera la información del sistema coordenado para la banda B2, coordinate reference system CRS (elemento raster b2)
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, columnsinformación del número de celdas y dimensiones, resolución espacial, Número de layers
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] 1se puede comparar 2 bandas, evaluando los valores similares de resolucion, proyeccion, número de columnas, entre otras
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)## [1] TRUE## [1] TRUESe puede crear un elemento tipo RasterStack, un formato que contiene varias bandas de información al simultaneo
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, 0.6924697Tambien se puede crear el formato tipo Stack, partiendo de buscar su localización
# 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"
## [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"## [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“————————————” COMPOSICIÓN DE MAPAS CON UNA SOLA BANDA “————————————”
Ploteamos las bandas, la función “par” hace referencia a los parametros del ploteo que serán modificados. “mflow = c(2,2)” se dividira el espacio de ploteo en pantalla en 2 filas y 2 columas.
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))
Se va a plotear un elemento tipo Stack, que se creara de la uno con la banda b4, b3 y b2 en RGB (Color verdadero)
landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
probamos con una imagen en Falso color (b5, b4, b3)
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
“————————————” SUBCONJUNTOS Y RENOMBRADO DE BANDAS “————————————”
# 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] 3landsat <- subset(landsat, 1:7)names(landsat)## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"“————————————” SPATIAL SUBSET OR CROP “————————————”
# 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)“————————————” SAVING RESULTS “————————————”
x <- writeRaster(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE)writeRaster(landsatcrop, filename="cropped-landsat.grd", overwrite=TRUE)“————————————” RELATION BETWEEN BANDS “————————————”
pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")
pairs(landsatcrop[[4:5]], main = "Red versus NIR")
“————————————” EXTRACT PIXEL VALUES “————————————”
# 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 red NIR SWIR1 SWIR2
## [1,] 0.1418727 0.1362776 0.14601481 0.18632990 0.3224774 0.3689297 0.2452520
## [2,] 0.1355836 0.1175188 0.10162266 0.10088532 0.1641013 0.2078428 0.1718000
## [3,] 0.1328078 0.1149815 0.09737211 0.09559382 0.1675928 0.2021393 0.1657061
## [4,] 0.1341740 0.1295331 0.13629928 0.18550581 0.3426891 0.3296556 0.1909491
## [5,] 0.1330029 0.1293379 0.13820769 0.18489859 0.3458120 0.3126317 0.1776336
## [6,] 0.1484003 0.1487039 0.16887231 0.22395587 0.3541830 0.4215410 0.2786925## 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“————————————” SPECTRAL PROFILES “————————————”
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# 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")