Hace unos días AEMET cerró el acceso público a sus servidores FTP. Esta decisión contrasta con las estrategias de open data defendidas desde la Comisión Europea, o con las posiciones adoptadas por otros organismos (por ejemplo, CM-SAF o NREL-MIDC).
Afortunadamente, los datos publicados por AEMET hasta el 28 de Octubre de este año han sido recopilados en la página web tiempodiario.com. Con mis limitados conocimientos de HTML y la inestimable ayuda del paquete XML, he aprendido a acceder a la información de radiación solar de manera automática con R. Expongo a continuación lo que he conseguido. Debo avisar de que no siempre funciona bien el acceso a las páginas, y es posible que se produzca un error en el bucle. Una posible solución es reiniciar el bucle a partir del enlace problemático, o añadir controles de error con try:
library(XML)
library(zoo)
mainURL <- 'http://www.tiempodiario.com'
radURL <- 'http://www.tiempodiario.com/radiacion/'
## Días disponibles en tiempodiario.com
days <- seq(as.POSIXct('2010-12-05'), as.POSIXct('2012-10-28'), by='day')
for (i in seq_along(days)){
day <- days[i]
dayURL <- paste(radURL, format(day, '%Y/%m/%d'), '/', sep='')
## Obtengo el contenido HTML de la página del día en cuestión
content <- htmlParse(dayURL, encoding='utf8')
## y extraigo los enlaces a cada estación para ese día
links <- xpathSApply(content, "//a/@href")
linksData <- links[grep('radiacion-estacion', links)]
## Los nombres de las estaciones están escritos en los encabezados H2
h2 <- xpathSApply(content, "//h2")
namesStations <- sapply(h2[-1], function(s){
val <- xmlValue(s)
substr(val, 10, nchar(val))
})
## pero sus códigos de identificación están en los enlaces correspondientes
idStation <- strsplit(linksData, '/')
idStations <- sapply(idStation, function(x)x[4])
dayStationURL <- paste(mainURL, linksData, sep='')
## Accedo a la página de cada estación para ese día
zList <- lapply(dayStationURL, function(ds){
print(ds)
## Obtengo el contenido HTML de la página
content <- htmlParse(ds, encoding='utf8')
## y extraigo las tablas
tabs <- readHTMLTable(content, header=TRUE)
## Cada tabla codifica el índice temporal en la primera
## columna. A pesar de usar header=TRUE, la primera fila
## contiene el nombre las columnas. De ahí que tenga que usar
## -1 en el indexado.
time <- tabs[[1]][-1,1]
dayTime <- format(as.POSIXct(paste(day, time)), '%Y-%m-%d %H:%M')
## El nombre de la variable de cada tabla está codificado en
## los encabezados H2
h2 <- getNodeSet(content, '//h2')
tabLabels <- sapply(h2, xmlValue)
## Sólo me quedaré con las radiación difusa, directa y global
radVars <- c('Radiación Difusa (DF)',
'Radiación Directa (DT)',
'Radiación Global (GL)')
whichTabs <- which(tabLabels %in% radVars)
## Extraigo la información de las tablas seleccionadas. Dado
## que la primera fila contiene los nombres de las columnas,
## todo el contenido es leído como un "factor". De ahí que
## tenga que usar la conversión as.character y después
## as.numeric.
dat <- sapply(tabs[whichTabs],
function(tab)as.numeric(as.character(tab[-1, 2])))
## Finalmente, si todo ha ido bien, construyo un objeto zoo
## para almacenar la serie temporal de esta estación para este
## día. Si algo no ha funcionado obtendré NULL.
if (length(dat) > 0 & !is.list(dat)) {
z <- zoo(dat, as.POSIXct(dayTime))
names(z) <- tabLabels[whichTabs]
} else {
z <- NULL
}
z
})
## El resultado de lapply será una lista de objetos zoo, uno
## para cada estación. Pongo los nombres de las estaciones a
## cada componente de esta lista y grabo el resultado como un
## fichero RData.
names(zList) <- make.names(namesStations)
save(zList, file=paste(day, 'RData', sep='.'))
## Este código se repetirá para cada día, de forma que obtendré
## un conjunto de ficheros RData, uno por día. Cada uno de ellos
## contiene una lista de zoo's, correspondiendo cada componente
## de la lista a una estación de medida.
}
Una vez ejecutado el código anterior obtendremos una colección de ficheros diarios. Si nos interesa tener el registro completo de una estación determinada, es posible obtenerlo recorriendo todos los ficheros:
## Ficheros que he obtenido del bucle anterior.
files <- dir(pattern='RData')
## Por ejemplo, para obtener el registro de la estación de
## Valladolid, recorro todos los ficheros con lapply, y extraigo
## de cada uno el componente que me interesa.
valladolid <- lapply(files, FUN=function(f){
load(f)
dat <- zList$Valladolid
})
## El resultado de lapply será una lista de objetos zoo, que puedo
## transformar en un único objeto zoo con do.call
valladolid <- do.call(rbind, valladolid)
Con otro bucle adicional se puede automatizar el paso anterior y obtener toda la colección de estaciones:
stations <- names(zList)
lapply(stations, FUN=function(name){
datStation <- lapply(files, FUN=function(f){
load(f)
dat <- zList[[name]]
})
datStation <- do.call(rbind, datStation)
write.zoo(datStation, sep=';', file=paste(name, '.txt', sep=''))
})
TiempoDiario.com dijo:
Hola Oscar, soy el webmaster de TiempoDiario.com. Si necesitas los datos de radiación o de cualquier otro tipo con gusto te los podemos dar en el formato que necesites. De todas formas si quieres datos de un día especifico en csv puedes descargartelos desde:
http://www.tiempodiario.com/rayos/aaaammdd/radiacion_datas.csv
Por ejemplo si quieres los datos de radiación del 5 de diciembre de 2012 tendrías que ir a:
http://www.tiempodiario.com/rayos/20101205/radiacion_datas.csv
Oscar Perpiñán Lamigueiro dijo:
¡Muchas gracias por la información!
Claramente esto facilita mucho el trabajo.
Añadiré texto y código para tener estos enlaces en cuenta.
Gracias por vuestra web.
Xavi dijo:
Oscar, eres un crack. Ya me hubiera gustado dominar R y tener copia de estos scripts tuyos hace años cuando recopilaba datos meteorológicos para la tesis.
Un abrazo, y gracias por compartir tus conocimientos, herramientas, truquis varios, etc.
longchamp pliage bags dijo:
Top submit. I anticipate reading additional. Cheers