I've pulled this data from the Weather Underground Station at Broad and Washington in Philadelphia, which is the closest one to my house. The weather station closest to my house doesn't appear to have been updated since 7pm. I'm now displaying the data from the weather station in University City.
Sys.time()
## [1] "2012-10-30 15:31:38 EDT"
library(ggplot2)
library(scales)
library(plyr)
library(reshape2)
library(circular)
### Load and save the data, since I believe the weather station only reports data from the current calendar day.
get_more_info <- function(data= NULL, url){
require(XML)
library(plyr)
if(is.null(data)){
out <- xmlToDataFrame(url)
}else{
new_data <- xmlToDataFrame(url)
new_data <- subset(new_data, !(observation_time %in% data$observation_time))
if(nrow(new_data)>0){
out <- rbind.fill(data, new_data)
cat("adding ")
cat(nrow(new_data))
cat(" rows of data\n")
}else{
out <- data
}
}
return(out)
}
if(file.exists("sandy.RData")){
load("sandy.RData")
sandy.df <- get_more_info(data = sandy.df, url = "http://api.wunderground.com/weatherstation/WXDailyHistory.asp?ID=KPAPHILA20&format=XML" )
}else{
sandy.df <- get_more_info(url = "http://api.wunderground.com/weatherstation/WXDailyHistory.asp?ID=KPAPHILA20&format=XML" )
}
save(sandy.df,file="sandy.RData")
sandy.df$Time <- as.POSIXct(sandy.df$observation_time_rfc822,
format = "%a, %d %b %Y %X", tz = "GMT")
sandy.df$Time <- as.POSIXct(format(sandy.df$Time, tz = "America/New_York"))
sandy.df$pressure_mb.n <- as.numeric(as.character(sandy.df$pressure_mb))
ggplot(sandy.df, aes(Time, pressure_mb.n)) +
geom_step() +
stat_smooth(span = 0.3)+
scale_x_datetime(labels = date_format("%I:%M %p"))+
ggtitle("barometric pressure")+
theme_bw()
As of now, the barometric pressure hasn't bottomed out.
The barometric pressure bottomed out around 9pm.
sandy.df$wind_mph.n <- as.numeric(as.character(sandy.df$wind_mph))
sandy.df$wind_gust_mph.n <- as.numeric(as.character(sandy.df$wind_gust_mph))
winds <- subset(melt(sandy.df[,c("Time","wind_mph.n","wind_gust_mph.n")], id = "Time"), value > -100)
ggplot(winds, aes(Time, value, color = variable))+
geom_point()+
ylim(0,60)+
scale_x_datetime(labels = date_format("%I:%M %p"))+
ylab("miles per hour")+
stat_smooth(span = 0.3)+
scale_color_brewer(palette = "Set1", labels = c("wind","gust"))+
theme_bw()+
ggtitle("Sustained wind speeds, and gusts")
sandy.df$precip_today_in.n <- as.numeric(as.character(sandy.df$precip_today_in))
midnight <- as.POSIXct("2012-10-29 23:59:59 EDT")
max_mon <- max(subset(sandy.df, Time <= midnight)$precip_today_in.n)
sandy.df$precip_today_in.n[sandy.df$Time >= midnight] <- sandy.df$precip_today_in.n[sandy.df$Time >= midnight] + max_mon
sandy.df$precip_1hr_in.n <- as.numeric(as.character(sandy.df$precip_1hr_in))
rain <- melt(sandy.df[,c("Time","precip_today_in.n", "precip_1hr_in.n")], id = "Time")
levels(rain$variable) <- c("cumulative","per-hour")
rain <- subset(rain, value > -1)
ggplot(rain, aes(Time, value))+
geom_step()+
ylab("inches")+
scale_x_datetime(labels = date_format("%I:%M %p"))+
facet_grid(variable~., scales = "free_y")+
theme_bw()+
ggtitle("Inches of rain (cumulative and per-hour)")
sandy.df$wind_degrees.n <- as.numeric(as.character(sandy.df$wind_degrees))
sandy.df <- subset(sandy.df, wind_degrees.n >= 0 & wind_degrees.n <= 360)
wind_dens <- density(as.circular(sandy.df$wind_degrees.n,
control.circular = list(units = "degrees")),
bw = 30)
wind.dens.df <- as.data.frame(wind_dens[c("x","y")])
ggplot(wind.dens.df, aes(x, y))+
geom_area()+
scale_x_continuous(breaks = c(0,90,180,180 + 90), labels = c("North", "East","South","West"))+
coord_polar()+
theme_bw()+
ggtitle("Wind direction")
sandy.df$Cut_Time <- cut(sandy.df$Time, round((((diff(as.numeric(range(round.POSIXt(sandy.df$Time, "hours")))))/60)/60)/4))
circ_dens <- dlply(sandy.df,
.(Cut_Time),
with,
density(as.circular(wind_degrees.n, control.circular =
list(units = "degrees")),
bw = 30))
circ_dens.df <- ldply(circ_dens, function(x)as.data.frame(x[c("x","y")]))
ggplot(circ_dens.df, aes(x, y))+
geom_area()+
scale_x_continuous(breaks = c(0,90,180, 180+90),labels = c("North","East","South","West"))+
facet_wrap(~Cut_Time) +
coord_polar()+
theme_bw()+
ggtitle("Wind direction in approximately 4 hour increments")
I believe this data is from the Atlantic City Airport.
sandy.nj <- ldply(c("http://www.wunderground.com/history/airport/KACY/2012/10/29/DailyHistory.html?format=1","http://www.wunderground.com/history/airport/KACY/2012/10/30/DailyHistory.html?format=1"), read.csv)
sandy.nj$DateUTC.br... <- gsub("<br />", "", sandy.nj$DateUTC.br...)
sandy.nj$Time <- as.POSIXct(sandy.nj$DateUTC.br..., tz = "GMT")
sandy.nj$Time <- as.POSIXct(format(sandy.nj$Time, tz = "America/New_York"))
## Also convert inHg to mb
sandy.nj$pressure_mb.n <- sandy.nj$Sea.Level.PressureIn * 33.8637526
ggplot(subset(sandy.nj, pressure_mb.n < 1025), aes(Time, pressure_mb.n)) +
geom_step() +
stat_smooth(span = 0.3)+
scale_x_datetime(labels = date_format("%I:%M %p"))+
ggtitle("barometric pressure")+
theme_bw()
sandy.nj$wind_mph.n <- as.numeric(as.character(sandy.nj$Wind.SpeedMPH))
sandy.nj$wind_gust_mph.n <- as.numeric(as.character(sandy.nj$Gust.SpeedMPH))
winds.nj <- subset(melt(sandy.nj[,c("Time","wind_mph.n","wind_gust_mph.n")], id = "Time"), value > -100)
ggplot(winds.nj, aes(Time, value, color = variable))+
geom_point()+
expand_limits(y = 0)+
scale_x_datetime(labels = date_format("%I:%M %p"))+
ylab("miles per hour")+
stat_smooth(span = 0.3)+
scale_color_brewer(palette = "Set1", labels = c("wind","gust"))+
theme_bw()+
ggtitle("Sustained wind speeds, and gusts")
Note: This rain data is not completely comparable to the data from the Philadelphia weather station.
sandy.nj$PrecipitationIn <- as.numeric(as.character(sandy.nj$PrecipitationIn))
sandy.nj$PrecipitationIn[is.na(sandy.nj$PrecipitationIn)] <- 0
sandy.nj$TotalPrecip <- cumsum(sandy.nj$PrecipitationIn)
rain.nj <- melt(sandy.nj[,c("Time","PrecipitationIn", "TotalPrecip")], id = "Time")
levels(rain.nj$variable) <- c("per-hour", "cumulative")
rain.nj$variable <- relevel(rain.nj$variable, "cumulative")
ggplot(rain.nj, aes(Time, value))+
geom_step()+
ylab("inches")+
scale_x_datetime(labels = date_format("%I:%M %p"))+
facet_grid(variable~., scales = "free_y")+
theme_bw()+
ggtitle("Inches of rain (cumulative and per-hour)")
wind_dens.nj <- density(as.circular(sandy.nj$WindDirDegrees,
control.circular = list(units = "degrees")),
bw = 30)
wind_dens.nj.df <- as.data.frame(wind_dens.nj[c("x","y")])
ggplot(wind_dens.nj.df, aes(x, y))+
geom_area()+
scale_x_continuous(breaks = c(0,90,180,180 + 90), labels = c("North", "East","South","West"))+
coord_polar()+
theme_bw()+
ggtitle("Wind direction")
sandy.nj$Cut_Time <- cut(sandy.nj$Time, round((((diff(as.numeric(range(round.POSIXt(sandy.nj$Time, "hours")))))/60)/60)/4))
circ_dens.nj <- dlply(sandy.nj,
.(Cut_Time),
with,
density(as.circular(WindDirDegrees, control.circular =
list(units = "degrees")),
bw = 30))
circ_dens.nj.df <- ldply(circ_dens.nj, function(x)as.data.frame(x[c("x","y")]))
ggplot(circ_dens.nj.df, aes(x, y))+
geom_area()+
scale_x_continuous(breaks = c(0,90,180, 180+90),labels = c("North","East","South","West"))+
facet_wrap(~Cut_Time) +
coord_polar()+
theme_bw()+
ggtitle("Wind direction in approximately 4 hour increments")
sandy.df$Location <- "Philly"
sandy.nj$Location <- "AC"
bar_comp <- rbind.fill(sandy.df[,c("Location","Time","pressure_mb.n")],
sandy.nj[,c("Location","Time","pressure_mb.n")])
ggplot(subset(bar_comp, pressure_mb.n < 1025), aes(Time, pressure_mb.n, color = Location))+
geom_step()+
stat_smooth(span = 0.3)+
scale_x_datetime(labels = date_format("%I:%M %p"))+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
ggtitle("barometric pressure")
winds$Location <- "Philly"
winds.nj$Location <- "AC"
winds_comp <- rbind.fill(winds, winds.nj)
levels(winds_comp$variable) <- c("wind","gust")
ggplot(winds_comp, aes(Time, value, color = Location))+
geom_point()+
expand_limits(y = 0)+
stat_smooth(span = 0.3)+
scale_x_datetime(labels = date_format("%I:%M"))+
ylab("miles per hour")+
scale_color_brewer(palette = "Dark2")+
facet_wrap(~variable)+
theme_bw()+
ggtitle("Sustained wind speeds, and gusts")
The cumulative rain figures either for Philadelphia or Atlantic City must be wrong, because the inches of rain per-hour don't look to be substantially different between the two.
rain$Location <- "Philly"
rain.nj$Location <- "AC"
rain_comp <- rbind.fill(rain, rain.nj)
ggplot(rain_comp, aes(Time, value, color = Location))+
geom_step()+
ylab("inches")+
scale_x_datetime(labels = date_format("%I:%M %p"))+
scale_color_brewer(palette= "Dark2")+
facet_grid(variable~., scales = "free_y")+
theme_bw()+
ggtitle("Inches of rain (cumulative and per-hour)")
wind.dens.df$Location <- "Philly"
wind_dens.nj.df$Location <- "AC"
wind_dens_comp <- rbind.fill(wind.dens.df, wind_dens.nj.df)
ggplot(wind_dens_comp, aes(x, y, fill = Location))+
geom_area(alpha = 0.8)+
scale_x_continuous(breaks = c(0,90,180,180 + 90), labels = c("North", "East","South","West"))+
scale_fill_brewer(palette = "Dark2")+
coord_polar()+
theme_bw()+
facet_grid(.~Location)+
ggtitle("Wind direction")
colnames(sandy.nj)[colnames(sandy.nj)=="WindDirDegrees"] <- "wind_degrees.n"
sandy.nj$Location <- "AC"
sandy.df$Location <- "Philly"
comp_to_cut <-rbind.fill(sandy.nj[,c("Time","Location","wind_degrees.n")],
sandy.df[,c("Time","Location","wind_degrees.n")])
comp_to_cut$Cut_Time <- cut(comp_to_cut$Time, round((((diff(as.numeric(range(round.POSIXt(comp_to_cut$Time, "hours")))))/60)/60)/6))
circ_dens_comp<- dlply(comp_to_cut,
.(Cut_Time,Location),
with,
density(as.circular(wind_degrees.n, control.circular =
list(units = "degrees")),
bw = 30))
circ_dens_comp.df <- ldply(circ_dens_comp, function(x)as.data.frame(x[c("x","y")]))
ggplot(circ_dens_comp.df, aes(x, y, fill = Location))+
geom_area(alpha = 0.8)+
scale_x_continuous(breaks = c(0,90,180, 180+90),labels = c("North","East","South","West"))+
scale_fill_brewer(palette = "Dark2")+
facet_grid(Cut_Time~Location) +
coord_polar()+
theme_bw()+
ggtitle("Wind direction in approximately 6 hour increments")
