En esta actividad se intentĆ³ mostrar las pautas de procesamiento de imĆ”genes digitales, utilizando los cĆ³digos que se ha compartido anteriormente en la siguiente url: https://rpubs.com/jhonathan1060/753925 en Markdown, despuĆ©s de lograr este proceso el ejercicio, la actividad se ha enfocado en generar los perfiles espectrales de los objetos mĆ”s representativos que se encuentran en la imagen, para ello se ubicaron las coordenadas de los objetos de los glaciares, Ć”rea urbana, cuerpos de agua y vegetaciĆ³n, con el fin de generar los perfiles espectrales de los objetos, luego nos basamos en los cĆ³digos publicados por (Robert J. Hijmans,2016) que se ubica en la siguiente url: https://rspatial.org/raster/rs/2-exploration.html, y despuĆ©s en el anĆ”lisis de Normalized Difference Snow Index (NDSI), logrando reclasificar los pixeles que corresponden a cobertura de glaciar. En tal sentido, el propĆ³sito de esta guĆ­a fue mostrar las pautas del proceso de anĆ”lisis de imĆ”genes, respuesta espectral de los objetos y lograr reclasificar las coberturas de glaciares mediante el principio de anĆ”lisis de Ć­ndices, utilizando el programa R Package y RStudio, donde se comparte los cĆ³digos para realizar evaluaciones similares, con las personas que vienen desarrollando habilidades en los programas de ingenierĆ­a, ciencias de la tierra y con la sociedad civil, teniendo en cuenta el alcance local, nacional e internacional.

library(sp)
library(raster)
library(rgeos)
## rgeos version: 0.5-2, (SVN revision 621)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-2 
##  Polygon checking: TRUE
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/Usuario/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/Usuario/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1
# importamos los datos raster 
#C:/Users/Usuario/Desktop/Nueva carpeta (2)
setwd("C:/Users/Usuario/Desktop/Nueva carpeta (2)")
files <- list.files(pattern='.tif')
files
## [1] "b_B2.tif" "b_B3.tif" "b_B4.tif" "b_B5.tif" "b_B6.tif" "b_B7.tif"
b2 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B2.tif")
b3 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B3.tif")
b4 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B4.tif")
b5 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B5.tif")
b6 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B6.tif")
b7 <- raster("C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B7.tif")
b2
## class      : RasterLayer 
## dimensions : 760, 955, 725800  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 220335, 248985, -1056045, -1033245  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/Usuario/Desktop/Nueva carpeta (2)/b_B2.tif 
## names      : b_B2 
## values     : 6669, 48598  (min, max)
plot(b4)

# Conversion de nivel digital a radiancia. 
rad <- function(x2,x3,x4,x5,x6,x7){
  rb2 <- 0.012452*x2+(-62.26080)
  rb3 <- 0.011475*x3+(-57.37280)
  rb4 <- 0.006760*x4+(-48.37997)
  rb5 <- 0.005921*x5+(-29.60613)
  rb6 <- 0.001472*x6+(-7.36278)
  rb7 <- 0.0004963*x7+(-2.48165)
  return(list(rb2,rb3,rb4,rb5,rb6,rb7))
}
rabs<- stack(rad(b2,b3,b4,b5,b6,b7))
rabs$b_B7
## class      : RasterLayer 
## dimensions : 760, 955, 725800  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 220335, 248985, -1056045, -1033245  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : b_B7 
## values     : -0.0378688, 25.27641  (min, max)
#CƔlculo de distancia sol a tierra unidades astronomicas.
dist<- function(j=201){
  da <- (1-0.01672*cos(0.9856*(j-4)*(pi/180)))
  return(da)}# ingresamos el valor del dia juliano en este caso es 217
d<-dist()
d
## [1] 1.016212
#CƔlculo de angulo cenital.
ac <- function(t=47.5053){
  a <- cos((pi/180)*(90-t))
  return(a)
}
z<- ac()
z
## [1] 0.7373398
# calculo de reflectancia 
#para determinar valores de Esun http://www.gisagmaps.com/landsat-8-atco/
ref <- function(){
  rf2<- ((pi)*rabs$b_B2*d^2)/(2067*z)
  rf3<- ((pi)*rabs$b_B3*d^2)/(1893*z)
  rf4<- ((pi)*rabs$b_B4*d^2)/(1603*z)
  rf5<- ((pi)*rabs$b_B5*d^2)/(972.6*z)
  rf6<- ((pi)*rabs$b_B6*d^2)/(245*z)
  rf7<- ((pi)*rabs$b_B7*d^2)/(79.72*z)
  return(list(rf2,rf3,rf4,rf5,rf6,rf7))
}
rfbs<-stack(ref())
rfbs$b_B3
## class      : RasterLayer 
## dimensions : 760, 955, 725800  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 220335, 248985, -1056045, -1033245  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : b_B3 
## values     : 0.006806428, 1.30095  (min, max)
plot(rfbs$b_B2, main = "Azul", col = gray(0:100/100))

plot(rfbs$b_B3, main = "Verde", col = gray(0:100/100))

plot(rfbs$b_B4, main = "Rojo", col = gray(0:100 / 100))

plot(rfbs$b_B5, main = "Infrarrojo I", col = gray(0:100 / 100))

plot(rfbs$b_B6, main = "Infrarojo II", col = gray(0:100 / 100))

library(ggplot2)
library(RStoolbox)
summary(rfbs$b_B2)
##                 b_B2
## Min.    3.196490e-02
## 1st Qu. 7.490512e-02
## Median  9.011975e-02
## 3rd Qu. 1.418070e-01
## Max.    1.311345e+00
## NA's    1.842000e+03
#mostrando los datos raster en el mapa.
plot(rfbs$b_B4)

#Combinacion de bandas.
ggRGB(rfbs,r=3,g=4,b=1,stretch = "lin") + ggtitle("Combinacion bandas 321")

ventana <-extent(240000, 250000,-1050000,-1040000)
ventana
## class      : Extent 
## xmin       : 240000 
## xmax       : 250000 
## ymin       : -1050000 
## ymax       : -1040000
ggRGB(rfbs,r=3,g=4,b=1,stretch = "lin",ext= ventana) + ggtitle("Combinacion bandas 321")

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union
## 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
#Importamos puntos de control
setwd("C:/Users/Usuario/Desktop/Nueva carpeta (2)")
list.files(path=".")
##  [1] "b_B2.tif"    "b_B3.tif"    "b_B4.tif"    "b_B5.tif"    "b_B6.tif"   
##  [6] "b_B7.tif"    "puntos.cpg"  "puntos.dbf"  "puntos.mshp" "puntos.prj" 
## [11] "puntos.shp"  "puntos.shx"
Points <- readOGR('puntos.shp', layer='puntos')
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Usuario\Desktop\Nueva carpeta (2)\puntos.shp", layer: "puntos"
## with 24 features
## It has 2 fields
Points$Name
##  [1] Glaciar     Glaciar     Glaciar     Glaciar     Glaciar     Glaciar    
##  [7] Agua        Agua        Agua        Agua        Agua        Agua       
## [13] Zona urbana Zona urbana Zona urbana Zona urbana Zona urbana Zona urbana
## [19] Vegetacion  Vegetacion  Vegetacion  Vegetacion  Vegetacion  Vegetacion 
## Levels: Agua Glaciar Vegetacion Zona urbana
crs(points)
## [1] NA
proj4string(Points) <-CRS("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ")
crs(Points)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
#
plotRGB(rfbs, r=4, g=3, b=2, axes=FALSE, stretch="lin", main="Landsat","lin")
plot(Points, add=TRUE, legend=FALSE)

# extraemos los valores de pixeles en los puntos de control
df = extract(rfbs, Points)
head(df)
##           b_B2      b_B3      b_B4      b_B5        b_B6        b_B7
## [1,] 0.9452398 0.9532560 0.6394874 0.8139745 0.013300149 0.014180903
## [2,] 0.5859202 0.6049754 0.3694918 0.4973882 0.013458764 0.014180903
## [3,] 0.9660737 0.9634980 0.6310819 0.8164120 0.023583672 0.025356938
## [4,] 0.5419463 0.5377357 0.3315837 0.4546106 0.010497955 0.010893834
## [5,] 0.7511606 0.7537241 0.4932545 0.6523465 0.013855301 0.015550515
## [6,] 0.5348956 0.5278938 0.3204135 0.4463604 0.008198041 0.008620278
#calculamos el ma media de los valores de cada tipo de cobertura
ms = aggregate(df, list(Points$Name), mean)

# mostramos los valores de la media de respuesta espectral de cada tipo de cobertura
rownames(ms) = ms[,1]
ms = ms[,-1]
ms
##                   b_B2       b_B3        b_B4      b_B5        b_B6        b_B7
## Agua        0.08065258 0.05114389 -0.02389529 0.0109192 0.004347227 0.003803808
## Glaciar     0.72087269 0.72351386  0.46421881 0.6135154 0.013815647 0.014797228
## Vegetacion  0.07589469 0.08145641  0.01276344 0.3890111 0.222988856 0.119093194
## Zona urbana 0.14260224 0.14838046  0.07076366 0.1843868 0.213939003 0.207077078
# creamos las colores de los objetos en el perfil especrtal 
mycolor = c('cyan', 'black', 'green',  'purple')
#transform ms from a data.frame to a matrix
ms = as.matrix(ms)
# creamos el plot
plot(0, ylim=c(0,0.9), xlim = c(1,6), type='n', xlab="Bandas", ylab = "Reflectancia")
# adicionamos las clases
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 1, lty = 1, col = mycolor[i])
}

# TĆ­tulo
title(main="Perfil espectral de la imagen Landsat 8", font.main = 2)
# Legend
legend("topright", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =2, bty = "n")

#Considerando los valores de respuesta espectral, podemos intentaridenticar cobertura de nieve  midiante calculo de indices  

index <- function(x,y){
  ndsi <- (x - y ) / (y + x) 
  return(ndsi)
}
NDSI<-index(rfbs$b_B3,rfbs$b_B6)


plot(NDSI, col = rev(terrain.colors(30)), main = 'NDSI con imagenes Landsat 8')

summary(NDSI)
##                 layer
## Min.      -0.64534031
## 1st Qu.   -0.41222847
## Median    -0.32248613
## 3rd Qu.   -0.03956041
## Max.       1.65998265
## NA's    1843.00000000
# reclasificando en base a umbrales 
glac <- reclassify(NDSI, c(-0.65,0,1, 0,0.5,2, 0.5,0.6,3, 0.6,Inf,4))
plot(glac,col = rev(terrain.colors(5)), main = 'NDVI base a umbrales')

#Identificando grupos de pixeles que no correponde a glaciares 
glac <- reclassify(NDSI, c(-0.65,0,NA, 0,0.5,NA, 0.5,0.9,NA, 0.9,Inf,2))
plot(glac,col = rev(terrain.colors(5)), main = 'NDVI base a umbrales')

# 
par(mfrow = c(2,2))
plotRGB(rfbs, r=4, g=3, b=2, axes=FALSE, stretch="lin", main="Landsat","lin",ext= ventana)
plotRGB(rfbs, r=4, g=3, b=2, axes=FALSE, stretch="lin", main="Landsat","lin",ext= ventana)
plot(glac, add=TRUE, legend=FALSE)
#

par(mfrow=c(1,1))

plotRGB(rfbs, r=4, g=3, b=2, axes=TRUE, stretch="lin", main="Cobertura Glaciar","lin")
plot(glac, add=TRUE, legend=FALSE)