Web pages of Lars Relund Nielsen

Animation of wind using Shiny

I have played a bit with Shiny the last days. Let us try to create an shiny web application which show a map and add a layer with wind directions and wind speed. First we need weather data for a given latitude and longitude. We will use the API from forecast.io to retrieve data. Your will need an API key for using the service. We define functions to retrieve the data

#' UNIX time (that is, seconds since midnight GMT on 1 Jan 1970)
#'
#' @param x A object to be conveted using \code{as.POSIXct}.
unixTime<-function(x) as.numeric(as.POSIXct(x), origin="1970-01-01")


#' Retrieve weather data for a specific time period
#'
#' Query for a specific time, past or future (for many places, 60 years in the past to 10 years in
#' the future). If \code{from} and \code{to} not are specified then retrive 7 day forecast.
#'
#' @param lat Latitude vector (decimal format)
#' @param lon Longitude vector (decimal format)
#' @param from Start date
#' @param to End date
#'
#' @return A data table of hourly wind data orderd by key unixTime.
#'
#' @note You must have set your \code{FORECASTIO_API_KEY}
getWeatherData<-function(lat, lon, from, to) {
   datH<-NULL
   days<-seq(as.Date(from), as.Date(to), "days")
   for (i in  1:length(days)) {
      #incProgress(detail = paste("  Get data for date:", format(days[i],"%Y-%m-%d"), "\n"))
      cat("  Get data for date:", format(days[i],"%Y-%m-%d"), "\n")
      for (j in 1:length(lat)) {
         cat("    (lat,lon)=(", lat[j],",",lon[j], ")\n", sep="")
         fioList <- get_forecast_for(lat[j],lon[j], timestamp = format(days[i],"%Y-%m-%dT%H:%M:%S"), 
                    units="si", exclude = "currently,minutely,daily,alerts,flags")
         datH<-rbind( datH, 
               cbind(lat=lat[j],lon=lon[j],fioList$hourly[,c('time','windSpeed','windBearing')]) )
      }
   }
   datH$unixTime<-unixTime(datH$time)
   datH$timeStr<-as.character(datH$time)
   datH<-as.data.table(datH)
   setkey(datH,"unixTime")
   #datD$unixTime<-unixTime(datD$time)
   #datD$timeStr<-as.character(datD$time)
   return(datH)
}


#' Calculate distance in kilometers between two points \code{(long1,lat1)} 
#' and \code{(long2,lat2)}
earthDist <- function (long1, lat1, long2, lat2)
{
   rad <- pi/180
   a1 <- lat1 * rad
   a2 <- long1 * rad
   b1 <- lat2 * rad
   b2 <- long2 * rad
   dlon <- b2 - a2
   dlat <- b1 - a1
   a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
   c <- 2 * atan2(sqrt(a), sqrt(1 - a))
   R <- 6378.145
   d <- R * c
   return(d)
}

#' Map object retrived using ggmap package
getMap<-function(cLon,cLat,zoomLevel){
   get_map(c(cLon,cLat), zoom = zoomLevel, maptype = "roadmap", source = "google")
}


#' Retrive weather data for specific date range.
#'
#' @return A data.table with the data 
mergeWeatherData<-function(fromDate=Sys.Date(),toDate=Sys.Date()) {
   # create grid of map points
   dat<-expand.grid(lon = seq(w$ll.lon+dLon/8, w$ur.lon-dLon/8, length.out=cols), 
                    lat = seq(w$ll.lat+dLat/8, w$ur.lat-dLat/8, length.out=cols))
   # retrive wind data
   if (diagDist<30) { # use just the center point to get wind data (faster)
      tmp<-getWeatherData(cLat,cLon, from = fromDate, to = toDate)
      tmp$lat<-NULL; tmp$lon<-NULL
      datH<-do.call("rbind", replicate(length(dat$lon), tmp, simplify = FALSE))
      datH$lon<-rep(dat$lon,each=length(tmp$windSpeed))
      datH$lat<-rep(dat$lat,each=length(tmp$windSpeed))
   } else { # get wind for each coordinate
      datH<-getWeatherData(dat$lat,dat$lon, from = fromDate, to = toDate)
   }
   setkey(datH,"unixTime")
   datH
   # convert wind to arrows
   length=diagDist/3000
   angle <- 90 - datH$windBearing # windBearing = 0 corresponds to north
   datH$lonEnd <- datH$lon+length*cos(angle*pi/180)
   datH$latEnd <- datH$lat+length*sin(angle*pi/180)
   datH$lonText <- datH$lon+dLon/40
   datH$latText <- datH$lat+dLat/40
   datH
}

We can now retrieve the data.

library(Rforecastio)    #install_github("hrbrmstr/Rforecastio")
library(data.table)
library(ggmap)
library(grid)

adr<-"Stauning, Denmark"
cord<-geocode(adr, source = "google")
cLat<-cord$lat
cLon<-cord$lon
zoomLevel = 13
cols<-4 # number of cols/rows on the wind layer
map<-getMap(cLon,cLat,zoomLevel)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=55.96518,8.370522&zoom=13&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
w<-attr(map,"bb")
dLat<-w$ur.lat-w$ll.lat  # delta lat (map width)
dLon<-w$ur.lon-w$ll.lon  # delta lon (map heigth)
diagDist<-earthDist(w$ll.lon,w$ll.lat,w$ur.lon,w$ur.lat) # diag distance in km
Sys.setenv(FORECASTIO_API_KEY = readLines("forecastioApi.txt"))
fromDate = format(Sys.Date(), "%Y-%m-%d", tz="Europe/Paris")
toDate = format(Sys.Date()+1, "%Y-%m-%d", tz="Europe/Paris")
datH<-mergeWeatherData(fromDate,toDate)
##   Get data for date: 2015-11-19 
##     (lat,lon)=(55.96518,8.370522)
##   Get data for date: 2015-11-20 
##     (lat,lon)=(55.96518,8.370522)
datH
##                     time windSpeed windBearing   unixTime             timeStr      lon      lat
##   1: 2015-11-19 00:00:00      9.00         253 1447887600 2015-11-19 00:00:00 8.329409 55.94206
##   2: 2015-11-19 00:00:00      9.00         253 1447887600 2015-11-19 00:00:00 8.356875 55.94206
##   3: 2015-11-19 00:00:00      9.00         253 1447887600 2015-11-19 00:00:00 8.384341 55.94206
##   4: 2015-11-19 00:00:00      9.00         253 1447887600 2015-11-19 00:00:00 8.411807 55.94206
##   5: 2015-11-19 00:00:00      9.00         253 1447887600 2015-11-19 00:00:00 8.329409 55.95743
##  ---                                                                                           
## 764: 2015-11-20 23:00:00      3.86          65 1448056800 2015-11-20 23:00:00 8.411807 55.97281
## 765: 2015-11-20 23:00:00      3.86          65 1448056800 2015-11-20 23:00:00 8.329409 55.98818
## 766: 2015-11-20 23:00:00      3.86          65 1448056800 2015-11-20 23:00:00 8.356875 55.98818
## 767: 2015-11-20 23:00:00      3.86          65 1448056800 2015-11-20 23:00:00 8.384341 55.98818
## 768: 2015-11-20 23:00:00      3.86          65 1448056800 2015-11-20 23:00:00 8.411807 55.98818
##        lonEnd   latEnd  lonText  latText
##   1: 8.326323 55.94112 8.332156 55.94360
##   2: 8.353789 55.94112 8.359622 55.94360
##   3: 8.381255 55.94112 8.387087 55.94360
##   4: 8.408721 55.94112 8.414553 55.94360
##   5: 8.326323 55.95649 8.332156 55.95897
##  ---                                    
## 764: 8.414731 55.97417 8.414553 55.97434
## 765: 8.332334 55.98954 8.332156 55.98972
## 766: 8.359799 55.98954 8.359622 55.98972
## 767: 8.387265 55.98954 8.387087 55.98972
## 768: 8.414731 55.98954 8.414553 55.98972

Next step is to add the layer

#' Draw wind layer at a specific unix time
drawWind<-function(uTime,dat) {
   p <- ggmap(map, extent = "device") +
      geom_segment(aes(x = lon, xend = lonEnd, y = lat, yend = latEnd), data=dat, size=1, 
                   color="black", arrow = arrow(length = unit(0.3,"cm"), ends = "first"), 
                   alpha=0.5, na.rm = TRUE) +
      geom_text(data=dat, mapping=aes(x=lonText, y=latText, label=round(windSpeed,0)), 
                size=8, hjust=0.45, vjust=0.4) +
      annotate('text', x=w$ur.lon-dLon/20, y=w$ur.lat-dLat/20, 
               label = format(as.POSIXct(uTime, origin="1970-01-01"), format="%d-%m-%Y %H:%M"), 
               colour = I('black'), size = 12, hjust=1) +
      theme(legend.position = "none")
   print(p)
}

drawWind(unixTime(Sys.Date()),datH[J(unixTime(Sys.Date()))])

 

We are now ready to make a shiny app that load the wind map for different hours using a slider. The user interface file becomes:

library(shiny)
shinyUI(fluidPage(
   tags$link(rel = 'stylesheet', type = 'text/css', href = 'progress.css'),

   titlePanel(h4("Wind direction and speed (today and tomorrow)", align="center")),

   fluidRow(
      column(12,
         uiOutput("ui")
      ), align="center"
   ),
   br(),
   fluidRow(
      column(12,
             sliderInput('myslider',
                         'Hour (next day +24)',
                         min=0, max=47, value=0, post = ":00",
                         animate=animationOptions(interval = 2000,loop=F)
             )
      ), align="center"
   )
))

and the server script:

library(shiny)
library(Rforecastio)    #install_github("hrbrmstr/Rforecastio")
library(data.table)
library(ggmap)
library(grid)
library(animation)
source("helpers.R", local=TRUE)

# center of map
adr<-"Stauning, Denmark"
cord<-geocode(adr, source = "google")
cLat<-cord$lat
cLon<-cord$lon
zoomLevel = 13
map<-getMap(cLon,cLat,zoomLevel)
w<-attr(map,"bb")
dLat<-w$ur.lat-w$ll.lat  # delta lat (map length lat)
dLon<-w$ur.lon-w$ll.lon  # delta lon (map length lon)
cols<-4
diagDist<-earthDist(w$ll.lon,w$ll.lat,w$ur.lon,w$ur.lat)
oopt <- animation::ani.options(interval = 0.5, nmax=48)


shinyServer(function(input, output, session) {

   observe({
      curHour<-as.integer(format(Sys.time(), "%H", tz="Europe/Paris") )
      fromDate = format(Sys.Date(), "%Y-%m-%d", tz="Europe/Paris")
      toDate = format(Sys.Date()+1, "%Y-%m-%d", tz="Europe/Paris")
      updateSliderInput(session, "myslider", value = curHour)  # use plot for current hour

      # should we remake the plots (done if current plots are 4 hours old)
      if (file.exists("lastRun.csv")) {
         lastGen<-read.csv("lastRun.csv")
         if (lastGen$fromDate!=fromDate || (lastGen$fromDate==fromDate & 
             curHour-lastGen$curHour>4) || length(list.files("www/img"))<23)
            rerun = TRUE else rerun <- FALSE
      } else {rerun <- TRUE }

      if (rerun) {
         withProgress(message = 'Making plot (be patient) ...', value = 0, {
            Sys.setenv(FORECASTIO_API_KEY = readLines("forecastioApi.txt"))
            datH<-mergeWeatherData(fromDate,toDate)
                        dir.create("www")
            dir.create("www/img")
            unlink("www/img/*", recursive = TRUE)
            incProgress(detail = paste0("Save pictures") )
            v = unique(datH$unixTime)
            j<-1
            lapply(v, function(i) {
               dat<-datH[J(i)]
               png(paste0("www/img/wind",j,".png"), width = 1280, height = 1280)
               drawWind(i,dat)
               dev.off()
               j<<-j+1
            })
            write.csv(data.frame(fromDate=fromDate,curHour=curHour), file="lastRun.csv")
            rerun<<-FALSE
         })
      }
   })

   imgurl <- reactive({
      i=input$myslider
      return(paste0("./img/wind",i+1,".png")) #the path of pictures
   })

   output$ui <- renderUI({
      tags$div(
         tags$img(src = imgurl(),  width="100%", height="100%", style="max-width:600px")
      )
   })

})

The final result can be see at shinyapps.io and the code is available at GitHub.

Leave a Reply