Maestría en Hidrología

En esta lección aprenderás a:
- Crear gráficos con ggplot
- Cuáles son las funciones básicas ggplot
- Manipular argumentos de funciones gráficas
- Crear gráficos en grillas de clasificación
- Crear un mapa y sus componentes
- Crear animaciones a partir de conjunto de mapas

Temario

A- Gráfico barplot
B- Gráfico boxplot
C- Unificar dos gráficos
D- Uso de facet wrap
E- Scatter plot
G- Gráfico de un mapa
H- Animaciones en R

Para trabajar en esta práctica necesitamos dos archivos de texto (resultados de anásilis previos) y por otra parte para la generación de mapas y animación descargar el archivo de texto (listado de estaciones y coordenadas) y archivo .tif (GEotiff). Los cuatro archivos se encuentran en este enlace.

Es posible que la práctica de mapas pueda efectuarse con algún archivo disponible del módulo de Geoestadística, en este caso Ud. necesita un archivo .tif (brick de mapas) y un archivo de texto de las coordenadas de las estaciones. De hecho, este ejercicio se sugiere para trabajar con el código proporcionado y modificarlo de acuerdo a las necesidades.

Gráfico de Barras

Dado el data frame df_events_results vamos a crear una tabla de conteo de los eventos con criterio de separación de año y vamos a convertir este resultado nuevamente a un data frame

library(ggplot2)
library(gridExtra)
library(raster)

## Loading required package: sp

df_events_results <- read.table("df_events_results.txt", sep = "\t")
table_events_by_year <- as.data.frame(table(df_events_results$year))

Ahora vamos a crear un objeto ggplot y a agregar elementos al resultado (imagine que estamos agregando capas a una plantilla en blanco).

graph_events_by_year <- ggplot(table_events_by_year, aes(x = Var1, y = Freq, fill = Var1)) +
  geom_bar(stat='identity',drop=FALSE) + xlab("Year") + ylab("Freq") +
  theme(legend.title=element_blank(),axis.text.x = element_text(angle = 20, hjust = 1))

## Warning: Ignoring unknown parameters: drop

ggsave(file='events_by_year.pdf',graph_events_by_year)

## Saving 7 x 5 in image

Eventos por año

Gráfico boxplot

pdf('events_duration_boxplot_group.pdf')
ggplot() + geom_boxplot(data=df_events_results, mapping=aes(x=month, y=duration, color=month))
ggplot() + geom_boxplot(data=df_events_results[df_events_results$duration < 500,], mapping=aes(x=month, y=duration, color=month))
ggplot() + geom_boxplot(data=df_events_results[df_events_results$duration < 100,], mapping=aes(x=month, y=duration, color=month))
dev.off()

## quartz_off_screen 
##                 2

Duración de eventos

Unificar dos gráficos

boxplot_a_SYNOP <- ggplot() +  #SYNOP
    geom_boxplot(data=df_events_results, mapping=aes(x=SYNOP, y=log(a), color=SYNOP)) +
    theme(legend.title=element_blank(),axis.text.x = element_text(angle = 20, hjust = 1))

boxplot_b_SYNOP <- ggplot() +
    geom_boxplot(data=df_events_results, mapping=aes(x=SYNOP, y=b, color=SYNOP)) +
    theme(legend.title=element_blank(),axis.text.x = element_text(angle = 20, hjust = 1))

join_boxplot_a_b_SYNOP <- arrangeGrob(boxplot_a_SYNOP,boxplot_b_SYNOP,ncol=1)
ggsave(file='boxplot_a_b_SYNOP.pdf',join_boxplot_a_b_SYNOP)

## Saving 7 x 5 in image

Boxplot a y b por SYNOP

Uso de facet

El comando facet_wrap permite geenrar varios gráficos del mismo tipo en base a una clasificación (subsetting) de los datos en un solo entorno gráfico.

pdf('events_duration_histogram_group.pdf')
ggplot(df_events_results, aes(x = duration, color=month)) + geom_histogram() + facet_wrap(~month) + theme(legend.position="none")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(df_events_results[df_events_results$duration < 500,], aes(x = duration, color=month)) + geom_histogram() + facet_wrap(~month) + theme(legend.position="none")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(df_events_results[df_events_results$duration < 100,], aes(x = duration, color=month)) + geom_histogram() + facet_wrap(~month) + theme(legend.position="none")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

dev.off()

## quartz_off_screen 
##                 2

Histogramas de duración

Scatter plot

Vamos a hacer uso de la incorporación de themes para mejorar la calidad del gráfico.

compile_df_events_results_details <- read.table("compile_df_events_results_details.txt", sep = "\t")
compile_df_events_results <- compile_df_events_results_details[compile_df_events_results_details$site=="ECSF",]

scatter_distrib_a_b_site_month_bw <- ggplot(compile_df_events_results_details, 
                                            aes(x = log(a),y = log(b), shape=site, colour=site, fill=site)) + 
  scale_shape_discrete(solid=F) +
  geom_jitter(size=2,alpha=0.8) +
  scale_colour_brewer(palette="Dark2") + scale_fill_brewer(palette="Dark2") +
  facet_wrap(~SYNOP,ncol=3) + 
  xlab(expression(paste("Coefficient ",italic("A"),", (log scale)"))) +
  ylab(expression(paste("Exponent ",italic("b"),", (log scale)"))) +
  theme_classic() + 
  theme(strip.background =  element_rect(fill = "#F1F1F1", colour = "#F1F1F1"),
        strip.text.x = element_text(colour = "black", angle=0, size = 12,
                                    hjust = 0.5, vjust = 0.5),
        legend.position=c(0.93,0.05), legend.title=element_blank(),
        legend.background=element_blank(),
        panel.margin.y= unit(1, "lines"),
        legend.key.size = unit(0.4, "cm"),
        axis.text = element_text(size=12),
        panel.margin.x = unit(1, "cm"))

## Warning: `panel.margin.x` is deprecated. Please use `panel.spacing.x`
## property instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

ggsave(file='scatter_distrib_a_b_site_month_colors.pdf',scatter_distrib_a_b_site_month_bw )

## Saving 7 x 5 in image

## Warning in log(b): NaNs produced

## Warning in log(b): NaNs produced

## Warning: Removed 12 rows containing missing values (geom_point).

Scatter a y b

Gráfico de un mapa

Abrimos el raster de varias capas y accedemos a la primera de ellas. A continuación generamos un data frame en base a las coordenadas y los valores para poder acer uso de la librería ggplot.

library(ggplot2)
library(gridExtra)
library(raster)
r <- raster('temp_max_dialy_data.tif', band=1)
val <- getValues(r)
xy <- as.data.frame(xyFromCell(r,1:ncell(r)))
xy <- cbind(xy,val)
xy

Usamos el comando geom_raster() para graficar objetos tipo raster.En este caso controlamos los datos NA posiblemente existentes en el raster a partir de la sentencia na.omit(xy), definimos la escala de color a n colores predefinidos por el usuario. Posteriormente incorporamos los puntos que representan las estaciones de monitoreo sobre el mapa. Agregamos etiquetas e incorporamos texto sobre los puntos para indicar los nombres de las estaciones.

file <- 'Meteo_Stations_Quinoas_Coords.txt'  
coords_mt_quinoas = read.table(file, header = T, sep="\t", dec=".")
coords_mt_quinoas

mapa <- ggplot(na.omit(xy)) + geom_raster(aes(x=x, y=y, fill=val)) + 
  coord_equal() +
  scale_fill_gradientn(colours=c("white","orange","red","dark red"), limits=c(0,40)) +
  geom_point(data=coords_mt_quinoas, aes(x=LONG, y=LAT), color="blue") + 
  labs(title="Mapa de Isotermas Cuenca del Quinoas") + 
  geom_text( data=coords_mt_quinoas, hjust=0.5, vjust=-0.5, 
             aes(x=LONG, y=LAT, label=LOCACION), colour="blue2", size=2)

ggsave(file='mapa.pdf',mapa )

Mapa

Mejore este mapa controlando los colores, el rango de valores y tamaño de fuente.

Animaciones En R

Primero instaleremos unos paquetes útiles para la generación de animaciones

install.packages('maptools')
install.packages('animation')
install.packages('RColorBrewer')

Cargamos en memoria las librerías que nos permiten generar una animación GIF mediante la secuencia de varios mapas.

library(ggplot2)
library(maptools)
library(animation)
library(rgdal)
library(raster)
library(RColorBrewer)
library(rasterVis)
library(gridExtra)
library(lubridate)

Cargamos el archivo de coordenadas y nombres de estaciones y a continuación generamos un objeto SpatialPoints para almacenar las coordenadas.

setwd("~/Documents/R_WORKSPACE/Paper_Metereology")
file <- 'data/Meteo_Stations_Quinoas_Coords.txt'  
coords_mt_quinoas = read.table(file, header = T, sep="\t", dec=".")   
head(coords_mt_quinoas)

##                 LOCACION     LONG     LAT
## 1     Matadero en Puente 707391.7 9686976
## 2    Matadero en Sayausi 714620.0 9681624
## 3           Oficinas PNC 709548.7 9685136
## 4 Llaviuco A.J Tomebamba 708354.5 9685555
## 5              Toreadora 697618.7 9692227
## 6           Virgen Cajas 701110.7 9692382

Contabilizamos el tiempo de ejecución de este script.

start_task <- Sys.time ()

Generamos un vector de fechas para poder visualizar en el encabezado del mapa.

fechas <- seq(as.Date('2014-02-01'),as.Date('2015-01-31'),by = 1)
fechas <- format(fechas, format="%d %B %Y")
fechas_str <- as.character(fechas)
step <- seq(1,365,100)

Definimos una paleta de colores aplicable para este ejemplo

myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")), space="Lab")

El comando para generar la animacion es saveGIF().

saveGIF({
  ani.options(nmax = 365) 
  for (i in step) {
    r <- raster('temp_max_dialy_data.tif', band=i)
    val <- getValues(r)
    xy <- as.data.frame(xyFromCell(r,1:ncell(r)))
    xy <- cbind(xy,val)
    
    p <- ggplot(na.omit(xy)) + geom_raster(aes(x=x, y=y, fill=val)) + coord_equal() +
      scale_size(range=c(50,50))+
      scale_fill_gradientn(colours = myPalette(10),limits=c(-2,20),name="Temperatura")  +
      geom_point(data=coords_mt_quinoas, aes(x=LONG, y=LAT), color="blue4",size=3) + 
      geom_text( data=coords_mt_quinoas, hjust=0.5, vjust=-0.5, 
                 aes(x=LONG, y=LAT, label=LOCACION), colour="blue4", size=4 ) + 
      labs(title=paste("MAPA DE ISOTERMAS - MICROCUENCA QUINOAS :: ", fechas_str[i]), x=NULL,y=NULL ) +
      theme(title = element_text(family="sans", size=16,color="blue4",face="bold"),
            legend.position = c(0.85, 0.8),
            legend.background = element_rect(fill = "transparent",colour = NA),
            legend.title = element_text(family="sans",size=10,color="blue4"),
            legend.text = element_text(family="sans",size=10,color="blue4"))
    
    q <- ggplot() + geom_histogram(aes(x=val),alpha=0.75,binwidth = 0.5) + labs(x=NULL,y=NULL) +
      xlim(-2,20) + ylim(0,60000) +
      theme_bw() +
      theme(axis.line = element_line(colour = "black"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            plot.margin = unit(c(0,0,0,0),"mm"),
            panel.border = element_rect(fill = "transparent",colour = NA),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA)) 
    
    ### clear plot area
    grid.newpage()
    ### define first plotting region (viewport)
    vp1 <- viewport(x = 0, y = 0, 
                    height = 1, width = 1,
                    just = c("left", "bottom"),
                    name = "lower left")
    ### enter vp1 
    pushViewport(vp1)
    ### show the plotting region (viewport extent)
    grid.rect(gp = gpar(fill = "white", alpha = 0.2,mar=c(0,0,0,0)))
    ### plot a plot - needs to be printed (and newpage set to FALSE)!!!
    print(p, newpage = FALSE)
    ### leave vp1 - up one level (into root vieport)
    upViewport(1)
    ### define second plot area
    vp2 <- viewport(x = 0.33, y = 0.1, 
                    height = 0.2, width = 0.2,
                    just = c("right", "bottom"),
                    name = "lower right")
    ### enter vp2
    pushViewport(vp2)
    ### show the plotting region (viewport extent)
    grid.rect(gp = gpar(fill =  "transparent", alpha = 0.9,mar=c(0,0,0,0),lwd = 0))
    ### plot another plot
    print(q, newpage = FALSE) 
    ani.pause() 
  }
}, interval = 0.2, movie.name = "bm_demo3day.gif", ani.width = 900, ani.height = 700)

Animación

Recuperamos la hora del sistema para contabilizar el tiempo de ejecución.

end_task <- Sys.time () - start_task
end_task

La animación contiene valores en gris, ¿por qué puede darse esto?. Encontrar la solución a este bug del script.