Etiquetas

, , , , ,

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=''))
})
About these ads