Analysis of a Philadelphia Weather Station Data during Hurricane Sandy

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.

Last Updated:

Sys.time()
## [1] "2012-10-30 15:31:38 EDT"

Load packages and load data

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")

Do date time conversions

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"))

Barometric pressure

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()

plot of chunk unnamed-chunk-5

As of now, the barometric pressure hasn't bottomed out.
The barometric pressure bottomed out around 9pm.

Windspeed

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")

plot of chunk unnamed-chunk-6

Rain

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)")

plot of chunk unnamed-chunk-7

Wind Direction

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")

plot of chunk unnamed-chunk-8


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")

plot of chunk unnamed-chunk-9

Comparison to Atlantic City

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...)

Date time conversions

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

Barometric Pressure

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()

plot of chunk unnamed-chunk-12

Windspeed

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")