imagen california

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

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

debo llamar los datos

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

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.

G:

getwd() 

definio la ruta

library (raster)

# Blue
b2 = raster('./rsdata/data/rs/LC08_044034_20170614_B2.tif')
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     : G:/Cursos/PR/rsdata/data/rs/LC08_044034_20170614_B2.tif 
names      : LC08_044034_20170614_B2 
values     : 0.0748399, 0.7177562  (min, max)

# Green
b3 = raster('./rsdata/data/rs/LC08_044034_20170614_B3.tif')
b3

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     : G:/Cursos/PR/rsdata/data/rs/LC08_044034_20170614_B3.tif 
names      : LC08_044034_20170614_B3 
values     : 0.04259216, 0.6924697  (min, max)

# Red
b4 = raster('./rsdata/data/rs/LC08_044034_20170614_B4.tif')
b4

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     : G:/Cursos/PR/rsdata/data/rs/LC08_044034_20170614_B4.tif 
names      : LC08_044034_20170614_B4 
values     : 0.02084067, 0.7861769  (min, max)

Near Infrared (NIR)

b5 = raster(‘./rsdata/data/rs/LC08_044034_20170614_B5.tif’) b2

# Near Infrared (NIR)
b5 = raster('./rsdata/data/rs/LC08_044034_20170614_B5.tif')
b5

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     : G:/Cursos/PR/rsdata/data/rs/LC08_044034_20170614_B5.tif 
names      : LC08_044034_20170614_B5 
values     : 0.0008457669, 1.012432  (min, max)

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     : G:/Cursos/PR/rsdata/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)

# 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

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

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

 [1] "./rsdata/data/rs/LC08_044034_20170614_B1.tif"  "./rsdata/data/rs/LC08_044034_20170614_B2.tif"  "./rsdata/data/rs/LC08_044034_20170614_B3.tif" 
 [4] "./rsdata/data/rs/LC08_044034_20170614_B4.tif"  "./rsdata/data/rs/LC08_044034_20170614_B5.tif"  "./rsdata/data/rs/LC08_044034_20170614_B6.tif" 
 [7] "./rsdata/data/rs/LC08_044034_20170614_B7.tif"  "./rsdata/data/rs/LC08_044034_20170614_B8.tif"  "./rsdata/data/rs/LC08_044034_20170614_B9.tif" 
[10] "./rsdata/data/rs/LC08_044034_20170614_B10.tif" "./rsdata/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,      

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

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

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

# 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

landsat <- subset(landsat, 1:7)

names(landsat)

[1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2" "LC08_044034_20170614_B3" "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"

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

x <- writeRaster(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE)

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

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

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

# load the polygons with land use land cover information
samp <- readRDS('./rsdata/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.1394872 0.1241549 0.1140490 0.1182128 0.1937900 0.2363822 0.1976068
[2,]  0.1363426 0.1197091 0.1095816 0.1129864 0.1844866 0.2381171 0.2014020
[3,]  0.1354318 0.1336969 0.1478582 0.2043296 0.3610792 0.3619901 0.2104669
[4,]  0.1269958 0.1180827 0.1185164 0.1564026 0.2825093 0.3155811 0.1909925
[5,]  0.1429570 0.1458413 0.1698265 0.2372280 0.3815729 0.4079219 0.2431484
[6,]  0.1380125 0.1214441 0.1093647 0.1094515 0.1923370 0.2201391 0.1791083

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

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