Y además vamos a analizar si de verdad llueve más los fines de semana en Madrid. Hace tiempo que me gustaría estudiar la influencia de la contaminación en algunos fenómenos atmosféricos. Por supuesto no tengo tiempo para elaborar un estudio de ese tipo. La base de este estudio iba a ser el paquete solaR. Por otro lado quería elaborar un monográfico sobre el paquete chron que contiene funciones muy interesantes para el manejo de fechas. Sin tiempo es imposible, por ello nos vamos a acercar a estos dos paquetes con un ejemplo y de propina os presento directlabels otro paquete muy interesante para añadir etiquetas a nuestros gráficos.
Nuestro trabajo va a comenzar con la descarga de datos agroclimáticos con solaR. Este paquete nos permite conectarnos con el SIAR (descripción) Sistema de Información Agroclimática para el Regadío:
library(solaR)
siar = readMAPA(prov=19, est=9, start='01/01/2000', end='22/05/2011')
#str(siar)
datos = data.frame(siar@data$Precipitacion)
names(datos)="lluvias"
La función readMAPA como Dora la Exploradora se conecta a MAPA y nos permite descargarnos una gran cantidad de datos. Por cierto, para ver MAPA usad IE. Esta función requiere la provincia (prov), la estación (est), fecha inicial y fecha final de la extracción. La estación más cercana a mi hogar está en Marchamalo Guadalajara. Para saber las estaciones usad los datos de la consulta que lanzáis en la web. Para el ejemplo sólo nos quedamos con las precipitaciones y creamos un data.frame con lluvias y fechas como nombres de filas. Ahora hemos de trabajar con las fechas y para ello empleamos chron:
library(chron)
datosdia = factor(weekdays(as.Date(row.names(datos),"%Y-%m-%d")))
datosmes = factor(as.yearmon(row.names(datos)))
datos$finde = factor(ifelse(is.weekend(as.Date(row.names(datos),"%Y-%m-%d")),"FINDE","RESTO"))
summary(datos)
Los nombres de fila de nuestro data.frame los transformamos en 3 factores. Con weekday (es una función del paquete Base) obtenemos el día de la semana a partir de una fecha. La función as.yearmon nos crea una variable MES-AAAA y la función is.weekend toma valores 1 y 0 si se trata de fin de semana. Ahora tenemos que estudiar si se produce una diferencia entre las precipitaciones los fines de semana y entre semana. Esto se hace con un análisis de la varianza, pero nosotros vamos a emplear lattice para crear un gráfico de densidades por meses con R:
library(lattice)
densityplot( ~ lluvias | mes , data = datos, groups=finde ,
layout = c(3, 4), xlab = "Estevita", bw = 2, plot.points = FALSE)
El resultado:
Echamos de menos un etiquetado para el gráfico de densidades que hemos creado con densityplot, para ello vamos a emplear directlabels:
#install.packages("directlabels")
library(directlabels)
densidades = densityplot( ~ lluvias | mes , data = datos, groups=finde ,
layout = c(3, 4), xlab = "Estevita", bw = 2, plot.points = FALSE)
direct.label(densidades)
El resultado es:
Este gráfico es claramente mejorable, además incluye un error pero mi hijo acaba de entrar por la puerta lleno de tierra con una espátula y… En fin, los resultados obtenidos son interesantes. Como ejercicio os dejo lo siguiente: réplica del gráfico con ggplot2 y análisis de la varianza.
Hola,
Gracias por el ejemplo. solaR tiene varios métodos para visualizar los datos que contienen sus objetos pero no se me había ocurrido incluir uno para dibujar densidades. ¡Buena idea! Lo incluiré en la siguiente versión.
También me ha servido para comprobar que debo actualizar el contenido de RedEstaciones (data(RedEstaciones)), porque precisamente esta estación #9 no está incluida en la lista.
Aprovecho para proponerte una forma diferente de hacer la tarea que relatas en tu entrada.
Saludos.
Oscar.
##############################################
library(solaR)
siar = readMAPA(prov=19, est=9, start=’01/01/2000′, end=’22/05/2011′)
##solaR tiene un método para extraer la información
##sin necesidad de referencias a sus «slots»
##El resultado de getData es un «zoo»
datos = getData(siar)
##Un objeto zoo tiene un índice temporal al que
##se puede acceder con index para obtener
##un POSIXct
idx <- index(datos)
##Para determinar si estamos o no en fin de semana
##el FAQ de zoo propone esta simple función
is.weekend <- function(x) ((as.numeric(x) – 2)%%7) < 2
datosfinde <- is.weekend(idx)
##A partir de aquí viene lo que haría el método
##densityplot para una «signature» definida por
##»formula, zoo», al estilo de las que
##ya están definidas en solaR para xyplot y levelplot
##Primero pasamos de «zoo» a «data.frame» con
##las variables que nos interesan
lluvia <- as.data.frame(datos)[c(‘Precipitacion’, ‘finde’)]
##luego añadimos variables basadas en el índice
##temporal utilizando algunas funciones sencillas
##incluidas en solaR
lluviadiaAnno=doy(idx)
lluviames=month(idx)
lluviaanno=year(idx)
##Ya podemos pintar
p <- densityplot( ~ Precipitacion | mes ,
data = lluvia,
groups=finde,
layout = c(3, 4),
xlab = «Estevita»,
bw = 2, plot.points = FALSE,
strip=strip.custom(
strip.levels=c(TRUE, TRUE))
)
##Y como alternativa al paquete direct.labels
##a mí me gustan las funcionalidades de glayer
##del paquete latticeExtra
##Primero definimos una función que etiqueta
label.densityplot <- function(x, finde, col.line){
d <- density(x, bw=2)
i <- which.max(dy)
lab <- ifelse(finde, ‘Fin’, ‘Lab’)
ltext(dx[i],d$y[i],lab,adj=c(0.2,0), pos=4, col=col.line, cex=0.7)
}
##y se la pasamos a glayer para que dibuje
p+glayer(label.densityplot(x, group.value, col.line))
Muchas gracias por el aporte. Voy a explorar las posibilidades gráficas de solaR y habrá nueva entrada al respecto.
Me han gustado más las posibilidades de latticeExtra que direct.labels a ver si puedo hacer un repaso rápido.