#CREACIÓN DE MÁSCARA DE NUBES PARA IMAGEN LANDSAT 7 ETM+ CON R

#REFERENCIAS:
##GOSLEE (2011): http://www.jstatsoft.org/article/view/v043i04/v43i04.pdf
##GUERRERO (2011): https://joseguerreroa.wordpress.com/2011/10/11/mascara-de-nubes-usando-lenguaje-r-y-qgis/

#IMAGEN FUENTE: ISLA DE BIOKO, GUINEA ECUATORIAL:
##ID: LE71870582003065SGS00
##DESCARGADA DESDE: http://earthexplorer.usgs.gov/
##NO ES NECESARIO DESCARGARLA DESDE EarthExplorer. SE INCLUYE EN EL SCRIPT EL VÍNCULO DE DESCARGA DE UN ZIP QUE CONTIENE SÓLO LAS BANDAS NECESARIAS. DE ESTA FORMA, SE REALIZA EL EJEMPLO DENTRO DEL PROPIO REALIZARLO DESDE EL PROPIO R

#CARGA DE PAQUETES
library(landsat)
library(rgdal)
library(raster)

#FIJACIÓN DE DIRECTORIO DE TRABAJO TEMPORAL Y DESCARGA DE DATOS
setwd(tempdir()) #FIJAR EL DIRECTORIO DE TRABAJO
temp <- tempfile() #ASIGNAR EL NOMBRE DE UN ARCHIVO TEMPORAL AL OBJETO temp
download.file("http://geografiafisica.org/r/nubes_etm.zip",temp) #DESCARGA EL ARCHIVO ZIP Y LO NOMBRA COMO TEMPORAL
unzip(temp) #DESCOMPRIME EN EL DIRECTORIO DE TRABAJO
unlink(temp) #DESVINCULA DEL ARCHIVO TEMPORAL Y LO BORRA
##NOTA: TAMBIÉN ES POSIBLE CONSEGUIR UN DEM POR CUADROS (TILES), MEDIANTE LA FUNCIÓN getData: getData('SRTM',lat=19,lon=-71). TIENE UN INCONVENIENTE: OBLIGA A DESCARGAR ARCHIVOS MUY GRANDES Y NO ESTÁN POR PAÍS, SINO POR TILES, Y RD SE CONSTRUYE CONCATENANDO DOS

#IMPORTANDO CAPAS A R
a03b1 <- readGDAL('LE71870582003065SGS00_B1.TIF')
a03b61 <- readGDAL('LE71870582003065SGS00_B6_VCID_1.TIF')
a03b62 <- readGDAL('LE71870582003065SGS00_B6_VCID_2.TIF')

#CONVERSIÓN A RASTER PARA REALIZAR CORTE
ra03b1 <- raster(a03b1)
ra03b61 <- raster(a03b61)
ra03b62 <- raster(a03b62)

#CORTE EN ÁMBITO DE INTERÉS
ra03b1ib <- crop(ra03b1,extent(430000,490000,350000,420000))
ra03b61ib <- crop(ra03b61,extent(430000,490000,350000,420000))
ra03b62ib <- crop(ra03b62,extent(430000,490000,350000,420000))

#QUITANDO OBJETOS GRANDES
rm(list=ls()[grepl('^a03',ls())])
rm(list=ls()[!grepl('ib',ls())])
ls()

#GRÁFICO EN PANEL
dev.new()
par(mfrow=c(1,3))
plot(ra03b1ib, col=gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))
plot(ra03b61ib, col=gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))
plot(ra03b62ib, col=gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))

#CONVERSIÓN DE RASTER A SpatialGridDataFrame
a03b1ib <- as(ra03b1ib,'SpatialGridDataFrame')
a03b61ib <- as(ra03b61ib,'SpatialGridDataFrame')
a03b62ib <- as(ra03b62ib,'SpatialGridDataFrame')

#EXPLORANDO LOS ESTADÍSTICOS DE CADA BANDA
summary(a03b1ib)
summary(a03b61ib)
summary(a03b62ib)

#CORRECCIÓN RADIOMÉTRICA
#DOS OPCIONES
#POR REFLECTANCIA APARENTE, PERO RESULTAN VALORES NEGATIVOS. LOS VALORES DE Grescale Y Brescale SE OBTIENEN DEL ARCHIVO DE METADATOS, Y TAMBIÉN PUEDEN CONSULTARSE EN CHANDER (2009) (VER: http://www.sciencedirect.com/science/article/pii/S0034425709000169) SÓLO PARA IMÁGENES DE SENSORES MSS, TM, ETM+ Y EO ALI
a03b1ibra <- radiocorr(a03b1ib, Grescale=0.779, Brescale=-6.97874, sunelev=56.67275067, edist=ESdist('2003-03-06'), Esun=1977, method = "apparentreflectance")
#POR MÉTODO DOS USANDO HAZE. BASADO EN LA AYUDA DE radiocorr (?radiocorr) (DISPONIBLE TAMBIÉN EN: http://www.inside-r.org/packages/cran/landsat/docs/radiocorr)
SHV <- table(a03b1ib@data) 
SHV <- min(as.numeric(names(SHV)[SHV > 1000]))
a03b1ibrdos <- radiocorr(a03b1ib, Grescale=0.779, Brescale=-6.97874, sunelev=56.67275067, edist=ESdist('2003-03-06'), Esun=1997, Lhaze=-4, method = "DOS")

#TÉRMICAS
a03b61ibt <- thermalband(a03b61ib, band = 61) #DA ERROR EN DETERMINADOS PIXELES, PERO AUN ASÍ GENERA LA CAPA
a03b62ibt <- thermalband(a03b62ib, band = 62)

#EN LA FUNCIÓN clouds() NO SE GENERA UNA MÁSCARA ADECUADA, Y TODOS LOS PIXELES SE CLASIFICAN COMO NA, PORQUE EL UMBRAL POR DEFECTO (level=0.0014) NO ES EL ADECUADO QUIZÁ POR LA ALTA CANTIDAD DE BRUMA. POR TAL RAZÓN, writeGDAL NO GENERA OBJETO Y PIDE QUE LE INTRODUZCAN UN OBJETO CON BANDA NUMÉRICA
writeGDAL(clouds(a03b1ibra,a03b61ibt), 'nubes01.tif', drivername = 'GTiff')
writeGDAL(clouds(a03b1ibra,a03b62ibt), 'nubes02.tif', drivername = 'GTiff')
writeGDAL(clouds(a03b1ibrdos,a03b61ibt), 'nubes03.tif', drivername = 'GTiff')
writeGDAL(clouds(a03b1ibrdos,a03b62ibt), 'nubes04.tif', drivername = 'GTiff')

#EDITANDO EL UMBRAL Y EL BUFFER POR DEFECTO, writeGDAL SÍ FUNCIONA. SE PRUEBAN VARIAS COMBINACIONES: LAS IMÁGENES CORREGIDAS RADIOMÉTRICAMENTE A PARTIR DE LA BANDA 1 USANDO LOS MÉTODOS REFLECTANCIA APARENTE Y DOS; LAS IMÁGENES TÉRMICAS GENERADAS A PARTIR DE LAS BANDAS 61 Y LA 62. TAMBIÉN SE PUEDE JUGAR CON EL UMBRAL, PERO AL PARECER 0.0007 ES APROPIADO
#PARA OTROS CASOS, SE RECOMIENDA PROBAR DISTINTOS VALORES DE UMBRAL (level) Y DE BUFFER (buffer). EVALUAR LOS RESULTADOS VISUALMENTE EN QGIS
writeGDAL(clouds(a03b1ibra,a03b61ibt,level=0.0007, buffer=5),'nubes05.tif',drivername = 'GTiff') #CONSERVADOR, ALGUNAS NUBES NO SON ENMASCARADAS
writeGDAL(clouds(a03b1ibra,a03b62ibt,level=0.0007, buffer=5),'nubes06.tif',drivername = 'GTiff')
writeGDAL(clouds(a03b1ibrdos,a03b61ibt,level=0.0007, buffer=5),'nubes07.tif',drivername = 'GTiff')
writeGDAL(clouds(a03b1ibrdos,a03b62ibt,level=0.0007, buffer=5),'nubes08.tif',drivername = 'GTiff')

#LOCALIZAR RUTA DE TRABAJO Y ABRIR LOS ARCHIVOS nubes05.tif, nubes06.tif, nubes07.tif Y nubes08.tif DESDE QGIS
getwd()