Watch Sandy in “R” (Including Forecast Cone)

As indicated in the code comments, Google took down the cone KML files. I’ll be changing the code to use the NHC archived cone files later tonight

NOTE: There is significantly updated code on github for the Sandy ‘R’ dataviz.

This is a follow-up post to the quickly crafted Watch Sandy in “R” post last night. I noticed that Google provided the KML on their crisis map and wanted to show how easy it is to add it to previous code.

I added comments inline to make it easier to follow along.

library(maps)
library(maptools)
library(rgdal)
 
# need this for handling "paste" extra spaces
trim.trailing <- function (x) sub("\\s+$", "", x)
 
# get track data from Unisys
wx = read.table(file="http://weather.unisys.com/hurricane/atlantic/2012/SANDY/track.dat", skip=3,fill=TRUE)
 
# join last two columns (type of storm) that read.table didn't parse as one
# and give the data frame real column names
wx$V7 = trim.trailing(paste(wx$V7,wx$V8," "))
wx$V8 = NULL
colnames(wx) = c("Advisory","Latitude","Longitude","Time","WindSpeed","Pressure","Status")
 
# annotate the name with forecast (+12/+24/etc) projections
wx$Advisory = unlist(strsplit(toString(wx$Advisory),", "))
wx$a = ""
wx$a[grep("\\+",wx$Advisory)] = wx$Advisory[grep("\\+",wx$Advisory)]
wx$Status = trim.trailing(paste(wx$Status,wx$a," "))
 
# cheap way to make past plots one color and forecast plots another
wx$color = "red"
wx$color[grep("\\+",wx$Advisory)] = "orange"
 
# only want part of the map
xlim=c(-90,-60)
 
# plot the map itself
map("state", interior = FALSE, xlim=xlim)
map("state", boundary = FALSE, col="gray", add = TRUE,xlim=xlim)
 
# plot the track with triangles and colors
lines(x=wx$Longitude,y=wx$Latitude,col="black",cex=0.75)
points(x=wx$Longitude,y=wx$Latitude,col=wx$color,pch=17,cex=0.75)
 
# annotate it with the current & projects strength status + forecast
text(x=wx$Longitude,y=wx$Latitude,col='blue',labels=wx$Status,adj=c(-.15),cex=0.33)
 
# get the forecast "cone" from the KML that google's crisis map provides
# NOTE: you need curl & unzip (with funzip) on your system
# NOTE: Google didn't uniquely name the cone they provide, so this will break
#       post-Sandy. I suggest saving the last KML & track files out out if you
#       want to save this for posterity
 
# this gets the cone and makes it into a string that getKMLcoordinates can process
coneKML = paste(unlist(system("curl -s -o - http://mw1.google.com/mw-weather/maps_hurricanes/three_day_cone.kmz | funzip", intern = TRUE)),collapse="")
 
# process the coneKML and make a data frame out of it
coords = getKMLcoordinates(textConnection(coneKML))
coords = data.frame(SpatialPoints(coords, CRS('+proj=longlat')))
 
# this plots the cone and leaves the track visible
polygon(coords$X1,coords$X2, pch=16, border="red", col='#FF000033',cex=0.25)
 
# no one puts Sandy in a box (well, except us)
box()

Cover image from Data-Driven Security
Amazon Author Page

18 Comments Watch Sandy in “R” (Including Forecast Cone)

  1. John Johnson

    Very nice.

    One thing, if you try to run this on Windows, you’ll need to replace system with shell. Otherwise the shell pipe won’t work.

    Reply
  2. Bob

    This does highlight for me what a pain it is that rgdal won’t install under cygwin.

    Reply
    1. hrbrmstr

      Do you have enough memory for a linux VM under VirtualBox? Also, you could also spin up a free EC2 instance and run
      R/RStudio there.

      Reply
  3. Ralph

    Cone is not working on my Windows Machine.. It’s throwing errors like :Error in !unlist(lapply(obj, is.numeric)) : invalid argument type! Also, no altitude values for KML object 1

    Any suggestion? Thanks.

    Reply
    1. hrbrmstr

      As predicted, Google removed the cone KML files right after Sandy was done being an ocean event. I’ll be changing the code to work with the NHC cone files : http://www.nhc.noaa.gov/gis/archiveforecastresults.php?id=al18&year=2012&name=Hurricane%20SANDY : but you can – in the interim – try grabbing the latest KMZ file from that page and running it through some R code to get the necessary polygon.

      Reply
  4. Bob

    Hi, I’m a beginner with R, I’m trying to test your code on windows with R2.15.2

    It seems I have a problem with funzip command. (curl is ok, funzip is not recognized)
    ..some tips for me?

    > coneKML5 = paste(unlist(system(“curl -s -o – http://www.nhc.noaa.gov/storm_graphics/api/AL182012_031adv_CONE.kmz | funzip”, intern = TRUE)),collapse=””)

    Warning message:
    running command ‘curl -s -o – http://www.nhc.noaa.gov/storm_graphics/api/AL182012_031adv_CONE.kmz | funzip’ had status 6

    >
    > # process the coneKML and make a data frame (lat/long/elevation coords) out of it
    > coords = getKMLcoordinates(textConnection(coneKML5),ignoreAltitude=TRUE)
    Errore in ifelse(ignoreAltitude, list(m[, 1:2]), list(m)) :
    indice fuori limite
    > coords = data.frame(SpatialPoints(coords[1], CRS(‘+proj=longlat’)))
    Errore in coordinates(coords) :
    error in evaluating the argument ‘obj’ in selecting a method for function ‘coordinates’: Errore: oggetto “coords” non trovato
    >

    Reply
    1. hrbrmstr

      It sounds like either the ‘unzip’ package is not installed or it doesn’t contain the the ‘funzip’ command. ‘funzip’ was developed primarily due to the fact that ‘unzip’ cannot read STDIN from a pipe. If you can’t add ‘unzip’ or your linux variant doesn’t have ‘funzip’ as part of the package, try building/installing ‘funzip’ from: http://www.info-zip.org/pub/infozip/

      Reply
  5. charlie

    hi hrbrmstr. thanks for this. trying to extract fema damage assessment data from either google crisis map (http://google.org/crisismap/2012-sandy-nyc) or from the fema site sourced by google (http://fema-services2.esri.com/arcgis/rest/services/2012_Sandy/ImageCat_NLT/MapServer)

    have not worked with these kinds of sources before. have been searching for a few hours for how to do it, but no luck. any advice? really want to get at the spatially referenced tabular data usually found in either a dbf or klm file

    any help, or even direction to possible help, greatly appreciated.

    Reply
  6. Pingback: Watching Hurricane Sandy Graze by Toronto in R -

  7. Pingback: Plot Me Like a Hurricane (a.k.a. animating historical North Atlantic basin tropical storm tracks) | rud.is

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.