Estadísticos descriptivos de la altura por municipios de República Dominicana

Script de R para calcular los estadísticos descriptivos de la altura por municipios de República Dominicana. La delimitación seguida es la de ONE, y el DEM proviene del SRTM, descargado desde EarthExplorer. Se incluye igualmente la tabla de resultados.

El resultado es una tabla que puede descargarse desde aquí. A continuación, una muestra de la tabla:

Muestra de la tabla de estadísticos descriptivos de la altura por municipios de República Dominicana

Muestra de la tabla de estadísticos descriptivos de la altura por municipios de República Dominicana

Los productos cartográficos que se pueden generar con estos estadísticos descriptivos son múltiples. A continuación, 4 mapas: media, mediana, máximo y mínimo.

Mapas de la media, mediana, máximo y mínimo de alturas por municipios de República Dominicana

Mapas de la media, mediana, máximo y mínimo de alturas por municipios de República Dominicana

A propósito del zika, estos estadísticos constituyen fuentes relevantes en la modelización de la distribución de su vector.

Nota: la tabla resultante es incoherente para municipios contiguos que comparten una misma cota máxima según el instrumento legal que los crea. Por ejemplo, Bohechío y San José de las Matas comparten los picos Duarte y La Pelona. Según el SRTM, el primero tiene una altura máxima de 3089 metros mientras que el segundo alcanza lo 3095 m (al final va a ser que La Pelona es más alta que el Duarte). Como se observa en el siguiente mapa, las máximas elevaciones del pico Duarte se inscriben en Bohechío, y las de La Pelona en San José de las Matas, obteniéndose como cotas máximas para cada uno los valores ya mencionados. Los límites del polígono, por lo tanto, no coinciden con los de las máximas alturas obtenidas por medio del SRTM. El resultado esperado tendría que ser que ambos municipios tuvieran la misma cota máxima (3095 m), pero el problema de alineación no lo permite.

Límites de los municipios Bohechío y San José de las Matas, con indicación de las alturas máximas en los picos Duarte y La Pelona

Límites de los municipios Bohechío y San José de las Matas, con indicación de las alturas máximas en los picos Duarte y La Pelona

#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)

Created by Pretty R at inside-R.org

Dr. José Ramón Martínez Batlle (Ph.D)

2 pensamientos en “Estadísticos descriptivos de la altura por municipios de República Dominicana

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *