Universidad Cooperativa de Colombia
# install.packages("raster")
# install.packages("sp")
# install.packages("rgdal")
# install.packages("readxl")
# install.packages("sf")
# install.packages("writexl")
library(raster)
## Loading required package: sp
library(sp)
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/User/OneDrive - Universidad Cooperativa de Colombia/Documents/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/User/OneDrive - Universidad Cooperativa de Colombia/Documents/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/User/OneDrive - Universidad Cooperativa de Colombia/Documents/R/win-library/4.1/rgdal/proj
library(readxl)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(writexl)
#se le da la lectura a las bandas
b2 <- raster ("C:/Users/User/OneDrive/PRACTICA_03/T18NVH_20170921T152631_B02.JP2")
b2
## class : RasterLayer
## dimensions : 10980, 10980, 120560400 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : T18NVH_20170921T152631_B02.jp2
## names : T18NVH_20170921T152631_B02
## values : 0, 65535 (min, max)
b3 <- raster ("C:/Users/User/OneDrive/PRACTICA_03/T18NVH_20170921T152631_B03.JP2")
b3
## class : RasterLayer
## dimensions : 10980, 10980, 120560400 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : T18NVH_20170921T152631_B03.jp2
## names : T18NVH_20170921T152631_B03
## values : 0, 65535 (min, max)
b4 <- raster ("C:/Users/User/OneDrive/PRACTICA_03/T18NVH_20170921T152631_B04.JP2")
b4
## class : RasterLayer
## dimensions : 10980, 10980, 120560400 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : T18NVH_20170921T152631_B04.jp2
## names : T18NVH_20170921T152631_B04
## values : 0, 65535 (min, max)
b8 <- raster ("C:/Users/User/OneDrive/PRACTICA_03/T18NVH_20170921T152631_B08.JP2")
b8
## class : RasterLayer
## dimensions : 10980, 10980, 120560400 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : T18NVH_20170921T152631_B08.jp2
## names : T18NVH_20170921T152631_B08
## values : 0, 65535 (min, max)
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 ="reg", col = gray(0:100 / 100))
plot (b8, main = "NIR", col = gray(0:100 / 100))
#se realiza la conbinacion de todas las bandas
sentinel <- stack(b2,b3,b4,b8)
sentinelRGB <- stack(b4,b3,b2)
# se plotea el sentinelRGB
plotRGB(sentinelRGB,r=3, g=2, b=1, axes= TRUE, stretch = "lin", main = "sentinel true color ")
extent(sentinel)
## class : Extent
## xmin : 399960
## xmax : 509760
## ymin : 190200
## ymax : 3e+05
area_recortada <- extent(399960, 469000, 238000, 3e+05)
sentinelcortada <- crop(sentinel,area_recortada)
plot(sentinelcortada)
plotRGB(sentinelcortada, r=3, g=2, b=1, axes = TRUE, stretch = "lin", main = "sentinel True Color")
#se le asignan nombres a las bandas
names(sentinel) <- c("blue","green","red","NIR")
sentinel
## class : RasterStack
## dimensions : 10980, 10980, 120560400, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : blue, green, red, NIR
## min values : 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535
preparacion <- read_xlsx("C:/Users/User/OneDrive/PRACTICA_03/Coordenadas.xlsx")
preparacion <- preparacion %>%
st_as_sf(coords = c("x","y"))
plot(preparacion$geometry)
df <- extract(sentinel, preparacion); df
## blue green red NIR
## [1,] 1314 1086 892 858
## [2,] 1192 993 823 564
## [3,] 1384 1199 1155 1001
## [4,] 1119 877 630 483
## [5,] 1185 968 775 681
## [6,] 1188 966 774 704
## [7,] 1171 945 733 617
## [8,] 1144 934 740 639
## [9,] 1204 1048 933 755
## [10,] 1209 1016 896 744
## [11,] 1565 1625 2015 3032
## [12,] 1381 1243 1067 2764
## [13,] 1383 1303 1415 2327
## [14,] 1212 1104 1039 1643
## [15,] 1544 1672 2113 2624
## [16,] 1369 1324 1386 2371
## [17,] 1185 1017 966 1523
## [18,] 1312 1220 1294 2539
## [19,] 1636 1683 2089 2830
## [20,] 1449 1481 1814 2552
## [21,] 1189 1101 802 2189
## [22,] 1340 1294 1311 2515
## [23,] 1123 966 717 3437
## [24,] 1136 900 762 852
## [25,] 797 602 399 904
## [26,] 1027 929 689 2560
## [27,] 1181 1047 1052 1956
## [28,] 1143 962 919 1468
## [29,] 1082 961 739 2269
## [30,] 1100 922 795 1691
## [31,] 1415 1401 1357 2514
## [32,] 1676 1581 1743 2082
## [33,] 1686 1595 1584 2574
## [34,] 1815 1700 1674 2775
## [35,] 2048 2111 2325 2839
## [36,] 1575 1506 1459 2709
## [37,] 1948 2032 2262 2934
## [38,] 1547 1489 1391 1800
## [39,] 1382 1286 1307 1898
## [40,] 1785 1779 1906 2276
ms <- aggregate(df, list(preparacion$clase), mean); ms
## Group.1 blue green red NIR
## 1 agua 1211.0 1003.2 835.1 704.6
## 2 construccion 1687.7 1648.0 1700.8 2440.1
## 3 suelo 1403.6 1367.2 1519.8 2420.5
## 4 vegetacion 1111.8 968.4 818.5 1984.1
write_xlsx(ms, "C:/Users/User/OneDrive/PRACTICA_03/datos_para_la_grafica.xlsx")
#fin de la practica 3