library(raster)
## Loading required package: sp
raslist <- paste0('C:/percepci/DATOS_TALLER/Lansat_8_F/Banda', 1:7, ".tif")
landsat <- stack(raslist)
landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]
par(mfrow = c(1,2))
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Verdadero color")
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Falso Color RGB 543")
par(mfrow = c(1,2))
B2 <- subset(landsat, 2)
B4 <- subset(landsat, 4)
B5 <- subset(landsat, 5)
B6 <- subset(landsat, 6)
ndvi <- (B5-B4)/(B5+B4)
BSI <- ((B6 + B4) - (B5 + B2)) / ((B6 + B4) + (B5 + B2))
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")
plot(BSI, col = rev(heat.colors(10)), main = "Índice de suelo desnudo (BSI)")
hist(ndvi,
main = "Distribution of NDVI values",
xlab = "NDVI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
library(sp)
library(rgdal)
## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/jcami/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/jcami/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1
samp <- readOGR(dsn = "C:/percepci/Shape", layer = "Uso")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\percepci\Shape", layer: "Uso"
## with 38 features
## It has 2 fields
## Integer64 fields read as strings: Id
# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
# add the land cover class to the points
ptsamp$Uso <- over(ptsamp, samp)$Uso
# extract values with points
df <- extract(landsat, ptsamp)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster
# To see some of the reflectance values
head(df)
## Banda1 Banda2 Banda3 Banda4 Banda5 Banda6
## [1,] 0.12461501 0.11716430 0.11497430 0.12973291 0.1886483 0.21738000
## [2,] 0.10200102 0.08355275 0.08636165 0.07695898 0.2059302 0.18131661
## [3,] 0.08224353 0.06389048 0.05027448 0.03237371 0.2187368 0.10595252
## [4,] 0.08148179 0.06305734 0.04703709 0.02944580 0.2333526 0.09828757
## [5,] 0.08155319 0.06315254 0.04970317 0.03080264 0.2486825 0.10464329
## [6,] 0.08169603 0.06281929 0.04975079 0.03070742 0.2256639 0.10269133
## Banda7
## [1,] 0.19919358
## [2,] 0.11825929
## [3,] 0.04518038
## [4,] 0.04201442
## [5,] 0.04263333
## [6,] 0.04501374
ms <- aggregate(df, list(ptsamp$Uso), mean)
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
## Banda1 Banda2 Banda3 Banda4 Banda5
## Arbustal 0.08236849 0.06701478 0.05377965 0.04484117 0.1265670
## Bosque denso 0.08529870 0.06645355 0.05178696 0.03392465 0.2234194
## Cuerpos de agua 0.07836344 0.06001040 0.03463511 0.01849586 0.0107833
## Ecosistema de paramo 0.08640014 0.07341368 0.06397692 0.06265435 0.1461983
## Nubes 0.46551502 0.46387255 0.45066121 0.46591970 0.5615650
## Pastos y cultivos 0.09633561 0.08500657 0.08701935 0.08356951 0.2692088
## Plantaciones forestales 0.08778990 0.07087464 0.05669685 0.04837966 0.1626684
## Suelo Desnudo 0.16578437 0.16492744 0.18163799 0.21463061 0.2547407
## Tejido urbano 0.11927543 0.10896374 0.10757268 0.11670903 0.2120419
## Tierras degradadas 0.12822432 0.11889307 0.13176222 0.14651191 0.2282823
## Vegetacion herbacea 0.10088398 0.08588997 0.07992879 0.06998348 0.2069441
## Banda6 Banda7
## Arbustal 0.11873537 0.070466390
## Bosque denso 0.10320725 0.044468082
## Cuerpos de agua 0.01194971 0.008498097
## Ecosistema de paramo 0.20071351 0.129760432
## Nubes 0.51966959 0.431832105
## Pastos y cultivos 0.23986527 0.142998906
## Plantaciones forestales 0.11719286 0.066356595
## Suelo Desnudo 0.32391575 0.222938292
## Tejido urbano 0.23312053 0.190871044
## Tierras degradadas 0.26899049 0.187226022
## Vegetacion herbacea 0.19888545 0.122077215
mycolor <- c('darkred', 'green', 'blue', 'cyan', 'burlywood', "red", "Orange", "yellow","purple", "gray", "black")
#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="Bandas", ylab = "Reflectancia")
# 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="Perfiles espectrales", font.main = 2)
# Legend
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
library(raster)
nlcd <- brick('C:/percepci/insumos/Usos_ras.tif')
names(nlcd) <- c("nlcd2011")
# The class names and colors for plotting
nlcdclass <- c("Cuerpos de agua", "Tejido urbano", "Ecosistema de paramo", "Pastos y cultivos", "Bosque denso", "Arbustal", "Tierras degradadas", "Vegetacion herbacea","Suelo Desnudo","Plantaciones forestales", "Nubes")
classdf <- data.frame(classvalue1 = c(1,2,3,4,5,6,7,8,9,10,11), classnames1 = nlcdclass)
# Hex codes of colors
classcolor <- c("#5475A8", "#B50000", "green", "#38814E", "#AF963C", "#D1D182", "#FBF65D", "#C8E6F8","black","yellow","gray")
# Now we ratify (RAT = "Raster Attribute Table") the ncld2011 (define RasterLayer as a categorical variable). This is helpful for plotting.
nlcd2011 <- nlcd[[1]]
nlcd2011 <- ratify(nlcd2011)
rat <- levels(nlcd2011)[[1]]
#
rat$landcover <- nlcdclass
levels(nlcd2011) <- rat
set.seed(99)
# Sampling
samp2011 <- sampleStratified(nlcd2011, size = 200, na.rm = TRUE, sp = TRUE)
## Warning in .local(x, size, ...): only 193 cells found for stratum 1
## Warning in .local(x, size, ...): only 145 cells found for stratum 9
## Warning in .local(x, size, ...): only 42 cells found for stratum 11
## Warning in .local(x, size, ...): fewer samples than requested found for strata:
## 1, 9, 11
samp2011
## class : SpatialPointsDataFrame
## features : 1980
## extent : 738591.6, 798561.6, 636142.9, 684802.9 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 2
## names : cell, nlcd2011
## min values : 9407, 1
## max values : 3265182, 11
table(samp2011$nlcd2011)
##
## 1 2 3 4 5 6 7 8 9 10 11
## 193 200 200 200 200 200 200 200 145 200 42
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
plt <- levelplot(nlcd2011, col.regions = classcolor, main = 'Distribution of Training Sites')
print(plt + layer(sp.points(samp2011, pch = 3, cex = 0.5, col = 1)))
sampvals <- extract(landsat, samp2011, df = TRUE)
sampvals <- sampvals[, -1]
sampdata <- data.frame(classvalue = samp2011@data$nlcd2011, sampvals)
library(rpart)
# Train the model
cart <- rpart(as.factor(classvalue)~., data=sampdata, method = 'class', minsplit = 5)
# print(model.class)
# Plot the trained classification tree
plot(cart, uniform=TRUE, main="Classification Tree")
text(cart, cex = 0.8)
pr2011 <- predict(landsat, cart, type='class')
pr2011
## class : RasterLayer
## dimensions : 2001, 2102, 4206102 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 737191.6, 800251.6, 634019.1, 694049.1 (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, 11 (min, max)
## attributes :
## ID value
## from: 1 1
## to : 11 11
pr2011 <- ratify(pr2011)
rat <- levels(pr2011)[[1]]
rat$legend <- classdf$classnames
levels(pr2011) <- rat
levelplot(pr2011, maxpixels = 1e6,
col.regions = classcolor,
scales=list(draw=FALSE),
main = "Clasificación supervisada")
# Extract the layer values for the locations
sampvals1 <- extract(landsat, ptsamp, df = TRUE)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster
# sampvals no longer has the spatial information. To keep the spatial information you use `sp=TRUE` argument in the `extract` function.
# drop the ID column
sampvals1 <- sampvals1[, -1]
# combine the class information with extracted values
sampdata1 <- data.frame(classvalue = ptsamp$Uso, sampvals1)
sam1<- readOGR("C:/percepci/Shape", layer = "Uso")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\percepci\Shape", layer: "Uso"
## with 38 features
## It has 2 fields
## Integer64 fields read as strings: Id
r <- raster(landsat)
rp <- rasterize(sam1, r, field = "Uso")