Practica 3

UNIVERSIDAD COOPERATIVA DE COLOMBIA

CAMPUS NEIVA

#install.packages("raster")
library(raster)

## Loading required package: sp

#install.packages("sp")
library(sp)
#install.packages("rgdal")
library(rgdal)

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/mcra9/OneDrive/Documentos/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/mcra9/OneDrive/Documentos/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/mcra9/OneDrive/Documentos/R/win-library/4.1/rgdal/proj

Se llaman las imagenes

b2 <- raster("C:/Users/Public/Andres/electiva_2/T18NVH_20170921T152631_B02.jp2")
b3 <- raster("C:/Users/Public/Andres/electiva_2/T18NVH_20170921T152631_B03.jp2")
b4 <- raster("C:/Users/Public/Andres/electiva_2/T18NVH_20170921T152631_B04.jp2")
b8 <- raster("C:/Users/Public/Andres/electiva_2/T18NVH_20170921T152631_B08.jp2")

Se plotean las bandas

par(mfrow = c(1,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(b8, main = "NIR", col = gray(0:100 / 100))

Combinación de bandas - color verdadero

sentinel <- stack(b2,b3,b4,b8)

plotRGB(sentinel, axes = TRUE, stretch = "lin", main = "Sentinel true color composite")

combinación de bandas - color falso

sentinelFCC <- stack(b8, b4, b3, b2)

plotRGB(sentinelFCC, axes = TRUE, stretch = "lin", main = "Sentinel false color composite") 

Recorte del area de estudio

extent(sentinel)

## class      : Extent 
## xmin       : 399960 
## xmax       : 509760 
## ymin       : 190200 
## ymax       : 3e+05

recorte <- extent(399960, 469000, 238000, 3e+05)

recortesentinel <- crop(sentinel,recorte)

recortesentinel

## class      : RasterBrick 
## dimensions : 6200, 6904, 42804800, 4  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 399960, 469000, 238000, 3e+05  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2021-09-28_205529_14564_91555.grd 
## names      : T18NVH_20170921T152631_B02, T18NVH_20170921T152631_B03, T18NVH_20170921T152631_B04, T18NVH_20170921T152631_B08 
## min values :                        154,                        127,                          1,                          1 
## max values :                      13166,                      13880,                      18724,                      18573

plotRGB(recortesentinel, r=3, g=2, b=1,  axes = TRUE, stretch = "lin", main = "Recorte Sentinel true color composite") 

Renombrar las bandas

names(sentinel) <- c("BLUE", "GREEN", "RED", "NIR" )

agregar el excel

#install.packages("readxl")
library(readxl)
#install.packages("sf")
library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

entrenamiento <- read_xlsx("C:/Users/Public/Andres/electiva_2/puntos.xlsx")

entrenamiento <- entrenamiento %>%
  st_as_sf(coords = c("X", "Y")) 

plot(entrenamiento$geometry)

df <- extract(sentinel, entrenamiento); df

##       BLUE GREEN  RED  NIR
##  [1,] 1122   878  631  478
##  [2,] 1376  1275 1275 1464
##  [3,] 1584  1546 1683 2110
##  [4,] 1121   934  702  813
##  [5,] 1207  1145  960 2835
##  [6,] 1307  1244 1147 1455
##  [7,] 1133   951  694  966
##  [8,] 1195  1040  795  860
##  [9,] 1342  1144 1004  824
## [10,] 1150   923  738  698
## [11,] 1575  1547 1894 2615
## [12,] 1186  1053 1091 2124
## [13,] 1295  1190 1173 2261
## [14,] 1244  1135 1290 2314
## [15,] 1014   975  710 2904
## [16,]  972   936  713 2553
## [17,]  973   915  783 3757
## [18,] 1056  1038  816 3099
## [19,] 1289  1266 1354 2157
## [20,]  932   946  678 2934
## [21,]  994   890  592 2664
## [22,]  994   851  587 2704
## [23,] 1004   861  640 2449
## [24,] 1057   942  653 3142
## [25,]  755   668  448 2553
## [26,] 2025  1830 1609 3221
## [27,] 1012   922  654 2833
## [28,] 1053   934  739 2470
## [29,]  952   858  583 3061
## [30,] 1104   935  755 1713
## [31,] 1688  1689 1748 2388
## [32,] 1486  1436 1279 2949
## [33,] 1943  1938 1893 2333
## [34,] 3667  4026 4706 5416
## [35,] 1336  1130 1105 1862
## [36,] 1652  1185 1107 2632
## [37,] 1541  1193 1140 2152
## [38,] 1676  1765 2128 2759
## [39,] 1215  1121 1577 2487
## [40,] 1236  1115 1577 2243

extraer valor promedio de bandas

ms <- aggregate(df, list(entrenamiento$clase), mean); ms

##        Group.1   BLUE  GREEN    RED    NIR
## 1         agua 1253.7 1108.0  962.9 1250.3
## 2 construccion 1744.0 1659.8 1826.0 2722.1
## 3        suelo 1153.6 1100.1 1050.2 2671.8
## 4   vegetacion 1095.0  969.1  726.0 2681.0

Exportar valores a excel

#install.packages("writexl")
library(writexl)
write_xlsx (ms, "C:/Users/Public/Andres/electiva_2/promedio.xlsx")