# This is an R-script.
# The R-environment, under which this script can be run, is freely available from CRAN:
# http://cran.r-project.org
#
# R.E. Benestad, Oslo, 15.04.2005
#
# in R, copy this script to working directory (ensure that this is the same as the working directory - type '? getwd' 
# or '?setwd') as 'hurricanes.R', then type: 'source("hurricanes.R")'
# R reads the data from source over URL...

rm(list=ls())   # Clear memory...
update <- TRUE

# Empirical ranking method:
#
# A formula for estimating the cumulative probability P corresponding to rank m
# Original references:
#  Jenkinson, A.F., 1977, U.K. Met.Office Synoptic Clim. Branch Memo 58
#  Beard, L.R., 1943, Trans. Amer. Meteor. Soc. Civ. Eng., 108, 1110-1160
#  Chegodaev, N.N., 1953 (in Russian) State Rail Transport Publishing House.
# Reference:
# Folland, C. and Anderson, C. (2002), J. Clim. 15, 2954-2960, equation (1)
#
# The formula applies to any series of continously distributed measurements
# provided it can be regarded as a random sample. Simple and accurate alternative
# to qbeta().
#
# Coded for R by R.E. Benestad, met.no, 27.06.2003.

empirical.ranking <- function(x) {
  N <- length(x)
  m <- rank(x)
  P <- 100*(m-0.31)/(N+0.38)
  sort <- order(x)
  P <- P[sort]
  y <- as.numeric(x[sort])
  results <- data.frame(x=y, P=P)
  results
}

# Function selects the highest storm category
max.category <- function(category,ii) {
  Categories <- c("Hurricane-Category-5","Hurricane-Category-4","Hurricane-Category-3","Hurricane-Category-2",
                  "Hurricane-Category-1","Tropical-Storm","Subtropical-Storm","Extratropical-Storm","other")
  category[!ii] <- "other"
  if (sum(ii) >= 1) {
  i <- 0
      iii <- rep(FALSE,length(ii))
    while ((sum(iii)==0) & (i < length(Categories))) {
      i <- i + 1
      iii <- iii | is.element(category,Categories[i])
    }
    ii <- iii & ii
  }
  ii
} 

Filename <- "http://www.aoml.noaa.gov/hrd/hurdat/hurdatTAB.txt"

head <- c("Storm.Name","Storm.No","Year","Month","Day","Hour",
          "Latitude","Longitude", "Storm.Direction","Storm.Speed.mph",
          "Storm.Speed.kph","Wind.Speed.mph","Wind.Speed.kph",
          "Pressure","Storm.category")

# If the data are not allready retrieved, get them from the URL.

if ((file.exists("hurricanes.Rdata")) & (!update)) {
  print("Reading locally stored data")
  ac <- read.table("hurricanes.txt",header=TRUE) 
} else {
  print(paste("Reading from",Filename))
  ac <- read.table(Filename,header=FALSE,sep="\t",skip=1)
  colnames(ac) <- head
  write.table(file="hurricanes.txt",ac)
}


levs <- levels(ac$Storm.category)
ac$Storm.category <- as.character(ac$Storm.category)
ac$Pressure <- as.numeric(as.character(ac$Pressure))

# A re-definition of category names is required due to inconsistent labelling
iES <- grep("Extratropical Storm",ac$Storm.category)
iSS <- grep("Subtropical Storm",ac$Storm.category)
iTS <- sort(c(grep("TropicalStorm",ac$Storm.category),
              grep("Tropical Storm",ac$Storm.category)))
iH1 <- grep("1",ac$Storm.category)
iH2 <- sort(c(grep("2",ac$Storm.category),
              grep("Hurricane - Category 2",ac$Storm.category)))
iH3 <- grep("3",ac$Storm.category)
iH4 <- grep("4",ac$Storm.category)
iH5 <- grep("5",ac$Storm.category)

ac$Storm.category <- as.character(ac$Storm.category)
ac$Storm.category[iES] <- "Extratropical-Storm"
ac$Storm.category[iSS] <- "Subtropical-Storm"
ac$Storm.category[iTS] <- "Tropical-Storm"
ac$Storm.category[iH1] <- "Hurricane-Category-1"
ac$Storm.category[iH2] <- "Hurricane-Category-2"
ac$Storm.category[iH3] <- "Hurricane-Category-3"
ac$Storm.category[iH4] <- "Hurricane-Category-4"
ac$Storm.category[iH5] <- "Hurricane-Category-5"

idefined <- sort(c(iES,iSS,iTS,iH1,iH2,iH3,iH4,iH5))

# The events not picked up are defined as 'other'
all <- 1:length(ac$Storm.category)
iother <- all[!is.element(all,idefined)]
print(table(ac$Storm.category[iother]))
ac$Storm.category[iother] <- "other"

storm <- c(NA); year <- storm; wind <- storm; psl <- storm; lat <- storm; lon <- storm; storm.name <- storm

storm.cat <- rownames(table(ac$Storm.category))

tot.storms <- diff(ac$Storm.No); tot.storms[tot.storms < 0] <- 1; 
tot.storms.1 <- sum(tot.storms) - 1

storm <- rep(NA,tot.storms.1); wind <- storm; psl <- storm; 
lon <- storm; lat <- storm; year <- storm; storm.name <- storm

print("record#   storm#   storm-name    storm-category")
i <- 1; istorm <- 1

while ((i <= length(ac$Storm.No)) & (istorm <= tot.storms.1)) {
  i0 <- is.element(ac$Storm.No,ac$Storm.No[i]) & is.element(ac$Year,ac$Year[i])
  ii <- i0 & max.category(ac$Storm.category,i0)
  iii <- ii & is.element(ac$Wind.Speed.kph, max(ac$Wind.Speed.kph[ii]))
  if (sum(iii) > 1) iii <- iii & is.element(ac$Pressure,min(ac$Pressure[iii]))
  storm[istorm] <-as.character(ac$Storm.category[iii][1])
  wind[istorm] <- as.numeric(as.character(ac$Wind.Speed.kph[iii][1]))
  psl[istorm] <-  as.numeric(as.character(ac$Pressure[iii][1]))
  lon[istorm] <-  as.numeric(as.character(ac$Longitude[iii][1]))
  lat[istorm] <-  as.numeric(as.character(ac$Latitude[iii][1]))
  year[istorm] <- ac$Year[iii][1]
  storm.name[istorm] <- as.character(ac$Storm.Name[iii][1])
  print(paste(i,istorm,year[istorm],storm.name[istorm],storm[istorm],wind[istorm]))
  i <- i + sum(i0); istorm <- istorm + 1
}
istorm <- istorm - 1
print(paste("Check: i=",i,"=",length(ac$Storm.No),"? istorm:",istorm,"=",
            tot.storms.1,"?"))
psl[psl < 800] <- NA
years <- as.numeric(rownames(table(year)))
ny <- length(years)

maxwind <- rep(NA,ny); minpsl <- maxwind; 
for (i in 1:ny) { 
  ii <- is.element(year,years[i])
  maxwind[i] <- max(wind[ii],na.rm=TRUE)
  minpsl[i] <- min(psl[ii],na.rm=TRUE)
}
minpsl[!is.finite(minpsl)] <- NA; maxwind[!is.finite(maxwind)] <- NA

plot(years,maxwind,pch=20,cex=1.7,ylim=c(50,350),xlim=c(1850,2020),
     ylab="Wind speed (kph)", 
     main="Maximum hurricane wind speed/min pressure",
     sub="data: http://www.aoml.noaa.gov/hrd/hurdat/hurdatTAB.txt")
grid()

p.min <- 2.75*(minpsl-min(minpsl,na.rm=TRUE))
lines(years+0.5,p.min+0.5,col="grey60",border="grey",lwd=5)
lines(years,p.min,col="grey60",lwd=5)
recent <- years > 1950
#abline(lm(p.min[recent] ~ years[recent]),lty=2,col="steelblue")

for (i in seq(min(minpsl,na.rm=TRUE),max(minpsl,na.rm=TRUE),by=10)) {
  lines(c(2017,2020),rep(2.5*(i-min(minpsl,na.rm=TRUE)),2),lwd=2)
  lines(c(2017,2020),rep(2.5*(i-min(minpsl,na.rm=TRUE)),2),col="grey60")
  text(2012,2.5*(i-min(minpsl,na.rm=TRUE)),i,cex=0.7)
}

abline(lm(maxwind ~ years),lty=2)
points(years,maxwind,pch=20,cex=1.7)
points(years+0.5,maxwind+0.5,pch=20,cex=1.5,col="grey20")
points(years+0.5,maxwind+0.5,pch=20,cex=1.0,col="grey40")
points(years+0.5,maxwind+0.5,pch=20,cex=0.8,col="grey60")

mtext("Sea level pressure (hPa)",font=2,side=4,col="grey60",cex=1.5)

text(1990,350,"Benestad (2005)",cex=0.8,col="grey90")

dev2bitmap("hurricane_speed.png",type="png256",res=200)

x11()


io0 <-  grep("other",storm)
is1 <- c(grep("Tropical Storm",storm),grep("TropicalStorm",storm),
         grep("Tropical-Storm",storm),grep("Subtropical-Storm",storm),
         grep("Extratropical-Storm",storm))
ih1 <- grep("1",storm)
ih2 <- grep("2",storm)
ih3 <- grep("3",storm)
ih4 <- grep("4",storm)
ih5 <- grep("5",storm)
i0 <- -sort(c(is1,ih1,ih2,ih3,ih4,ih5))
print(table(storm))

yrs <- seq(1850,2005,by=5)

other <- hist(year[io0],breaks=yrs)
s1 <- hist(year[is1],breaks=yrs)
h1 <- hist(year[ih1],breaks=yrs)
h2 <- hist(year[ih2],breaks=yrs)
h3 <- hist(year[ih3],breaks=yrs)
h4 <- hist(year[ih4],breaks=yrs)
h5 <- hist(year[ih5],breaks=yrs)

all <- s1$counts + h1$counts + h2$counts + h3$counts + h4$counts + h5$counts + other$counts
pentads <- s1$mids

plot(pentads,all,type="n",main="Most severe hurricane category",xlab="count",ylim=c(0,70),
     sub="data: http://www.aoml.noaa.gov/hrd/hurdat/hurdatTAB.txt")
grid()
for (i in 1:length(pentads)) {
  polygon(c(-4,4,4,-4,-4)/2+pentads[i],c(0,0,all[i],all[i],0),col="grey60",border="grey30",lwd=2)
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,h5$counts[i],h5$counts[i],0),col="red",lwd=2)
  i0 <- h5$counts[i]
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,h4$counts[i],h4$counts[i],0)+i0,col="blue")
  i0 <- i0 + h4$counts[i]
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,h3$counts[i],h3$counts[i],0)+i0,col="steelblue")
  i0 <- i0 + h3$counts[i]
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,h2$counts[i],h2$counts[i],0)+i0,col="pink")
  i0 <- i0 + h2$counts[i]
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,h1$counts[i],h1$counts[i],0)+i0,col="lightblue")
  i0 <- i0 + h1$counts[i]
  polygon(c(-3,3,3,-3,-3)/2+pentads[i],c(0,0,s1$counts[i],s1$counts[i],0)+i0,col="grey90",border="grey30",lwd=2)

}


lines(pentads,40*(h3$counts + h4$counts + h5$counts)/(s1$counts + h1$counts + h2$counts),lwd=5,col="white")
lines(pentads,40*(h3$counts + h4$counts + h5$counts)/(s1$counts + h1$counts + h2$counts),lwd=2,col="grey60")

mtext("# major hurricanes/# weaker cyclones (x 40)",font=2,side=4,col="grey60",cex=1.25)

text(1990,70,"Benestad (2005)",cex=0.8,col="grey90")

legend(1850,70,c("Tropical Storm        ","Category 1","Category 2","Category 3","Category 4","Category 5"),
       col=c("grey90","lightblue","pink","steelblue","blue","red"),lwd=5,cex=0.8)
dev2bitmap("hurricane_category.png",type="png256",res=200)

#print(cbind(storm.name[is.element(year,2004)],storm[is.element(year,2004)]))

# Check sums / control...

tot.storms.2 <- sum(all)

x11()
n.year <- table(year)
edf <- empirical.ranking(n.year)
plot(edf$x,edf$P,type="l",lwd=4,main="empirical distribution function",
     sub="Jenkinson (1977)",xlab="Number of storms per year",ylab="Probability")
y <- as.numeric(attributes(edf)$row.names)
i2004 <- is.element(y,2004)
lines(rep(edf$x[i2004],2),c(0,edf$P[i2004]),lty=2,col="red")
lines(c(min(edf$x),edf$x[i2004]),rep(edf$P[i2004],2),lty=3,col="red")
text(edf$x[i2004],edf$P[i2004]+5,paste("2004: n>",round(edf$P[i2004]),"%"))
text(19,0,"Benestad (2005)",cex=0.8,col="grey90")
dev2bitmap("hurricane_edf.png",type="png256",res=200)

print(c(tot.storms.1,tot.storms.2,sum(table(storm)),istorm))
if (tot.storms.1 != tot.storms.2) stop("Check-sum error: total 1")
if (tot.storms.1 != sum(table(storm))) stop("Check-sum error: total 2")
if (tot.storms.1 != istorm) stop("Check-sum error: total 3")

if (length(is1) != sum(s1$counts)) stop("Check-sum error: S1")
if (length(ih1) != sum(h1$counts)) stop("Check-sum error: H1")
if (length(ih2) != sum(h2$counts)) stop("Check-sum error: H2")
if (length(ih3) != sum(h3$counts)) stop("Check-sum error: H3")
if (length(ih4) != sum(h4$counts)) stop("Check-sum error: H4")
if (length(ih5) != sum(h5$counts)) stop("Check-sum error: H5")

print("--------------------------------------------")
print(" Check-sums OK! ")