Apoyo de Código - Reporte N°3

Objetivo

1. Mostrar el resultado del perfíl espectral generado por los píxeles de las diferentes categorías de clasificación y la correlación entre las bandas de las imagenes empleadas para el análisis multitemporal.

2. Generar los procesos de cálculo de los diferentes parámetros definidos por Pontius, Quantity, Exchange y Shift a partír de la matriz de confusión, para los modelos de clasificación CART y Random Forest.

library(knitr)

library(sp)
library(raster)
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/Elkin/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/Elkin/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.4-1 

library(diffeR)

Loading required package: ggplot2
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Loading required package: reshape2

Área de estudio

library(leaflet)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

tenjo <- data.frame(
    lat = c(4.871973,4.876513,4.904042,4.890152,4.812393,4.841831,4.857221,4.815281,4.871385),
    lng = c(-74.145082,-74.128094,-74.128364,-74.147039,-74.175483,-74.137124,-74.129088,-74.119443,-74.083082))
 
tenjo_popup = popup = c("Parque Principal","Casa","Peña de Juaica","Cerro Churugüaco Alto","Quebrada Socha","Quebrada Socha","Río Chicú","Serranía del Majuy","Alto de la Cruz")
 
tenjo %>% 
    leaflet() %>%
    addTiles() %>%
    addMarkers(popup = tenjo_popup, clusterOptions = markerClusterOptions())%>%
      addProviderTiles(providers$Esri.WorldStreetMap) %>%
      addMiniMap(
        tiles = providers$Esri.WorldStreetMap,
        toggleDisplay = TRUE)

NA

library(raster)
MOT <- shapefile("MOT.shp")
proj4string(MOT)

[1] "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"

MOTr <- spTransform(MOT, CRS("+init=epsg:4326"))
MOTr

class       : SpatialPolygonsDataFrame 
features    : 38 
extent      : -74.21967, -74.08091, 4.749226, 4.904683  (xmin, xmax, ymin, ymax)
crs         : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 10
names       : objectid, tipo_avalu, sector_cat,                                    reporte,    shape_Leng,    shape_Area,                                CATEGORIA,       SUELO,       AREA_ha,        area 
min values  :        0,         00,          2, Area de Actividad Agropecuaria Tradicional,             0, 28.1855385081, Area de Actividad Agropecuaria Intensiva,  PROTECCION, 0.01814860016,  0.02013345 
max values  :    70061,         00,          2,                                         ZU, 72441.4088517, 11767763.4469,                              Zona Urbana, ZONA URBANA, 7111.09671597, 53607182.05 

Polígono del MODELO DE ORDENAMIENTO TERRITORIAL, si se hace click sobre cualquier zona del polígono, se mostrará la categoría de uso de suelo a la que corresponde

par(mfrow = c(1,2))
leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%addPolygons(data = MOTr, popup = MOTr$SUELO,
    stroke = FALSE, fillOpacity = 0.3, smoothFactor = 0,
   )

Tabla de atributos del MOT

head(MOT)

#Composición en color verdadero
tenjo2 <- stack('tenjo/olitirs/TENJO_2020_ALLBANDS.tif')
tenjo1 <- stack('tenjo/olitirs/TENJO_2018_ALLBANDS.tif')
#tenjolist <- stack(bcr1,bcr2,bcr3,bcr4,bcr5,bcr6,bcr7,bcr9,bcr10,bcr11)
sen1 <- tenjo1[[c(3,2,1)]]
sen2 <- tenjo2[[c(3,2,1)]]
par(mfrow = c(1,2))
plotRGB(sen1, axes = TRUE, stretch = "lin", main = "Imagen Color Verdadero 2018")
plotRGB(sen2, axes = TRUE, stretch = "lin", main = "Imagen Color Verdadero 2020")

tenjo2

class      : RasterStack 
dimensions : 1721, 1538, 2646898, 6  (nrow, ncol, ncell, nlayers)
resolution : 10, 10  (x, y)
extent     : 586530, 601910, 524990, 542200  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      :    B2,    B3,    B4,    B8,   B11,   B12 
min values :     0,     0,     0,     0,     0,     0 
max values : 65535, 65535, 65535, 65535, 65535, 65535 

tenjo1

class      : RasterStack 
dimensions : 1721, 1538, 2646898, 6  (nrow, ncol, ncell, nlayers)
resolution : 10, 10  (x, y)
extent     : 586530, 601910, 524990, 542200  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      :    B2,    B3,    B4,    B8,   B11,   B12 
min values :     0,     0,     0,     0,     0,     0 
max values : 65535, 65535, 65535, 65535, 65535, 65535 

Correlación entre bandas

Correlación

pairs(tenjo1[[1:6]], main = "Correlación entre bandas Sentinel 2 - Reflectancia BOA 2018")

pairs(tenjo2[[1:6]], main = "Correlación entre bandas Sentinel 2 - Reflectancia BOA 2019")

Extracción de píxeles para perfil espectral - imagen de 2018

lc_samples <-shapefile('usodesuelo.shp')
class(lc_samples)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

proj4string(lc_samples)

[1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

ptsamp <- spsample(lc_samples, 300, type='regular')
ptsamp$code <- over(ptsamp, lc_samples)$code
df <- extract(tenjolist, ptsamp)
head(df)

       B2   B3   B4   B8  B11  B12
[1,] 2764 2822 2652 2236 1886 2389
[2,] 4488 4608 4508 4145 4019 4038
[3,] 3946 4144 4062 3654 3801 4121
[4,] 4090 4112 4140 3846 4735 4690
[5,] 1530 2250 2812 3208 4030 3346
[6,] 1750 2530 2964 3553 4289 3604

ms <- aggregate(df, list(ptsamp$code), mean)
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms

Perfil espectral - imagen de 2018

mycolor <- c('green', 'red', 'burlywood', 'cyan', 'black', 'pink')
ms <- as.matrix(ms)
plot(0, ylim=c(0,6000), xlim = c(1,6), type='n', xlab="Bandas", ylab = "Reflectancia")
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

title(main="Perfil Espectral de Sentinel 2_SR - Tenjo", font.main = 2)
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

Análisis matriz de confusión

Clasificaciones

lc18CART <- raster('CART_2018.tif')
lc20CART <- raster('CART_2020.tif')
lc18RF <- raster('RF_2018.tif')
lc20RF <- raster('RF_2020.tif')

plot(stretch(lc18CART, minq=0, maxq=6), main = "CART 2018", ylab="Coordenada Y", xlab="Coordenada X", col = grey(0:6 / 6))

plot(stretch(lc20CART, minq=0, maxq=6), main = "CART 2020", ylab="Coordenada Y", xlab="Coordenada X", col = grey(0:6 / 6))

plot(stretch(lc18RF, minq=0, maxq=6), main = "RF 2018", ylab="Coordenada Y", xlab="Coordenada X", col = grey(0:6 / 6))

plot(stretch(lc20RF, minq=0, maxq=6), main = "RF 2020", ylab="Coordenada Y", xlab="Coordenada X", col = grey(0:6 / 6))

Cálculos matriz de confusión

ctmatCompRefCART <- crosstabm(lc20CART, lc18CART)
ctmatCompRefRF <- crosstabm(lc20RF, lc18RF)

La matriz de confusión resultante de los datos descargados de Google Earth Engine cuenta con una columna de datos vacíos o ceros, que se evidencia en los mapas previos con el color negro, procedemos a corregír esta información

ctmatCompRef2CART <- ctmatCompRefCART[-1,-1]
ctmatCompRef2RF <- ctmatCompRefRF[-1,-1]

matriz de confusión para modelo CART

ctmatCompRef2CART

       1     2      3     4     5      6
1 165994  1310  14616  1872  1160  32549
2   6822 19766  12656  6753  1805  42153
3   6990  5334  78265  5751   884  54901
4   2863  6549   8203 33919  9327  21196
5   1788  1167   1130  7414 28577  10781
6  65409 32113 111804 19877 11146 311154

matriz de confusión para modelo RF

ctmatCompRef2RF

       1     2      3     4     5      6
1 162493  1216  11649   765  1497  28370
2   5899 18367  12851  5257   668  42125
3   5567  5796  65905  3282   399  55304
4   1505  5148   4828 40183  7275  19063
5   2617   492    706  5443 27341  10165
6  60799 34556 114344 17679 10974 353470

differenceMR(lc20CART, lc18CART, eval='original', percent=TRUE)

differenceMR(lc20RF, lc18RF, eval='original', percent=TRUE)

Métricas CART

diffTablej(ctmatCompRef2CART)

NA

Métricas RF

diffTablej(ctmatCompRef2RF)

Diagrama de componentes CART

categoryComponentsPlot(ctmatrix = ctmatCompRef2CART, units = "pixels")

Diagrama de componentes RF

categoryComponentsPlot(ctmatrix = ctmatCompRef2RF, units = "pixels")