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()
16 Comments
Nice. Can I add this to the graph gallery ?
Absolutely. The best one to use is prbly the last one I generated, available via GitHub https://github.com/hrbrmstr/sandy
Nice. Will do shortly and let you know.
Done. http://gallery.r-enthusiasts.com/graph/SuperStorm_Sandy,170
Done, see http://gallery.r-enthusiasts.com/graph/SuperStorm_Sandy,170
Awesome. Again, much obliged. Keep spreading the goodness of ‘R’!
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.
Nice catch. Thx. I’m mostly in Linux & Mac OS land. I’ll see if I can put an OS-test wrapper around it.
This does highlight for me what a pain it is that rgdal won’t install under cygwin.
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.
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.
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.
I just posted (literally) some alternate code up on the github repository : https://github.com/hrbrmstr/sandy … give it a try (it won’t be perfect, but you’ll get some results)
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
>
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/
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.
2 Trackbacks/Pingbacks
[…] because with something like this it is just a matter of who has an interesting approach. Bob Rudis showed how he was tracking Sandy using R, and when he plotted his map it still looked like Sandy was […]
[…] have some history with hurricane data and thought I’d jump on the bandwagon using the same data and making some […]