Maximum Windspeeds from Hurricane Katrina using H*Wind Analyses

Rob Hodges, PhD

http://robhodges.co

This is my second post using AOML/NOAA's H*Wind data for analysis using R.

My goal here is to show the cumulative windfield maxima from Hurricane Katrina (2005).

First, we library in the appropriate R packages.

library(rgdal)
## Loading required package: sp
## rgdal: version: 0.7.25, (SVN revision 400) Geospatial Data Abstraction
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
## 1.9.0, released 2011/12/29 Path to GDAL shared files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to
## PROJ.4 shared files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/proj
library(colorspace)
library(maps)
library(mapdata)
library(maptools)
## Loading required package: foreign
## Loading required package: grid
## Loading required package: lattice
## Checking rgeos availability: TRUE

We also set the working directory.

setwd("/Users/robhodges17/GDrive/HWind/Katrina2005")

Here are the urls for each analyzed windfield from Hurricane Katrina, available every 3 hours starting with its position just off the Florida coast. Note: If you wanted to examine other storms, you can copy-and-paste in the corresponding URLs.

urls = c("ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0825/2100/AL122005_0825_2100shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/0000/AL122005_0826_0000shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/0300/AL122005_0826_0300shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/0600/AL122005_0826_0600shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/0900/AL122005_0826_0900shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/1200/AL122005_0826_1200shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/1500/AL122005_0826_1500shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/1800/AL122005_0826_1800shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0826/2100/AL122005_0826_2100shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/0000/AL122005_0827_0000shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/0300/AL122005_0827_0300shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/0600/AL122005_0827_0600shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/0900/AL122005_0827_0900shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/1200/AL122005_0827_1200shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/1500/AL122005_0827_1500shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/1800/AL122005_0827_1800shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0827/2100/AL122005_0827_2100shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/0000/AL122005_0828_0000shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/0300/AL122005_0828_0300shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/0600/AL122005_0828_0600shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/0900/AL122005_0828_0900shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/1200/AL122005_0828_1200shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/1500/AL122005_0828_1500shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/1800/AL122005_0828_1800shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0828/2100/AL122005_0828_2100shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/0000/AL122005_0829_0000shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/0300/AL122005_0829_0300shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/0600/AL122005_0829_0600shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/0900/AL122005_0829_0900shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/1200/AL122005_0829_1200shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/1500/AL122005_0829_1500shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/1800/AL122005_0829_1800shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0829/2100/AL122005_0829_2100shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0830/0000/AL122005_0830_0000shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0830/0300/AL122005_0830_0300shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0830/0600/AL122005_0830_0600shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0830/0900/AL122005_0830_0900shape.tar.gz", 
    "ftp://ftp.aoml.noaa.gov/hrd/pub/hwind/PostAnalysis/2005/AL122005/0830/1200/AL122005_0830_1200shape.tar.gz")

You'll notice that each of the above URLs is a link to a zipped (“tarred”) file. I'll use a for-loop to download, unzip, and read in the corresponding spatial data from a temporary directory. Your download speed will influence the duration of this code chunk.

for (i in 1:length(urls)) {
    file = basename(urls)[i]
    download.file(urls[i], paste(getwd(), "/", file, sep = ""), quiet = T)
    untar(file)
}

To produce a windswath of Hurricane Katrina's windfields, we first combine all of the spatial data. We read in the first downloaded vector data file, specify which columns we want to keep, and update it with the remaining shapefiles.

files = list.files(pattern = "*.shp$", recursive = TRUE, full.names = TRUE)
spdf.data = readOGR(dsn = getwd(), layer = gsub("^.*/(.*).shp$", "\\1", files[1]))
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/robhodges17/GDrive/HWind", layer: "al12.2005_0825_21_00"
## with 25921 features and 7 fields
## Feature type: wkbPoint with 2 dimensions

columns.to.keep <- c("ID", "LONGITUDE", "LATITUDE", "SFC_SPD", "SFC_DIR")
spdf.data <- spdf.data[columns.to.keep]

for (i in 2:length(files)) {
    temp.data <- readOGR(dsn = getwd(), layer = gsub("^.*/(.*).shp$", "\\1", 
        files[i]), verbose = F)
    temp.data <- temp.data[columns.to.keep]
    spdf.data <- spRbind(spdf.data, temp.data)
}
## Error: Cannot open layer

Of interest here is the geo-reference wind speeds (“SFC_SPD”) in meters per second. To avert problems from overlap, we apply a tessellation of hexagons to the study area and identify the maximum wind speed observed per region.

hpt = spsample(spdf.data, type = "hexagonal", n = 20000)
hpg = HexPoints2SpatialPolygons(hpt)
Wind.hexid = overlay(spdf.data, hpg)
Wind.split = split(spdf.data@data, Wind.hexid)
names(Wind.split) = sapply(hpg@polygons, function(x) x@ID)[as.numeric(names(Wind.split))]
Wind.max = sapply(Wind.split, function(x) max(x$SFC_SPD))
Wind.max = data.frame(Wind.max)
Wind.spdf = SpatialPolygonsDataFrame(hpg[rownames(Wind.max)], Wind.max)

Now, we can plot a representation of the maximum windspeeds from all windfields from Hurricane Katrina.

box = bbox(Wind.spdf)
cl = map("worldHires", xlim = c(box[1], box[3]), ylim = c(box[2], box[4]), plot = FALSE)
clp1 = map2SpatialLines(cl)  #, proj4string=CRS(ll))
cl = map("state", xlim = c(box[1], box[3]), ylim = c(box[2], box[4]), interior = T, 
    boundary = F, plot = FALSE)
clp2 = map2SpatialLines(cl)  #, proj4string=CRS(ll))

spplot(Wind.spdf, "Wind.max", col = "transparent", at = c(17, 33, 43, 50, 58, 
    70, 100), lwd = 0.1, col.regions = c("darkslategray1", rev(heat_hcl(5))), 
    colorkey = list(space = "bottom", labels = list(at = c(17, 33, 43, 50, 58, 
        70, 100)), cex = 0.9), xlim = c(box[1] * 0.97, box[3]), ylim = c(box[2], 
        box[4] * 0.8), sub = list(label = expression(paste("Maximum observed wind speed [m", 
        s^-1, "]")), cex = 1), panel = function(x, y, ...) {
        panel.polygonsplot(x, y, ...)
        sp.lines(clp1, col = "black", lwd = 0.2)
        sp.lines(clp2, col = "black", lwd = 0.3)
    })

plot of chunk plotFields