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.
