#CARGA DE PAQUETES
library(raster) #PARA GENERACIÓN Y MANEJO DE OBJETOS RÁSTER
library(rgdal) #PARA MANEJAR ARCHIVOS DE LA BIBLIOTECA GDAL
library(maptools)
 
#FIJACIÓN DE DIRECTORIO DE TRABAJO TEMPORAL Y DESCARGA DE DATOS
#PARA FACILITAR LA LECTURA Y EL PROCESAMIENTO DE LOS VALORES DE ALTURA, ESTE SCRIPT CONSIDERA UN DEM CON VALOR MÍNIMO 0. DE ESTA MANERA SE EVITAN LOS VALORES NEGATIVOS QUE SE OBTIENEN EN MUNICIPIOS COSTEROS. TIENE COMO DESVENTAJA QUE LOS MUNICIPIOS POR DEBAJO DEL NIVEL DEL MAR NO REFLEJARÁN VALORES NEGATIVOS. DADO QUE LA FINALIDAD ES ESTABLECER CORRELACIONES DE LA ALTURA RESPECTO DEL VECTOR DEL ZIKA Y OTRAS ENFERMEDADES, SU IMPACTO ES MÍNIMO. POR LO TANTO, LOS VALORES NEGATIVOS NO SON CONSIDERADOS. PARA GENERACIÓN DE UN DEM CON VALORES NEGATIVOS, CONSULTAR SCRIPT SEPARADO EN LA RUTA (DESCARGA DOS DEM, LOS UNE Y LOS ENMASCARA EN EL ÁMBITO DE RD): https://www.geografiafisica.org/2015/09/28/r-generacion-modelo-digital-de-elevaciones-de-rd-uniendo-dos-imagenes-separadas/
setwd(tempdir()) #FIJAR EL DIRECTORIO DE TRABAJO
getwd() #MUESTRA DÓNDE ESTÁ EL DIRECTORIO DE TRABAJO
temp <- tempfile() #ASIGNAR EL NOMBRE DE UN ARCHIVO TEMPORAL AL OBJETO temp
download.file("http://geografiafisica.org/r/maestria_geom/practica_01.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 DE ELLOS. TIENEN LA VENTAJA DE QUE OFRECEN VALORES BRUTOS, POR LO QUE SE TRABAJA CON ALTURAS "VERDADERAS", AUNQUE ALGUNOS ERRORES SE HAN VERIFICADO
 
#LECTURA MUNICIPIOS DE LA ONE
temp <- tempfile() #ASIGNAR EL NOMBRE DE UN ARCHIVO TEMPORAL AL OBJETO temp
download.file("http://geografiafisica.org/r/divisionrd/municipios.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
 
#LECTURA DEL DEM DEL SRTM
dem<-raster(readGDAL("demesp.tif")) #LEE EL TIF DEL DEM DE LA ESPAÑOLA CON LA FUNCIÓN readGDAL DEL PAQUETE rgdal, EL CUAL SE ENCONTRABA EN EL COMPRIMIDO ZIP, Y LO ASIGNA AL OBJETO dem
dem #MUESTRA LAS CARACTERÍSTICAS DEL raster DENOMINADO dem
dev.new();dev.new() #CREA UNA VENTANA GRÁFICA NUEVA
plot(dem) #VISUALIZA EL RASTER dem EN UN MAPA
 
#LECTURA DE MUNICIPIOS
mu <- shapefile('MUNCenso2010.shp')
plot(mu,add=T) #VISUALIZA LOS MUNICIPIOS SOBRE EL RASTER dem
 
#HAY TRES ALTERNATIVAS PARA GENERAR LOS ESTADÍSTICOS DESCRIPTIVOS:
##1) ESTADÍSTICOS DESCRIPTIVOS DE LA ALTURA POR MUNICIPIO MEDIANTE zonal A RESOLUCIÓN NATIVA. EN ESTE CASO, SE USA LA RESOLUCIÓN ORIGINAL DEL DEM, SIN REMUESTREO, PARA GENERAR LOS ESTADÍSTICOS. UNA ALTERNATIVA SERÍA HACER REMUESTREO Y TRABAJAR CON MENOS CELDAS, A SABIENDAS DE LA PÉRDIDA DE PRECISIÓN. PARA ESTUDIOS DE ESCALA NACIONAL, UN REMUESTREO NO TENDRÍA IMPACTO SIGNIFICATIVO
mu$ENLACE2 <- as.numeric(paste('999', mu$ENLACE,sep='')) #CREACIÓN DE UN CAMPO NUMÉRICO DENOMINADO 'ENLACE2' QUE SERÁ USADO COMO ATRIBUTO DEL RASTER A CREAR.
rmu <- rasterize(mu,dem,'ENLACE2') #RASTERIZACIÓN DEL SPATIALPOLYGONSDATAFRAME mu. PODRÍA TARDAR
edh<-as.data.frame(cbind(id=sort(unique(na.omit(getValues(rmu)))),sapply(c('mean','median','min','max','sd','var','cv','count'),function(x) zonal(dem,rmu,x)[,2]))) #GENERACIÓN DE UNA TABLA DE ESTADÍSTICOS DESCRIPTIVOS MEDIANTE UNA ÚNICA SENTENCIA. ESTA ES LA PARTE EFICIENTE DE ESTA ALTERNATIVA, DADO QUE MEDIANTE EXTRACT SE VUELVE MÁS LENTO
edh$ENLACE <- substr(edh$id,4,100) #CREACIÓN DE UN CAMPO DE ENLACE PARA UNIR ESTOS ESTADÍSTICOS A LOS ATRIBUTOS DEL SPATIALPOLYGONSDATAFRAME
###ANEXANDO LOS ESTADÍSTICOS DESCRIPTIVOS GENERADOS AL OBJETO SP MEDIANTE spCbind
o <- match(mu$ENLACE,edh$ENLACE) #OBTENER UN CAMPO DE ORDENACIÓN
edho <- edh[o,] #CON EL CAMPO DE ORDENACIÓN, GENERAR UN NUEVO OBJETO DE ESTADÍSTICOS DESCRIPTIVOS DEBIDAMENTE ORDENADO
rownames(edho) <- rownames(slot(mu,'data')) #ASIGNAR LOS MISMO NOMBRES DE FILA AL OBJETO DE ESTADÍSTICOS DESCRIPTIVOS, TOMÁNDOLO DE LOS NOMBRES DE FILA DEL OBJETO SP
match(mu$ENLACE, edho$ENLACE) #CONFIRMACIÓN DE QUE EL ORDEN ESTÁ SINCRONIZADO ENTRE AMBAS FUENTES
edhsp <- spCbind(mu, edho) #UNIÓN Y GENERACIÓN DE NUEVO OBJETO DE ESTADÍSTICOS DESCRIPTIVOS DE TIPO SPATIALPOLYGONSDATAFRAME
edhsp <- edhsp[,!(names(edhsp) %in% c('ENLACE2','id','ENLACE.1'))]
slot(edhsp,'data')
dev.new();spplot(edhsp['mean'], col.regions = grey(seq(.9,.5,length.out=17)), col='grey')
dev.new();spplot(edhsp[c('mean','median','min','max')], col.regions = grey(seq(.9,.5,length.out=17)), col='grey')
###EXPORTAR A ARCHIVO DE VALORES SEPARADOS POR COMA
write.csv(slot(edhsp,'data'),file='tabla_estadisticos_descriptivos.csv',row.names=F) 
###COMO DATA FRAME
edht <- merge(slot(mu,'data'),edh)
write.csv(edht,file='tabla_estadisticos_descriptivos.csv',row.names=F) 
 
##2)ÍDEM ANTERIOR, PERO CON UN DEM REMUESTREADO PARA AGILIZAR LOS PROCEDIMIENTOS. EN EL EJEMPLO, SE USARÁ UN FACTOR DE 4, ES DECIR, LA RESOLUCIÓN DEL DEM ORIGINAL, QUE ERA DE CA. 90, QUEDARÁ EN CA. 360 METROS
dem4 <- aggregate(dem,4,fun=mean,na.rm=T) #REMUESTREO AGREGANDO
mu$ENLACE2 <- as.numeric(paste('999', mu$ENLACE,sep='')) #A PARTIR DE ESTA SENTENCIA SE REPITEN LOS MISMOS PASOS QUE EN LA ALTERNATIVA 1
rmu4 <- rasterize(mu,dem4,'ENLACE2')
edh4<-as.data.frame(cbind(id=sort(unique(na.omit(getValues(rmu4)))),sapply(c('mean','median','min','max','sd','var','cv','count'),function(x) zonal(dem4,rmu4,x)[,2])))
edh4$ENLACE <- substr(edh4$id,4,100)
o <- match(mu$ENLACE,edh4$ENLACE) #OBTENER UN CAMPO DE ORDENACIÓN
edho4 <- edh4[o,]
rownames(edho4) <- rownames(slot(mu,'data'))
match(mu$ENLACE, edho4$ENLACE)
edhsp4 <- spCbind(mu, edho4)
edhsp4 <- edhsp4[,!(names(edhsp4) %in% c('ENLACE2','id','ENLACE.1'))]
slot(edhsp4,'data')
dev.new();spplot(edhsp4['mean'], col.regions = grey(seq(.9,.5,length.out=17)), col='grey')
dev.new();spplot(edhsp4[c('mean','median','min','max')], col.regions = grey(seq(.9,.5,length.out=17)), col='grey')
write.csv(slot(edhsp,'data'),file='tabla_estadisticos_descriptivos.csv',row.names=F) 
 
##NOTA. UNA COMPARACIÓN MEDIANTE PRUEBA ESTADÍSTICA APAREADA REVELA QUE LOS VALORES MEDIOS DE ALTURA OBTENIDOS NO SON SIGNIFICATIVAMWNTE DIFERENTES
comparar <- read.csv('http://geografiafisica.org/r/divisionrd/contraste_medias.csv')
wilcox.test(log(comparar$mean),log(comparar$mean4), paired=T)
 
##3)ESTADÍSTICOS DESCRIPTIVOS DE LA ALTURA POR MUNICIPIO MEDIANTE extract
mediah <- extract(dem,mu, na.rm=T, fun=mean, df=T, sp=T) #CALCULA LA MEDIA DE LA ALTURA PARA CADA MUNICIPIO. PODRÍA TARDAR
colnames(mediah@data)[which(colnames(mediah@data)=='band1')] <- 'mediah' #CAMBIA EL NOMBRE A LA COLUMNA
medianah <- extract(dem,mu, na.rm=T, fun=median, df=T, sp=T) #CALCULA LA MEDIANA DE LA ALTURA PARA CADA MUNICIPIO. PODRÍA TARDAR
colnames(medianah@data)[which(colnames(medianah@data)=='band1')] <- 'medianah' #CAMBIA EL NOMBRE A LA COLUMNA
minh <- extract(dem,mu, na.rm=T, fun=min, df=T, sp=T) #CALCULA EL VALOR MÍNIMO DE LA ALTURA PARA CADA MUNICIPIO. PODRÍA TARDAR
colnames(minh@data)[which(colnames(minh@data)=='band1')] <- 'minh' #CAMBIA EL NOMBRE A LA COLUMNA
maxh  <- extract(dem,mu, na.rm=T, fun=max, df=T, sp=T) #CALCULA EL VALOR MÁXIMO DE LA ALTURA PARA CADA MUNICIPIO. PODRÍA TARDAR
colnames(maxh@data)[which(colnames(maxh@data)=='band1')] <- 'maxh' #CAMBIA EL NOMBRE A LA COLUMNA
###UNIR RESULTADOS
t1 <- merge(as.data.frame(mediah), as.data.frame(medianah))
t2 <- merge(as.data.frame(minh), as.data.frame(maxh))
t3 <- merge(t1, t2)
###EXPORTAR A ARCHIVO DE VALORES SEPARADOS POR COMA
write.csv(t3,file='tabla_estadisticos_descriptivos.csv',row.names=F)