Práctica 1, asignatura Geomorfología, licenciatura en Geografía, UASD. Estadísticos descriptivos y gráficos de la pendiente y la altura por geomórfonos

Script de R sobre estadísticos descriptivos y gráficos de la pendiente y la altura por geomórfonos (Jasiewicz et al, 2013), para la asignatura «Geomorfología» de la licenciatura en Geografía, UASD. Una copia TXT se encuentra alojada aquí.

Los datos de ejemplo consisten en un modelo digital de elevaciones (DEM) del SRTM (NASA et al, 2004), resolución 90m x 90m, del Bahoruco Oriental, situado al sudoeste de República Dominicana, dentro del cuadrante dado por los siguientes límites (coordenadas UTM, WGS84):

2019565: falso norte
1997785: falso sur
282171: falso este
251751: falso oeste

#CONSOLIDADO DE PRÁCTICAS DE LA ASIGNATURA "GEOMORFOLOGÍA"

#PRÁCTICA 1. ESTADÍSTICOS DESCRIPTIVOS Y GRÁFICOS DE GEOMÓRFONOS, MODELO DIGITAL DE ELEVACIONES Y PENDIENTES CON R. UN MODELO COMPLETO DE LA ESPAÑOLA ESTÁ DISPONIBLE EN: http://geografiafisica.org/r/maestria_geom/practica_01.zip
#IMPORTANTE: EN ESTE SCRIPT, CADA GRÁFICO SE CREA EN UNA VENTANA INDEPENDIENTE. ES NECESARIO CERRAR CADA GRÁFICO CUANDO SE PREVEE QUE NO SE UTILIZARÁ MÁS

#PAQUETES
library(raster) #PARA GENERACIÓN Y MANEJO DE OBJETOS RÁSTER
library(rgdal) #PARA MANEJAR ARCHIVOS DE LA BIBLIOTECA GDAL
library(dplyr) #MANIPULACIÓN DE DATOS EN TABLAS, TABLAS RESUMEN
library(Hmisc) #CARGA EL PAQUETE HARRELL MISCELANEOUS, PARA REALIZAR OPERACIONES VARIADAS. EN ESTE CASO SÓLO SE USARÁ PARA LA FUNCIÓN describe, QUE OBTIENE MEDIDAS DE TENDENCIA CENTAL, DE DISPERSIÓN Y DE POSICIÓN
library(RgoogleMaps) #PARA GENERAR, EN VENTANA GRÁFICA DE R, MAPAS DE GOOGLE
library(ggplot2) #GRÁFICOS ESTILIZADOS

##FIJACIÓN DE DIRECTORIO DE TRABAJO TEMPORAL Y DESCARGA DE DATOS
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/sem_2015_02/geo114/dembo.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

#CARGAR DATOS BASE. SE TRATA DE UN MODELO DIGITAL DE ELEVACIONES (SRTM, 90X90 m) DEL BAHORUCO ORIENTAL, REPÚBLICA DOMINICANA
dem <- raster(readGDAL("dembo.tif")) #ESTE DEM PROVIENE DE OTRO DE ÁMBITO NACIONAL, QUE ABARCA LA ESPAÑOLA COMPLETA, ALOJADO EN http://geografiafisica.org/r/maestria_geom/practica_01.zip. SE RECOMIENDA CORTAR EL DEM EN UN SIG CON INTERFAZ GRÁFICO. DE TODAS FORMAS, R OFRECE LA POSIBILIDAD DE ELEGIR ÁREAS MEDIANTE RECTÁNGULOS CON HERRAMIENTAS GRÁFICAS
pen <- terrain(dem,opt='slope',unit='degrees')

#CON EL DEM, SE GENERÓ IMAGEN DE "GEOMÓRFONOS". ESTO SE LOGRA SUBIENDO EL DEM http://sil.uc.edu/geom/app EL SITIO DEVUELVE EL ARCHIVO (HAY QUE ACTIVAR LOS POP-UPS EN EL NAVEGADOR), EL CUAL ES EL QUE SE CARGA EN LA SENTENCIA SITUADA AQUÍ DEBAJO. EN ESTE SCRIPT, EL ARCHIVO TIF DE GEOMÓRFONOS YA FUE GENERADO Y SE ENCONTRABA EN EL COMPRIMIDO DESCARGADO EN LOS PASOS ANTERIORES, POR LO QUE YA FUE COLOCADO EN LA CARPETA DE TRABAJO Y SÓLO HAY QUE CARGARLO EN R
geom <- raster(readGDAL('geom_dembo.tif.tif'))
geom #NOTAR QUE NO PRESENTA NI NOMBRE DE BANDA (DEBE APARECER EL NOMBRE POR DEFECTO, "band1") TABLA DE ATRIBUTOS
mar <- shapefile("mar.shp")
marr <- rasterize(mar,geom)
marr[marr$layer==1] <- 0
marr[marr$layer!=0] <- NA
geom <- cover(marr,geom)
names(geom) <- "geomorfono" #ASIGNA NOMBRE A LA BANDA
geom<-as.factor(geom) #CONVIERTE EL OBJETO geom A UNA IMAGEN CUYOS VALORES SON DE TIPO FACTORIAL
geom #DEBE PRESENTAR ATRIBUTOS
g<-levels(as.factor(geom))[[1]] #EXTRAE LA CANTIDAD DE NIVELES DEL OBJETO geom Y LOS COLOCA EN UNA TABLA
g$geomorfono<-c("00. sin clasificar", "01. llanura", "02. pico", "03. interfluvio", "04. hombrera", "05. espolón/gajo", "06. vertiente", "07. vaguada", "08. piedemonte", "09. valle", "10. sima")
g
levels(geom)<-g #ESTABLECE LOS NIVELES DEL RASTER DE GEOMÓRFONOS
geom #MUESTRA INFORMACIÓN DEL RASTER DE GEOMÓRFONOS, ACTUALIZANDO LA INFORMACIÓN SOBRE LOS VALORES

#GENERACIÓN DE PALETA DE COLORES DE GEOMÓRFONOS
geompal<-c("#000000","#dcdcdc","#380000","#c80000","#ff5014","#fad23c","#ffff3c","#b4e614","#3cfa96","#0000ff","#000038") #GENERA UNA PALETA DE COLORES PARA CONSEGUIR UN MAPA CON LA COLORACIÓN CONVENCIONAL ESTABLECIDA POR EL AUTOR DE GEOMÓRFONOS

#MAPA DEL RÁSTER DE GEOMÓRFONOS SIN CONFIGURACIÓN DE LEYENDA
dev.new() #GENERA UNA NUEVA VENTANA GRÁFICA
plot(geom,col=geompal) #MUESTRA UN MAPA DE LOS GEOMÓRFONOS, CON LOS COLORES CONVENCIONALES

#MAPA DEL RÁSTER DE GEOMÓRFONOS CON CONFIGURACIÓN DE LEYENDA
dev.new()
par(xpd = FALSE)
plot(geom,col=geompal,legend=F)
par(xpd = TRUE)
legend("right",legend=g$geomorfono,fill=geompal,inset=-0.4)

#CÁLCULO DE ÁREA EN M2 DE CADA GEOMÓRFONO, Y GENERACIÓN DE GRÁFICOS DE BARRAS DE LAS ÁREAS DE CADA GEOMÓRFONO
tfgeom<-table(as.data.frame(geom)$geomorfono)*xres(geom)*yres(geom) #GENERAR UNA TABLA DE FRECUENCIAS DE CADA CLASE DE LOS GEOMÓRFONOS. EN LA SIGUIENTE SENTENCIA SE OBSERVA QUÉ SIGNIFICA CADA NÚMERO
tfgeom
dev.new()
barplot(tfgeom,cex.names=0.9,col=geompal,ylab='área en metros cuadrados',main='Área por tipo de geomórfono')

#SUAVIZACIÓN DEL RÁSTER DE GEOMÓRFÓNOS
mf<-matrix(rep(1,25),nrow=5) #MATRIZ DE FILTRAJE, 5X5, SIN PESOS EN LA PONDERACIÓN
mf #MUESTRA LA MATRIZ
geomf<-focal(geom,mf,modal) #GENERA UN RÁSTER DE GEOMÓRFONOS FILTRADO MEDIANTE LA FUNCIÓN modal. EL ALGORITMO RECUPERA LA CELDA MÁS REPETIDA EN LA VENTANA MÓVIL 5X5, SUAVIZANDO ASÍ EL RESULTADO
dev.new()
par(mfrow=c(1,2)) #LA VENTANA GRÁFICA ADQUIERE LA FORMA DE UN PANEL DE 1 FILA Y 2 COLUMNAS
plot(geom,col=geompal) #MUESTRA UN MAPA DE LOS GEOMÓRFONOS, CON LOS COLORES CONVENCIONALES, EN EL PRIMER PANEL (IZQUIERDA)
plot(geomf,col=geompal) #MUESTRA UN MAPA DE LOS GEOMÓRFONOS FILTRADOS, CON LOS COLORES CONVENCIONALES, EN EL SEGUNDO PANEL (DERECHA)

#RECLASIFICACIÓN, SÓLO VALLES
rcl<-matrix(c(0,9,0,9,10,1,10,11,0),nrow=3,byrow=T) #GENERA LA MATRIZ DE RECLASIFICACIÓN
rcl #MUESTRA LA MATRIZ DE RECLASIFICACIÓN
geomfvalles<-reclassify(geomf,rcl, right=F) #RECLASIFICA EL GEOMÓRFONO FILTRADO
dev.new()
plot(geomfvalles,col=c('white','black'))

#ASIGNANDO NOMBRE A CAPA DE INFORMACIÓN Y COLOCANDO ATRIBUTOS
names(geomfvalles) <- "valles"
geom<-as.factor(geom)
geomfvalles
gvalles<-levels(as.factor(geomfvalles))[[1]] #EXTRAE LA CANTIDAD DE NIVELES DEL OBJETO geomfvalles Y LOS COLOCA EN UNA TABLA
gvalles
gvalles$geomorfono<-c("0. no valles", "1. valles")
gvalles
levels(geomfvalles)<-gvalles #ESTABLECE LOS NIVELES
geomfvalles #MUESTRA INFORMACIÓN DEL RASTER

#MAPA DE VALLES CON LEYENDA
dev.new()
par(xpd = FALSE)
plot(geomfvalles,col=c('white','black'),legend=F)
par(xpd = TRUE)
legend("right",legend=gvalles$geomorfono,fill=c('white','black'),inset=-0.3)

#AGRUPANDO LAS MANCHAS DE VALLES QUE TIENEN CONTINUIDAD TERRITORIAL, ASIGNÁNDOLES UN IDENTIFICADOR ÚNICO
geomfvallesgrp<-clump(geomfvalles)
dev.new()
plot(geomfvallesgrp)

#CONTEO DE CELDAS, ÁREAS
estareavalles <- as.data.frame(zonal(dem,geomfvallesgrp,'count'))
estareavalles$aream2<-estareavalles$count*xres(geomfvallesgrp)*yres(geomfvallesgrp)
estareavalles$areaha<-estareavalles$count*xres(geomfvallesgrp)*yres(geomfvallesgrp)/10000
estareavalles$areakm2<-estareavalles$count*xres(geomfvallesgrp)*yres(geomfvallesgrp)/1000000
head(estareavalles)

#ESTADÍSTICOS DE LA ALTURA
estaltvalles<-as.data.frame(cbind(zone=unique(na.omit(getValues(geomfvallesgrp))),sapply(c('mean','min','max','sd','var','cv'),function(x) zonal(dem,geomfvallesgrp,x)[,2])))
head(estaltvalles)
names(estaltvalles) <- c("zone", paste("alt_",names(estaltvalles)[!names(estaltvalles)=="zone"],sep=""))
head(estaltvalles)

#ESTADÍSTICOS DE LA PENDIENTES
estpenvalles<-as.data.frame(cbind(zone=unique(na.omit(getValues(geomfvallesgrp))),sapply(c('mean','min','max','sd','var','cv'),function(x) zonal(pen,geomfvallesgrp,x)[,2])))
head(estpenvalles)
names(estpenvalles) <- c("zone", paste("pen_",names(estpenvalles)[!names(estpenvalles)=="zone"],sep=""))
head(estpenvalles)

#UNIENDO TABLAS
estaralpevalles <- merge(estareavalles,estaltvalles)
estaralpevalles <- merge(estaralpevalles,estpenvalles)
head(estaralpevalles)

#VALLES DE 2 KM O MÁS
geomfvallesgrp2km <- geomfvallesgrp
geomfvallesgrp2km[geomfvallesgrp2km$clumps %in% (estaralpevalles %>% filter(areakm2<2))$zone] <- NA
dev.new()
par(mfrow=c(1,2))
plot(geomfvallesgrp)
plot(geomfvallesgrp2km)

#NORMALIDAD DE ALTURAS Y PENDIENTES DE VALLES DE 2 KM O MÁS
zonal(dem,geomfvallesgrp2km,function(x,...) shapiro.test(na.omit(x))$p.value)
zonal(pen,geomfvallesgrp2km,function(x,...) shapiro.test(na.omit(x))$p.value)

#QQNORM DE ALTURAS Y PENDIENTES DE VALLES 2 KM O MÁS
dev.new()
par(mfrow=c(3,5))
zonal(dem,geomfvallesgrp2km,function(x,...) qqnorm(x))
dev.new()
par(mfrow=c(3,5))
zonal(pen,geomfvallesgrp2km,function(x,...) qqnorm(x))

#BOXPLOT DE ALTURAS Y PENDIENTES DE VALLES 2 KM O MÁS VALLES DE 2 KM O MÁS
dev.new()
par(mfrow=c(1,2))
boxplot(dem,geomfvallesgrp2km)
boxplot(pen,geomfvallesgrp2km)