PRACTICA 3
Laura Dussan y Cristian Guarnizo
23/9/2021
#PRACTICA 3
###Instalar Paquetes
#install.packages("raster")#install.packages("sp")#install.packages("rgdal")#Paquetes instalados
library(sp)
library(raster)
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/ASUS/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/ASUS/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/ASUS/Documents/R/win-library/4.1/rgdal/proj#Instalar los objetos a mis variables
b2 <-raster("C:/Users/ASUS/Downloads/PRACTICA 3/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/ASUS/Downloads/PRACTICA 3/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/ASUS/Downloads/PRACTICA 3/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/ASUS/Downloads/PRACTICA 3/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)#Ploteo composicion a color verdadero
sentinel <-stack(b2, b3, b4, b8)
sentinelRGB <-stack(b4, b3, b2)
plotRGB(sentinelRGB, r = 3, g = 2, b = 1, axes = TRUE, stretch = "hist", main = "Sentinel True Color Composite")
###Area de recorte y ploteo a color verdadero
area_recorte <- extent(400000, 5000000, 230000, 300000)
#recorte de mi imagen o stack sentinel
recorte <-crop(sentinel, area_recorte)
recorte## class : RasterBrick
## dimensions : 7000, 10976, 76832000, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 4e+05, 509760, 230000, 3e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2021-09-27_193018_12208_85821.grd
## names : T18NVH_20170921T152631_B02, T18NVH_20170921T152631_B03, T18NVH_20170921T152631_B04, T18NVH_20170921T152631_B08
## min values : 154, 127, 1, 1
## max values : 15367, 13880, 18724, 18743plotRGB(recorte, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "Sentinel Subset True Color Composite")
#Renombrar las bandas
names(recorte) <- c("blue", "green", "red", "NIR")#Instalar paquetes faltantes
#install.packages("readxl")
#install.packages("sf")#Integrar el Excel
library(readxl) # leer Excel
library(sf) #Manejar datos geograficos## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1datos <- read_xlsx("C:/Users/ASUS/Documents/TATI 2021-1/ELECTIVA 2/PUNTOS ELECTIVA 2.xlsx")
datos <- datos %>%
st_as_sf(coords = c("X", "Y"))
datos## Simple feature collection with 60 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 407473.1 ymin: 238716.4 xmax: 467379.8 ymax: 297330.3
## CRS: NA
## # A tibble: 60 x 4
## NUMERO `ID. CLASE` CLASE geometry
## <dbl> <dbl> <chr> <POINT>
## 1 1 1 AGUA (423794.6 239479.1)
## 2 2 1 AGUA (418243.3 238716.4)
## 3 3 1 AGUA (415632.4 240996.6)
## 4 4 1 AGUA (439158.5 264713.7)
## 5 5 1 AGUA (442228.7 292964.7)
## 6 6 1 AGUA (451472.6 294446.4)
## 7 7 1 AGUA (463268.9 295218.5)
## 8 8 1 AGUA (443453.3 242467.5)
## 9 9 1 AGUA (432302.4 255623.3)
## 10 10 1 AGUA (434004.4 260376.1)
## # ... with 50 more rowsplot(datos$geometry)
df <- extract(recorte, datos); df## blue green red NIR
## [1,] 1257 1129 960 1307
## [2,] 1130 943 759 909
## [3,] 1142 1058 835 2976
## [4,] 1234 1163 1130 2739
## [5,] 1599 1677 1945 2562
## [6,] 1137 882 662 467
## [7,] 1197 1019 822 1328
## [8,] 1029 970 742 2427
## [9,] 1160 961 766 882
## [10,] 1194 970 784 715
## [11,] 1179 958 754 636
## [12,] 1196 984 839 796
## [13,] 1225 1170 1122 2819
## [14,] 1604 1601 1540 2955
## [15,] 1169 1018 897 923
## [16,] 1086 980 786 2358
## [17,] 1063 964 836 2192
## [18,] 1324 1227 1113 2227
## [19,] 1148 1009 855 1990
## [20,] 1025 945 624 3734
## [21,] 1322 1215 1364 2415
## [22,] 1366 1271 1412 2187
## [23,] 1529 1512 1579 2451
## [24,] 1362 1379 1650 2312
## [25,] 1429 1367 1392 2825
## [26,] 1135 1031 940 2139
## [27,] 952 879 746 3130
## [28,] 1473 1454 1635 2688
## [29,] 1444 1458 1604 2590
## [30,] 1370 1401 1487 2879
## [31,] 1043 976 674 3137
## [32,] 1143 1040 854 2770
## [33,] 1134 1001 758 2604
## [34,] 1008 868 608 3157
## [35,] 901 789 533 2763
## [36,] 766 691 430 3004
## [37,] 809 693 473 2201
## [38,] 996 884 641 2265
## [39,] 1056 896 650 2437
## [40,] 1109 986 736 2754
## [41,] 1134 1040 704 3128
## [42,] 1175 1015 936 1819
## [43,] 1126 975 777 2073
## [44,] 903 786 537 3048
## [45,] 1077 909 792 1550
## [46,] 1266 1159 977 2436
## [47,] 1417 1366 1417 1757
## [48,] 1671 1757 1841 2574
## [49,] 1519 1484 1541 1977
## [50,] 1335 1305 1276 2607
## [51,] 1353 1199 1310 1973
## [52,] 1381 1239 1247 2139
## [53,] 1162 1088 975 2443
## [54,] 1211 1081 1003 1417
## [55,] 1226 1139 1056 1727
## [56,] 1023 939 828 2345
## [57,] 1597 1477 1399 2409
## [58,] 1166 1151 988 2799
## [59,] 1770 1691 1769 2253
## [60,] 2411 2337 2231 3767#Valores promedios por cada banda
ms <- aggregate(df, list(datos$CLASE), mean)
ms## Group.1 blue green red NIR
## 1 AGUA 1230.133 1100.2000 970.4667 1629.400
## 2 SUELO 1268.533 1206.1333 1201.5333 2541.133
## 3 URBANO 1433.867 1360.8000 1323.8667 2308.200
## 4 VEGETACION 1025.333 903.2667 673.5333 2580.667#Exportar valores a Excel
#install.packages("writexl")
library(writexl)
#write_xlsx(ms, "C:/Users/ASUS/Documents/TATI 2021-1/ELECTIVA 2/puntos finales.xlsx" )