Storm warning validity
Nam-Young Kang and James B. Elsner
Abstract
When facing the potential of a storm disaster the warning process and associated administrative activities across the western North Pacific are confounded by various tropical cyclone (TC) classifications possibly made worse by global warming. In this study, the influence global warming has had on the warning categories in the TC classifications from Japan Meteorological Agency, U.S. Joint Typhoon Warning Center and Saffir-Simpson’s hurricane intensity scale is quantitatively reviewed and compared over a 30-year period (1986–2015). TC classifications modeled with logistic regression show that the validity of current warning categories suffers from the warming environment. Cumulative proportions of storms falling into all warning categories are increasing. Enlarging proportions are correctly interpreted as an increasing threat. At the same time, the frequency of warnings is potentially making the warnings less effective. The results from this study may change as time goes by, but the feature of the response to a warming planet would remain as long as the warming continues. The findings here should help inform and improve the warning procedures toward more efficient preparedness against the increasingly often catastrophic storms.
Figure 1. P-P plot of the TC intensities from JMA and JTWC
# Sources
#-------------------------------------------------------------------------------
AB=read.table('Data/ABtrack2015.dat',T)
Year=c(1986,2015)
LMI=round(seq(40,150,10))
# best tracks
#-------------------------------------------------------------------------------
PLEVjm=NULL ; PLEVjt=NULL
for(i in Year[1]:Year[2])
{imsi=AB$lmw[AB$bst=='JM' & AB$yr==i & AB$mo>=1 & AB$mo<=12]
Fn=ecdf(-imsi) ; PLEVjm=cbind(PLEVjm, Fn(-LMI))
imsi=AB$lmw[AB$bst=='JT' & AB$yr==i & AB$mo>=1 & AB$mo<=12]
Fn=ecdf(-imsi) ; PLEVjt=cbind(PLEVjt, Fn(-LMI))
}
Mjm=NULL; for(i in 1:length(LMI)) Mjm=c(Mjm,mean(PLEVjm[i,]))
Mjt=NULL; for(i in 1:length(LMI)) Mjt=c(Mjt,mean(PLEVjt[i,]))
# Data
#-------------------------------------------------------------------------------
Mjm=NULL; for(i in 1:length(LMI)) Mjm=c(Mjm,mean(PLEVjm[i,]))
Mjt=NULL; for(i in 1:length(LMI)) Mjt=c(Mjt,mean(PLEVjt[i,]))
Data=data.frame(LMI=LMI,JMA=Mjm,JTWC=Mjt)
#Plot
#-------------------------------------------------------------------------------
par(mfrow=c(1,1),mar=c(5,5.5,3,3))
plot(-999,-999,xlim=c(1+0.03,-0.03),ylim=c(1+0.03,-0.03), xlab='',ylab='', axes=FALSE)
abline(a=0,b=1,lty=2)
points(Mjt,Mjm,pch=16,cex=2.3,col='#424242')
#axis
axis(1,at=seq(1,0,-0.1),labels=c('100','','80','','60','','40','','20','','0'),cex.axis=1.2)
axis(2,at=seq(1,0,-0.1),labels=c('100','','80','','60','','40','','20','','0'),cex.axis=1.2,las=1)
mtext('Cumulative proportion (%) in JTWC',side=1,line=3.2, cex=1.3)
mtext('Cumulative proportion (%) in JMA',side=2,line=3.2, cex=1.3)
textcol='#33339993'
text(Mjt[1:8],Mjm[1:8],paste(LMI[1:8],'kt'),cex=0.9,adj=1.4,col=textcol)
text(Mjt[9],Mjm[9]-0.03,paste(LMI[9],'kt'),cex=0.9,adj=1.3,col=textcol)
text(Mjt[10],Mjm[10]-0.05,paste(LMI[10],'kt'),cex=0.9,adj=0.8,col=textcol)
text(Mjt[11],Mjm[11]-0.05,paste(LMI[11],'kt'),cex=0.9,adj=0.4,col=textcol)
text(Mjt[12],Mjm[12]-0.05,paste(LMI[12],'kt'),cex=0.9,adj=0.2,col=textcol)
box()
Fig. 1. P-P plot of the LMIs at 10-kt intervals in JTWC and JMA best-track data. Each black dot represents the 30-yr (1986–2015) mean of annual probability levels. Probability level is calculated by the cumulative proportion inversely from the highest LMI. The dashed diagonal line draws the matching probability level of the TC intensities between JTWC and JMA.
Figure 2. Correlation coefficients of the TC intensities with time and GMSST
# Sources
#-------------------------------------------------------------------------------
gmsst=read.table('Data/gmsst_12m_1854to2015.dat',T)[133:162,1]
AB=read.table('Data/ABtrack2015.dat',T)
Year=c(1986,2015)
LMI=round(seq(40,150,10))
# best tracks
#-------------------------------------------------------------------------------
PLEVjm=NULL ; PLEVjt=NULL
for(i in Year[1]:Year[2])
{imsi=AB$lmw[AB$bst=='JM' & AB$yr==i & AB$mo>=1 & AB$mo<=12]
Fn=ecdf(-imsi) ; PLEVjm=cbind(PLEVjm, Fn(-LMI))
imsi=AB$lmw[AB$bst=='JT' & AB$yr==i & AB$mo>=1 & AB$mo<=12]
Fn=ecdf(-imsi) ; PLEVjt=cbind(PLEVjt, Fn(-LMI))
}
# Data
#-------------------------------------------------------------------------------
imsi=NULL;imsi2=NULL
for (i in 1:12) {ct=cor.test(1:30,PLEVjm[i,]);imsi=c(imsi,ct$estimate);imsi2=rbind(imsi2,ct$conf.int)} ; corr1=imsi ; conf.int1=imsi2## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zeroimsi=NULL;imsi2=NULL
for (i in 1:12) {ct=cor.test(gmsst,PLEVjm[i,]);imsi=c(imsi,ct$estimate);imsi2=rbind(imsi2,ct$conf.int)} ; corr2=imsi ;conf.int2=imsi2## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zeroimsi=NULL;imsi2=NULL
for (i in 1:12) {ct=cor.test(1:30,PLEVjt[i,]);imsi=c(imsi,ct$estimate);imsi2=rbind(imsi2,ct$conf.int)} ; corr3=imsi ; conf.int3=imsi2
imsi=NULL;imsi2=NULL
for (i in 1:12) {ct=cor.test(gmsst,PLEVjt[i,]);imsi=c(imsi,ct$estimate);imsi2=rbind(imsi2,ct$conf.int)} ; corr4=imsi ; conf.int4=imsi2
Sensitivity=data.frame(LMI=LMI, jm4time=corr1, jm4gmsst=corr2,jt4time=corr3, jt4gmsst=corr4)
Sensitivity_conf.int=data.frame(LMI=LMI, jm4time=conf.int1, jm4gmsst=conf.int2, jt4time=conf.int3, jt4gmsst=conf.int4)
# Plot
#-------------------------------------------------------------------------------
par(mfrow=c(1,2),mar=c(5,5.5,4,6))
# JMA 110kt => upper 5% level of 30-year LMIs
plot(-999,-999,xlim=c(-1,1),ylim=c(35,110),xaxs='i', xlab='',ylab='', axes=FALSE)
abline(v=0,lty=2,col='#333333')
polygon(cbind(c(Sensitivity_conf.int[1:8,2],Sensitivity_conf.int[8:1,3]),c(LMI[1:8],LMI[8:1])),border=NA,col='#33333333')
polygon(cbind(c(Sensitivity_conf.int[1:8,4],Sensitivity_conf.int[8:1,5]),c(LMI[1:8],LMI[8:1])),border=NA,col='#ff660033')
lines(Sensitivity[1:8,2],LMI[1:8],col='#333333',lwd=2) ; points(Sensitivity[1:8,2],LMI[1:8],col='#333333',pch=16,cex=1.8)
lines(Sensitivity[1:8,3],LMI[1:8],col='#ff6600',lwd=2) ; points(Sensitivity[1:8,3],LMI[1:8],col='#ff6600',pch=16,cex=1.8)
axis(1,at=seq(-1,1,0.2),labels=c('-1.0','','-0.6','','-0.2','','0.2','','0.6','','1.0'),cex.axis=1.2)
axis(2,at=LMI,labels=LMI,cex.axis=1.2,las=1)
axis(3,at=seq(-1,1,0.2),labels=rep('',11))
axis(4,at=LMI, labels=round(100*Data$JMA),cex.axis=1.2,las=1)
mtext('(a) JMA best track',side=3,line=1.5, cex=1.7, adj=-0.1)
mtext('Correlation coefficient',side=1,line=3.2, cex=1.5)
mtext('Threshold LMI (kt)',side=2,line=3.2, cex=1.5)
text(par("usr")[2] + 0.32, 96.7, srt=-90, adj = 0, labels = 'Cumulative proportion (%)', cex=1.5,xpd = TRUE)
box()
# JTWC 150kt => upper 3% level of 30-year LMIs
plot(-999,-999,xlim=c(-1,1),ylim=c(35,150),xaxs='i',xlab='',ylab='', axes=FALSE)
abline(v=0,lty=2,col='#333333')
polygon(cbind(c(Sensitivity_conf.int[1:12,6],Sensitivity_conf.int[12:1,7]),c(LMI[1:12],LMI[12:1])),border=NA,col='#33333333')
polygon(cbind(c(Sensitivity_conf.int[1:12,8],Sensitivity_conf.int[12:1,9]),c(LMI[1:12],LMI[12:1])),border=NA,col='#ff660033')
lines(Sensitivity[1:12,4],LMI[1:12],col='#333333',lwd=2) ; points(Sensitivity[1:12,4],LMI[1:12],col='#333333',pch=16,cex=1.8)
lines(Sensitivity[1:12,5],LMI[1:12],col='#ff6600',lwd=2) ; points(Sensitivity[1:12,5],LMI[1:12],col='#ff6600',pch=16,cex=1.8)
axis(1,at=seq(-1,1,0.2),labels=c('-1.0','','-0.6','','-0.2','','0.2','','0.6','','1.0'),cex.axis=1.2)
axis(2,at=LMI,labels=LMI,cex.axis=1.2,las=1)
axis(3,at=seq(-1,1,0.2),labels=rep('',11))
axis(4,at=LMI, labels=round(100*Data$JTWC),cex.axis=1.2,las=1)
mtext('(b) JTWC best track',side=3,line=1.5, cex=1.7, adj=-0.1)
mtext('Correlation coefficient',side=1,line=3.2, cex=1.5)
mtext('Threshold LMI (kt)',side=2,line=3.2, cex=1.5)
text(par("usr")[2] + 0.32, 130, srt=-90, adj = 0, labels = 'Cumulative proportion (%)', cex=1.5,xpd = TRUE)
box()
#-Ref 1. GMSST change over the 30 years (1986--1985)
#############################################################
model=lm(gmsst~c(1:30))
summary(model)##
## Call:
## lm(formula = gmsst ~ c(1:30))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.091214 -0.051742 0.009687 0.038627 0.147444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.156232 0.024083 753.89 < 2e-16 ***
## c(1:30) 0.011002 0.001357 8.11 7.89e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06431 on 28 degrees of freedom
## Multiple R-squared: 0.7014, Adjusted R-squared: 0.6907
## F-statistic: 65.77 on 1 and 28 DF, p-value: 7.886e-09summary(model)$coef[2,1]*30## [1] 0.3300504summary(model)$coef[2,2]*30## [1] 0.04069743# From a linear perspective, the increase of global ocean warmth is significant as $+$0.33 $\pm$ 0.041 (s.e.) C/30 yr.Fig. 2. Correlation coefficients of the annual probability levels of LMI with time (linked black dots) and GMSST (linked orange dots), respectively from (a) JMA and (b) JTWC best-track data. 95% confidence intervals are shaded in gray and orange colors for time and GMSST, for each. The average of the annual probability levels is labeled in the right ordinate.
Figure 3. Global warming influence on the warning categories in the TC classifications
library(akima)## Warning: package 'akima' was built under R version 3.3.2# Sources
#--------------------------------------------------------------------------------
gmsst=read.table('Data/gmsst_12m_1854to2015.dat',T)[133:162,1]
AB=read.table('Data/ABtrack2015.dat',T)
Year=c(1986,2015) ; time=1986:2015
prob=round(seq(0.01,0.99,0.01),2)
CLASSjm=c(64,85,105) # JMA scale
CLASSjt=c(64,130) # JTWC sclale
CLASSss=c(64,83,96,113,137) # Saffir-Siimpson scale
# Prediction values of GMSST
#-------------------------------------------------------------------------------
pred_gmsst= predict(lm(gmsst~c(Year[1]:Year[2])))
# Annual variation of the probability level of threshold LMI
#-------------------------------------------------------------------------------
PLEVjm=NULL ; PLEVkm=NULL ; PLEVjt=NULL ; PLEVss=NULL
for(i in Year[1]:Year[2])
{imsi=AB$lmw[AB$bst=='JM' & AB$yr==i & AB$mo>=1 & AB$mo<=12] #JMA best track
Fn=ecdf(-imsi)
PLEVjm=rbind(PLEVjm, Fn(-CLASSjm))
imsi=AB$lmw[AB$bst=='JT' & AB$yr==i & AB$mo>=1 & AB$mo<=12] #JTWC best track
Fn=ecdf(-imsi)
PLEVjt=rbind(PLEVjt, Fn(-CLASSjt))
PLEVss=rbind(PLEVss, Fn(-CLASSss))
}
#Model fitting and data
#-------------------------------------------------------------------------------
PLEVall=cbind(PLEVjm,PLEVjt,PLEVss)
CLASS=c('jm64kt', 'jm85kt', 'jm105kt'
,'jt64kt', 'jt130kt'
,'ss64kt', 'ss83kt', 'ss96kt','ss113kt','ss137kt')
pred_p_by_time=NULL ; pred_p_by_gmsst= NULL
p.val_p_by_time=NULL ; p.val_p_by_gmsst=NULL
pred_p_mean= NULL ; pred_p_mean2=NULL
for (i in 1:10)
{
Success= round(1000*PLEVall[,i]) ; Failure=1000-Success ; Trials = cbind(Success, Failure)
p=Success/(Success+Failure) ; Trials = cbind(Success, Failure)
model = glm(Trials ~ time, family = binomial(link="logit")) # by time
model2 = glm(Trials ~ gmsst, family = binomial(link="logit")) # by gmsst
pred_gmsst=predict(lm(gmsst~time))
p.val_p_by_time=c(p.val_p_by_time, summary(model)$coef[2,4])
p.val_p_by_gmsst=c(p.val_p_by_gmsst, summary(model2)$coef[2,4])
pred_p_by_time= cbind(pred_p_by_time,predict(model,type="response"))
pred_p_by_gmsst= cbind(pred_p_by_gmsst,predict(model2,data.frame(gmsst=pred_gmsst),type="response"))
pred_p_mean= c(pred_p_mean,mean(predict(model2,type="response")))
pred_p_mean2= c(pred_p_mean2,mean(predict(model2,data.frame(gmsst=pred_gmsst),type="response")))
}
PLEVall=as.data.frame(PLEVall) ; names(PLEVall)=CLASS
pred_p_by_time=as.data.frame(pred_p_by_time) ; names(pred_p_by_time)=CLASS
pred_p_by_gmsst=as.data.frame(pred_p_by_gmsst) ; names(pred_p_by_gmsst)=CLASS
p.val=data.frame(CLASS=CLASS,time=round(p.val_p_by_time,3),gmsst=round(p.val_p_by_gmsst,3))
p.val## CLASS time gmsst
## 1 jm64kt 0.000 0.000
## 2 jm85kt 0.000 0.000
## 3 jm105kt 0.000 0.000
## 4 jt64kt 0.000 0.297
## 5 jt130kt 0.000 0.000
## 6 ss64kt 0.000 0.297
## 7 ss83kt 0.834 0.000
## 8 ss96kt 0.000 0.000
## 9 ss113kt 0.000 0.000
## 10 ss137kt 0.000 0.000# P-value
#All GLM results shows significant response to the increasing global ocean warmth at 95% confidence level, except for the thresholds of 'Typhoon' in JTWC and 'Hurricane category 1' in SS.
#Plot
#-------------------------------------------------------------------------------
par(mfrow=c(1,1),mar=c(5,5,3,2.5))
b=c(1,30,59,88) #border
plot(-999,-999,xlim=c(1,88),ylim=c(1,0),xlab='',ylab='',xaxs='i',yaxs='i', axes=FALSE)
collist=c('#99663305','#99663315','#99663340', '#99663360',
# '#ff990005','#ff990011','#ff990022','#ff990055',
'#66660011','#66660033','#66660055',
'#00009900','#00009905','#00009911','#00009922','#00009922','#00009922')
#JM
polygon(c(b[1],b[1],b[2],b[2],b[1]),c(0,1,1,0,0),border=NA,col=collist[1])
for(i in 1:3) polygon(c(b[1],c(b[1]:b[2]),b[2],b[1]),c(0,pred_p_by_gmsst[,i],0,0),border=NA,col=collist[i+1])
for(i in 1:3) lines(b[1]:b[2], pred_p_by_gmsst[,i], col=1,lwd=1)
#for(i in 1:3) lines(b[1]:b[2], pred_p_by_time[,i], col='#00640040',lwd=3)
abline(v=b[2])
#JT
polygon(c(b[2],b[2],b[3],b[3],b[2]),c(0,1,1,0,0),border=NA,col=collist[5])
for(i in 4:5) polygon(c(b[2],c(b[2]:b[3]),b[3],b[2]),c(0,pred_p_by_gmsst[,i],0,0),border=NA,col=collist[i+2])
for(i in 4:5) lines(b[2]:b[3], pred_p_by_gmsst[,i], col=1,lwd=1)
#for(i in 4:5) lines(b[2]:b[3], pred_p_by_time[,i],col='#00640040',lwd=3)
abline(v=b[3])
#SS
polygon(c(b[3],b[3],b[4],b[4],b[3]),c(0,1,1,0,0),border=NA,col=collist[9])
for(i in 6:10) polygon(c(b[3],c(b[3]:b[4]),b[4],b[3]),c(0,pred_p_by_gmsst[,i],0,0),border=NA,col=collist[i+3])
for(i in 6:10) lines(b[3]:b[4], pred_p_by_gmsst[,i], col=1,lwd=1)
#for(i in 6:10) lines(b[3]:b[4], pred_p_by_time[,i], col='#00640040',lwd=3)
abline(v=b[4])
#axis
mtext('Each 30-year period (1986-2015)',side=1,line=3.3, cex=1.5)
mtext('Cumulative proportion (%)',side=2,line=3.3, cex=1.5)
axis(1,at=c(5,15,25,34,44,54,63,73,83),labels=rep(c(1990,2000,2010),3),cex.axis=0.9)
#axis(2,at=seq(1,0,-0.1),labels=c(100,'',80,'',60,'',40,'',20,'',0),cex.axis=1.4,las=1)
axis(2,at=seq(1,0,-0.2),labels=c(100,80,60,40,20,0),cex.axis=1.4,las=1)
#axis(4,at=seq(1,0,-0.2),labels=rep('',6),cex.axis=1.4,las=1)
text(15,-0.04,labels='JMA',cex=1.4,col=1,xpd=TRUE)
text(44,-0.04,labels='JTWC',cex=1.4,col=1,xpd=TRUE)
text(73,-0.04,labels='SS',cex=1.4,col=1,xpd=TRUE)
#note proportion values
val1=as.numeric(pred_p_by_gmsst[1,])
val2=as.numeric(pred_p_by_gmsst[30,])
text(cbind(c(b[1],b[1],b[1],b[2],b[2],b[3],b[3],b[3],b[3],b[3])-2.3,val1-0.005),labels=round(100*val1),cex=0.8,col='dark blue',xpd=TRUE)
text(cbind(c(b[2],b[2],b[2],b[3],b[3],b[4],b[4],b[4],b[4],b[4])+2.3,val2+0.005),labels=round(100*val2),cex=0.8,col='#EE7600',xpd=TRUE)
#Category names
text(cbind(c(15,15,15,15),c(0.03,0.19,0.45,0.83)),labels=c('Violent typhoon','Very strong typhoon','Strong typhoon','Typhoon'),
cex=1.0,col=colors()[188],xpd=TRUE) # JMA
text(cbind(c(44,44,44),c(0.08,0.45,0.83)),labels=c('Super typhoon','Typhoon','Tropical storm'),
cex=1.0,col=colors()[188],xpd=TRUE) #JTWC
text(cbind(c(73,73,73,73,73,73),c(0.05,0.19,0.33,0.43,0.58,0.83)),labels=c('Cat. 5','Cat. 4','Cat. 3', 'Cat. 2', 'Cat. 1','Tropical storm'),
cex=1.0,col=colors()[188],xpd=TRUE) #JTWC
box()
Fig. 3. Response of the probability level of storm warning categories in the TC classifications to global warming over the 30 years (1986–2015). Binomial logistic regression is employed as a GLM to model the probability level of the threshold LMIs. Based on a linear assumption of global warming, the modeled GMSST by OLS method is used for the single explanatory variable. All GLM results show significant responses to GMSST increase at 95% confidence level, except for the thresholds of ‘Typhoon’ in JTWC and ‘Hurricane category 1’ in SS. The storm warning categories for JTWC and SS are referenced by JTWC best-track data with 1-min average of the winds, which is different from JMA with 10-min average. Probability level indicated by the cumulative proportion, however, enables the comparison between the best-track data sets following different observational procedures. All values are rounded off to the nearest whole number.
Table 1. 30-yr mean of annual cumulative proportions
p_mean=NULL; for (i in 1:10) p_mean=c(p_mean, mean(PLEVall[,i]))
proportion=round(100*p_mean) # original mean vlaue
cbind(CLASS, proportion)## CLASS proportion
## [1,] "jm64kt" "55"
## [2,] "jm85kt" "31"
## [3,] "jm105kt" "7"
## [4,] "jt64kt" "64"
## [5,] "jt130kt" "18"
## [6,] "ss64kt" "64"
## [7,] "ss83kt" "46"
## [8,] "ss96kt" "36"
## [9,] "ss113kt" "29"
## [10,] "ss137kt" "12"Tab. 1. TC classifications for JMA, JTWC, and SS. Each 30-yr mean of the cumulative proportions at threshold LMI is denoted in parentheses. JMA is referenced by JMA best-track data following 10-min period for averaging winds, while JTWC and SS by JTWC best-track data following 1-min period. Values are rounded off to the nearest whole number.