# install.packages("raster")
# install.packages("sp")
# install.packages("rgdal")
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/NICOLAS/Documents/R/win-library/4.0/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/NICOLAS/Documents/R/win-library/4.0/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/NICOLAS/Documents/R/win-library/4.0/rgdal/proj
setwd("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3")
# Llamar imagines satelitales
# Azul banda 2
b2 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B02.jp2")
# Verde banda 3
b3 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B03.jp2")
# Rojo banda 4
b4 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B04.jp2")
b8 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B08.jp2")
# Informacion de la imagen
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)
crs(b2) # consultar el sistema de referencia
## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
ncell(b2) # Consultar numero de celdas
## [1] 120560400
dim(b2) # Dimensiones, es decir alto x ancho
## [1] 10980 10980 1
res(b2) # Resolucion
## [1] 10 10
nlayers(b2) # Numero de capas
## [1] 1
compareRaster(b2,b3) # Comparar 2 imagenes
## [1] TRUE
s <- stack(b8,b4,b3)
s
## class : RasterStack
## dimensions : 10980, 10980, 120560400, 3 (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 : T18NVH_20170921T152631_B08, T18NVH_20170921T152631_B04, T18NVH_20170921T152631_B03
## min values : 0, 0, 0
## max values : 65535, 65535, 65535
par(mfrow = c(2,2))
plot(b2, main = "azul", col = gray(0:100 / 100))
plot(b3, main = "verde", col = gray(0:100 / 100))
plot(b4, main = "rojo", col = gray(0:100 / 100))
plot(b8, main = "infra rojo", col = gray(0:100 / 100))
# Combinacion para lograr el color verdadero
landsatRGB <- stack(b4,b3,b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat color verdadero")
filenames <- paste0("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B0", 2:8,".jp2")
filenames
## [1] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B02.jp2"
## [2] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B03.jp2"
## [3] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B04.jp2"
## [4] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B05.jp2"
## [5] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B06.jp2"
## [6] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B07.jp2"
## [7] "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/T18NVH_20170921T152631_B08.jp2"
landsat <- stack(b4,b3,b2,b8)
landsat
## 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 : T18NVH_20170921T152631_B04, T18NVH_20170921T152631_B03, T18NVH_20170921T152631_B02, T18NVH_20170921T152631_B08
## min values : 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535
names(landsat)
## [1] "T18NVH_20170921T152631_B04" "T18NVH_20170921T152631_B03"
## [3] "T18NVH_20170921T152631_B02" "T18NVH_20170921T152631_B08"
# usando la funcion extent
extent(landsat)
## class : Extent
## xmin : 399960
## xmax : 509760
## ymin : 190200
## ymax : 3e+05
# limitar el recorte
recorte <- extent(399960, 469000, 238000, 3e+05)
# se realiza el recorte
landsatcortada <- crop(landsat, recorte)
plot(landsatcortada)
plotRGB(landsatcortada, axes = TRUE, stretch = "lin", main = "Landsatcortada color verdadero")
# REnombre de bandas
names(landsatcortada) <- c("rojo","verde", "azul", "NIR")
#se usa otro metodo llamando un csv
# install.packages("dplyr")
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
puntos_de_control <- read.csv("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/Puntos de control 2.csv")
# usamos el paquete para manejar datos geograficos
# install.packages("sf")
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
puntos_de_control <- puntos_de_control %>%
st_as_sf(coords = c("x", "y"))
plot(puntos_de_control$geometry)
df <- extract(landsatcortada, puntos_de_control); df
## rojo verde azul NIR
## [1,] 658 902 1129 513
## [2,] 666 896 1124 586
## [3,] 816 980 1158 659
## [4,] 1841 1618 1602 2453
## [5,] 692 930 1147 580
## [6,] 925 1025 1193 905
## [7,] 936 1017 1193 914
## [8,] 911 998 1205 884
## [9,] 880 1016 1159 1667
## [10,] 941 1037 1179 1121
## [11,] 861 960 1061 2092
## [12,] 1145 1124 1220 1992
## [13,] 1119 1045 1155 2142
## [14,] 1180 1120 1199 1930
## [15,] 1871 1528 1496 2733
## [16,] 988 977 1174 1282
## [17,] 1090 1064 1168 2139
## [18,] 1686 1536 1568 2205
## [19,] 1313 1247 1362 2430
## [20,] 1371 1499 1531 3091
## [21,] 737 1100 1088 3974
## [22,] 965 1302 1303 1358
## [23,] 1581 1378 1456 2119
## [24,] 1720 1406 1420 2256
## [25,] 1739 1435 1401 3236
## [26,] 1060 1078 1220 1858
## [27,] 733 999 1127 3120
## [28,] 1213 1212 1304 2434
## [29,] 1908 1492 1498 3178
## [30,] 1167 1299 1340 3321
## [31,] 1278 1384 1520 2149
## [32,] 1438 1460 1573 2097
## [33,] 1807 1767 1832 2543
## [34,] 1141 1211 1267 2036
## [35,] 1523 1209 1240 2398
## [36,] 1302 1315 1319 2481
## [37,] 1950 1783 1671 2618
## [38,] 2597 2273 2311 3153
## [39,] 1350 1434 1495 1965
## [40,] 1405 1376 1605 2070
ms <- aggregate(df, list(puntos_de_control$Clase), mean)
ms
## Group.1 rojo verde azul NIR
## 1 Agua 926.6 1041.9 1208.9 1028.2
## 2 Construccion 1579.1 1521.2 1583.3 2351.0
## 3 Suelo 1262.4 1210.0 1293.4 2203.6
## 4 Vegetacion 1282.3 1270.1 1315.7 2685.4
# install.packages("writexl")
library(writexl)
# write_xlsx(ms, "D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 3/Ultimo excel.xlsx")