Code by Scott Miles @ http://resilscience.com/scott-miles/

In [81]:
from pandas import *
import matplotlib.pyplot as plt
import datetime
from IPython.core.display import Image, HTML
from pandas.core.datetools import *
from numpy import arange,array,ones#,random,linalg
from scipy import stats
figsize(12,10)
pandas.__version__
#matplotlib.__version__
Out[81]:
'0.13.1'
In [82]:
pandas.set_option('display.height', 500)
pandas.set_option('display.max_rows', 500)
pandas.set_option('display.max_columns', 500)
pandas.set_option('display.width', 500)
pandas.set_option('display.mpl_style', False)
In [83]:
xkcd()
matplotlib.matplotlib_fname()
rcParams['lines.linewidth'] = 3.0
rcParams['legend.shadow'] = True
rcParams['legend.fancybox'] = True
rcParams['legend.borderpad'] = 0.25

Import EMDAT data from Excel files

EMDAT data obtained here: http://www.emdat.be

EMDAT defines a disaster occurrence as

  • 10 or more people reported killed and/or
  • 100 or more people reported affected and/or
  • Call for international assistance/declaration of a state of emergency

In EMDAT "storms" include tropical storms, winter storms, and local/convective storms

In [84]:
eq_data = read_excel('eq_disasters_by_year&region_dollarss.xlsx', 0, index_col=None, parse_dates=[1])
eq_data.drop(['id','Affected','Impacted'],axis=1,inplace=True)
eq_data.head()
---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-84-7a968d18e2c8> in <module>()
----> 1 eq_data = read_excel('eq_disasters_by_year&region_dollarss.xlsx', 0, index_col=None, parse_dates=[1])
      2 eq_data.drop(['id','Affected','Impacted'],axis=1,inplace=True)
      3 eq_data.head()

/Users/geomando/Library/Enthought/Canopy_32bit/User/lib/python2.7/site-packages/pandas/io/excel.pyc in read_excel(io, sheetname, **kwds)
    101     engine = kwds.pop('engine', None)
    102 
--> 103     return ExcelFile(io, engine=engine).parse(sheetname=sheetname, **kwds)
    104 
    105 

/Users/geomando/Library/Enthought/Canopy_32bit/User/lib/python2.7/site-packages/pandas/io/excel.pyc in __init__(self, io, **kwds)
    134 
    135         if isinstance(io, compat.string_types):
--> 136             self.book = xlrd.open_workbook(io)
    137         elif engine == "xlrd" and isinstance(io, xlrd.Book):
    138             self.book = io

/Users/geomando/Library/Enthought/Canopy_32bit/User/lib/python2.7/site-packages/xlrd/__init__.pyc in open_workbook(filename, logfile, verbosity, use_mmap, file_contents, encoding_override, formatting_info, on_demand, ragged_rows)
    392         peek = file_contents[:peeksz]
    393     else:
--> 394         f = open(filename, "rb")
    395         peek = f.read(peeksz)
    396         f.close()

IOError: [Errno 2] No such file or directory: 'eq_disasters_by_year&region_dollarss.xlsx'
In [ ]:
storm_data = read_excel('storm_disasters_by_year&region_dollarss.xlsx', 0, index_col=None, parse_dates=[1])
storm_data.drop(['id','Affected','Impacted'],axis=1,inplace=True)
#storm_data.drop(['id','Affected','Impacted', 'Homeless','Killed','Injured'],axis=1,inplace=True)
storm_data.head()

General earthquake statistics

In [ ]:
eq_dmg = eq_data.ix[eq_data['Damage'] > 1000000] # Number over $1billion
print "Number of earthquake disasters over $1billion since 1900 = ", len(eq_dmg)
eq_kill = eq_data.ix[eq_data['Killed'] > 10000] # Number over 10,000
print "Number of earthquake disasters with over 10,000 dead since 1900 = ", len(eq_kill)
eq_inj = eq_data.ix[eq_data['Injured'] > 10000] # Number over 10,000
print "Number of earthquake disasters with over 10,000 injured since 1900 = ", len(eq_inj)
eq_home = eq_data.ix[eq_data['Homeless'] > 10000] # Number over 10,000
print "Number of earthquake disasters with over 10,000 homeless since 1900 = ", len(eq_home)

General storm statistics

In [ ]:
storm_dmg = storm_data.ix[storm_data['Damage'] > 1000000] # Number over $1billion
print "Number of storm disasters over $1billion since 1900 = ", len(storm_dmg)
storm_kill = storm_data.ix[storm_data['Killed'] > 10000] # Number over 10,000
print "Number of storm disasters with over 10,000 dead since 1900 = ", len(storm_kill)
storm_inj = storm_data.ix[storm_data['Injured'] > 10000] # Number over 10,000
print "Number of storm disasters with over 10,000 injured since 1900 = ", len(storm_inj)
storm_home = storm_data.ix[storm_data['Homeless'] > 10000] # Number over 10,000
print "Number of storm disasters with over 10,000 homeless since 1900 = ", len(storm_home)

Regional earthquake statistics

In [85]:
eq_region = eq_data.groupby('Region')
eq_region_max = eq_region.max()

eq_max_occur_region = eq_region_max['Occurrence'].idxmax()
eq_max_occur_count = eq_region['Occurrence'].max().max()
print 'Region with maximum earthquake disaster occurrences:', eq_max_occur_region, '=', eq_max_occur_count

eq_max_kill_region = eq_region_max['Killed'].idxmax()
eq_max_kill_count = eq_region['Killed'].max().max()
print 'Region with maximum earthquake disaster deaths:', eq_max_kill_region, '=', eq_max_kill_count

eq_max_inj_region = eq_region_max['Injured'].idxmax()
eq_max_inj_count = eq_region['Injured'].max().max()
print 'Region with maximum earthquake disaster injuries:', eq_max_inj_region, '=', eq_max_inj_count

eq_max_home_region = eq_region_max['Homeless'].idxmax()
eq_max_home_count = eq_region['Homeless'].max().max()
print 'Region with maximum earthquake disaster homeless population:', eq_max_home_region, '=', eq_max_home_count

eq_max_dmg_region = eq_region_max['Damage'].idxmax()
eq_max_dmg_count = eq_region['Damage'].max().max()
print 'Region with maximum earthquake disaster damage ($):', eq_max_dmg_region, '=', eq_max_dmg_count
Region with maximum earthquake disaster occurrences: Eastern Asia = 14
Region with maximum earthquake disaster deaths: Eastern Asia = 242000
Region with maximum earthquake disaster injuries: Eastern Asia = 369060
Region with maximum earthquake disaster homeless population: Southern Asia = 5150000
Region with maximum earthquake disaster damage ($): Eastern Asia = 210346400000

Regional storm statistics

In [86]:
storm_region = storm_data.groupby('Region')
storm_region_max = storm_region.max()

storm_max_occur_region = storm_region_max['Occurrence'].idxmax()
storm_max_occur_count = storm_region['Occurrence'].max().max()
print 'Region with maximum storm disaster occurrences:', storm_max_occur_region, '=', storm_max_occur_count

storm_max_kill_region = storm_region_max['Killed'].idxmax()
storm_max_kill_count = storm_region['Killed'].max().max()
print 'Region with maximum storm disaster deaths:', storm_max_kill_region, '=', storm_max_kill_count

storm_max_inj_region = storm_region_max['Injured'].idxmax()
storm_max_inj_count = storm_region['Injured'].max().max()
print 'Region with maximum storm disaster injuries:', storm_max_inj_region, '=', storm_max_inj_count

storm_max_home_region = storm_region_max['Homeless'].idxmax()
storm_max_home_count = storm_region['Homeless'].max().max()
print 'Region with maximum storm disaster homeless population:', storm_max_home_region, '=', storm_max_home_count

storm_max_dmg_region = storm_region_max['Damage'].idxmax()
storm_max_dmg_count = storm_region['Damage'].max().max()
print 'Region with maximum storm disaster damage ($):', storm_max_dmg_region, '=', storm_max_dmg_count
Region with maximum storm disaster occurrences: Western Europe = 41
Region with maximum storm disaster deaths: Southern Asia = 300317
Region with maximum storm disaster injuries: Southern Asia = 600000
Region with maximum storm disaster homeless population: Eastern Asia = 11002376
Region with maximum storm disaster damage ($): Northern America = 158230000000

Group data by year

In [87]:
eq_yr = eq_data.groupby('Year')
eq_yr_sum = eq_yr.sum()

storm_yr = storm_data.groupby('Year')
storm_yr_sum = storm_yr.sum()

#storm_eq_yr_sum = concat([eq_yr_sum, storm_yr_sum])
storm_eq_yr_sum = concat([eq_yr_sum.rename(columns={'Occurrence': 'Occurrence (Earthquake)', 
                        'Killed': 'Killed (Earthquake)', 
                        'Injured': 'Injured (Earthquake)',
                        'Damage': 'Damage (Earthquake)',
                        'Homeless': 'Homeless (Earthquake)'}),
                         storm_yr_sum.rename(columns={'Occurrence': 'Occurrence (Storm)', 
                        'Killed': 'Killed (Storm)', 
                        'Injured': 'Injured (Storm)',
                        'Damage': 'Damage (Storm)',
                        'Homeless': 'Homeless (Storm)'})], axis=1)

Monotonic correlation (Spearman) between occurrence and losses

In [88]:
print 'Earthquake Spearman correlation = ', eq_yr_sum['Occurrence'].corr(eq_yr_sum['Damage'], method='spearman')
print 'Storm Spearman correlation = ', storm_yr_sum['Occurrence'].corr(storm_yr_sum['Damage'], method='spearman')

print 'Earthquake Spearman correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='spearman')
print 'Storm Spearman correlation from 1980 = ', storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='spearman')
Earthquake Spearman correlation =  0.871554100714
Storm Spearman correlation =  0.912862080289
Earthquake Spearman correlation from 1980 =  0.330575499317
Storm Spearman correlation from 1980 =  0.748376270905

...with rolling mean

In [89]:
print 'Earthquake Spearman correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(eq_yr_sum['Damage'], 3, center=True), method='spearman')
print 'Storm Spearman correlation (rolling mean) = ', rolling_mean(storm_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['Damage'], 3, center=True), method='spearman')

print 'Earthquake Spearman correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='spearman')
print 'Storm Spearman correlation from 1980 (rolling mean) = ', rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='spearman')
Earthquake Spearman correlation (rolling mean) =  0.850376849192
Storm Spearman correlation (rolling mean) =  0.939086822093
Earthquake Spearman correlation from 1980 (rolling mean) =  0.280183656033
Storm Spearman correlation from 1980 (rolling mean) =  0.742506215628

Calcuate linear correlation (Pearson) between occurrence and losses

In [90]:
print 'Earthquake Pearson correlation = ', eq_yr_sum['Occurrence'].corr(eq_yr_sum['Damage'], method='pearson')
print 'Storm Pearson correlation = ', storm_yr_sum['Occurrence'].corr(storm_yr_sum['Damage'], method='pearson')

print 'Earthquake Pearson correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='pearson')
print 'Storm Pearson correlation from 1980 = ', storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='pearson')
Earthquake Pearson correlation =  0.403172662868
Storm Pearson correlation =  0.670543627718
Earthquake Pearson correlation from 1980 =  0.185048723441
Storm Pearson correlation from 1980 =  0.576873804276

...with rolling mean

In [91]:
print 'Earthquake Pearson correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(eq_yr_sum['Damage'], 3, center=True), method='pearson')
print 'Storm Pearson correlation (rolling mean) = ', rolling_mean(storm_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['Damage'], 3, center=True), method='pearson')

print 'Earthquake Pearson correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='pearson')
print 'Storm Pearson correlation from 1980 (rolling mean) = ', rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='pearson')
Earthquake Pearson correlation (rolling mean) =  0.537345456962
Storm Pearson correlation (rolling mean) =  0.772914594633
Earthquake Pearson correlation from 1980 (rolling mean) =  0.0431994288626
Storm Pearson correlation from 1980 (rolling mean) =  0.646522054775

Calculate monotonic correlation coefficient (Spearman) between storm and earthquake indicators

In [92]:
print 'Occurrence Spearman correlation = ', eq_yr_sum['Occurrence'].corr(storm_yr_sum['Occurrence'], method='spearman')
print 'Damage Spearman correlation = ', eq_yr_sum['Damage'].corr(storm_yr_sum['Damage'], method='spearman')

print 'Occurrence Spearman correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], method='spearman')
print 'Damage Spearman correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Damage'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='spearman')
Occurrence Spearman correlation =  0.845976093708
Damage Spearman correlation =  0.758181805369
Occurrence Spearman correlation from 1980 =  0.625568924763
Damage Spearman correlation from 1980 =  0.430099312452

...with rolling mean

In [93]:
print 'Occurrence Spearman correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['Occurrence'], 3, center=True), method='spearman')
print 'Damage Spearman correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Damage'], 3, center=True).corr(rolling_mean(storm_yr_sum['Damage'], 3, center=True), method='spearman')

print 'Occurrence Spearman correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True), method='spearman')
print 'Damage Spearman correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='spearman')
Occurrence Spearman correlation (rolling mean) =  0.916570629247
Damage Spearman correlation (rolling mean) =  0.837228748134
Occurrence Spearman correlation from 1980 (rolling mean) =  0.823162619013
Damage Spearman correlation from 1980 (rolling mean) =  0.600806451613

Calculate linear correlation coefficient (Pearson) between storm and earthquake indicators

In [94]:
print 'Occurrence Pearson correlation = ', eq_yr_sum['Occurrence'].corr(storm_yr_sum['Occurrence'], method='pearson')
print 'Damage Pearson correlation = ', eq_yr_sum['Damage'].corr(storm_yr_sum['Damage'], method='pearson')

print 'Occurrence Pearson correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], method='pearson')
print 'Damage Pearson correlation from 1980 = ', eq_yr_sum['1980-01-01':'2013-01-01']['Damage'].corr(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], method='pearson')
Occurrence Pearson correlation =  0.914548164346
Damage Pearson correlation =  0.369449872073
Occurrence Pearson correlation from 1980 =  0.637817097687
Damage Pearson correlation from 1980 =  0.214650955373

...with rolling mean

In [95]:
print 'Occurrence Pearson correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['Occurrence'], 3, center=True), method='pearson')
print 'Damage Pearson correlation (rolling mean) = ', rolling_mean(eq_yr_sum['Damage'], 3, center=True).corr(rolling_mean(storm_yr_sum['Damage'], 3, center=True), method='pearson')

print 'Occurrence Pearson correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Occurrence'], 3, center=True), method='pearson')
print 'Damage Pearson correlation from 1980 (rolling mean) = ', rolling_mean(eq_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True).corr(rolling_mean(storm_yr_sum['1980-01-01':'2013-01-01']['Damage'], 3, center=True), method='pearson')
Occurrence Pearson correlation (rolling mean) =  0.967032011206
Damage Pearson correlation (rolling mean) =  0.575345078295
Occurrence Pearson correlation from 1980 (rolling mean) =  0.788610735142
Damage Pearson correlation from 1980 (rolling mean) =  0.322886963221

Plot earthquake and storm indicators

In [205]:
fig, ax1 = subplots()


ax1.plot(eq_yr_sum.index, eq_yr_sum['Occurrence'], label='Earthquake Occurrences', color='g', linewidth=3.0)
ax1.plot(storm_yr_sum.index, storm_yr_sum['Occurrence'], label='Storm Occurrences', color='m', linewidth=3.0)

ax1.tick_params(labelsize=20)
ax1.set_xlabel('Year', size='24')
ax1.set_ylabel('Occurrences (EMDAT criteria)', size='24')

frame = legend(loc='upper left', fontsize='20').get_frame()
frame.set_facecolor('1.0')

ax2 = ax1.twinx()
ax2.plot(eq_yr_sum.index, eq_yr_sum['Damage'] / 1000000000.0, label='Earthquake Losses', color='g', linestyle='--', linewidth=3.0)
ax2.plot(storm_yr_sum.index, storm_yr_sum['Damage']  / 1000000000.0, label='Storm Losses', color='m', linestyle='--', linewidth=3.0)
ax2.ticklabel_format(style='plain', axis='y')
ax2.tick_params(labelsize=20)
ax2.set_ylabel('Losses ($billion)', size='24')
#xlim(['1980-01-01', '2015-01-01'])
xlim(['1900-01-01', '2015-01-01'])

frame = legend(loc='upper right', fontsize='20').get_frame()
frame.set_facecolor('1.0')

...with rolling mean

In [142]:
fig, ax1 = subplots()


ax1.plot(eq_yr_sum.index, rolling_mean(eq_yr_sum['Occurrence'], 3), label='Earthquake Occurrences', color='g')
ax1.plot(storm_yr_sum.index, rolling_mean(storm_yr_sum['Occurrence'], 3), label='Storm Occurrences', color='m')

ax1.tick_params(labelsize=20)
ax1.set_xlabel('Year', size='24')
ax1.set_ylabel('Occurrences (EMDAT criteria)', size='24')

frame = legend(loc='upper left', fontsize='20').get_frame()
frame.set_facecolor('1.0')

ax2 = ax1.twinx()
ax2.plot(eq_yr_sum.index, rolling_mean(eq_yr_sum['Damage'] / 1000000000.0, 3), label='Earthquake Losses', color='g', linestyle='--')
ax2.plot(storm_yr_sum.index, rolling_mean(storm_yr_sum['Damage'] / 1000000000.0, 3), label='Storm Losses', color='m', linestyle='--')
ax2.ticklabel_format(style='plain', axis='y')
ax2.tick_params(labelsize=20)
ax2.set_ylabel('Losses ($billion)', size='24')
xlim(['1980-01-01', '2015-01-01'])

frame = legend(loc='upper right', fontsize='20').get_frame()
frame.set_facecolor('1.0')

Difference between storm and earthquake occurrences and losses

In [98]:
storm_eq_occur_diff = (storm_yr_sum['Occurrence'] - eq_yr_sum['Occurrence'] )
storm_eq_dmg_diff = (storm_yr_sum['Damage'] - eq_yr_sum['Damage'] ) 

...with rolling mean

In [99]:
storm_eq_occur_diff = rolling_mean(storm_yr_sum['Occurrence'] - eq_yr_sum['Occurrence'], 3, center=True) 
storm_eq_dmg_diff = rolling_mean(storm_yr_sum['Damage'] - eq_yr_sum['Damage'], 3, center=True ) 

Plot storm/earthquake damage differences

In [210]:
fig, ax1 = subplots()

ax1.plot(storm_eq_occur_diff.index, storm_eq_occur_diff, label='Occurrences', color='g')

ax1.tick_params(labelsize=20)
ax1.set_xlabel('Year', size='24')
ax1.set_ylabel('Storm - Earthquake Occurrences (EMDAT criteria)', size='24')

frame = legend(loc='upper left', fontsize='20').get_frame()
frame.set_facecolor('1.0')

ax2 = ax1.twinx()
ax2.plot(storm_eq_occur_diff.index, storm_eq_dmg_diff / 1000000000.0, label='Damage', color='m')

ax2.ticklabel_format(style='plain', axis='y')
ax2.tick_params(labelsize=20)
ax2.set_ylabel('Storm - Earthquake Losses ($billion)', size='24')

xlim(['1980-01-01', '2015-01-01'])
#xlim(['1900-01-01', '2015-01-01'])

frame = legend(loc='upper right', fontsize='20').get_frame()
frame.set_facecolor('1.0')

...now with rolling mean

In [196]:
fig, ax1 = subplots()

ax1.plot(storm_eq_dmg_diff.index, storm_eq_occur_diff, label='Occurrences', color='g')

ax1.tick_params(labelsize=20)
ax1.set_xlabel('Year', size='24')
ax1.set_ylabel('Storm - Earthquake Occurrences (EMDAT criteria)', size='16')

frame = legend(loc='upper left', fontsize='20').get_frame()
frame.set_facecolor('1.0')

ax2 = ax1.twinx()
ax2.plot(storm_eq_dmg_diff.index, storm_eq_dmg_diff / 1000000000.0, label='Damage ($)', color='m')
ax2.ticklabel_format(style='plain', axis='y')
ax2.tick_params(labelsize=20)
ax2.set_ylabel('Storm - Earthquake Losses ($billion)', size='24')
xlim(['1980-01-01', '2015-01-01'])

frame = legend(loc='upper right', fontsize='20').get_frame()
frame.set_facecolor('1.0')

Differences between storms and earthquakes for other indicators (rolling mean only)

In [102]:
storm_eq_kill_diff = rolling_mean(storm_yr_sum['Killed'] - eq_yr_sum['Killed'], 3, center=True ) 
storm_eq_inj_diff = rolling_mean(storm_yr_sum['Injured'] - eq_yr_sum['Injured'], 3, center=True )
storm_eq_occur_diff = rolling_mean(storm_yr_sum['Occurrence'] - eq_yr_sum['Occurrence'], 3, center=True) 
storm_eq_home_diff = rolling_mean(storm_yr_sum['Homeless'] - eq_yr_sum['Homeless'], 3, center=True )

Determine linear trend lines for differences (on rolling means)

In [172]:
occur_diff_roll = rolling_mean(storm_eq_occur_diff, 3, center=True).dropna()
dmg_diff_roll = rolling_mean(storm_eq_dmg_diff, 3, center=True).dropna()
kill_diff_roll = rolling_mean(storm_eq_kill_diff, 3, center=True).dropna()
inj_diff_roll = rolling_mean(storm_eq_inj_diff, 3, center=True).dropna()
home_diff_roll = rolling_mean(storm_eq_home_diff, 3, center=True).dropna()
years = occur_diff_roll.index.year

occur_slope, occur_int, occur_r, occur_p, occur_std = stats.linregress(years,occur_diff_roll)
occur_trend = occur_slope*years+occur_int

dmg_slope, dmg_int, dmg_r, dmg_p, dmg_std = stats.linregress(years,dmg_diff_roll)
dmg_trend = dmg_slope*years+dmg_int

kill_slope, kill_int, kill_r, kill_p, kill_std = stats.linregress(years,kill_diff_roll)
kill_trend = kill_slope*years+kill_int

inj_slope, inj_int, inj_r, inj_p, inj_std = stats.linregress(years,inj_diff_roll)
inj_trend = inj_slope*years+inj_int

home_slope, home_int, home_r, home_p, home_std = stats.linregress(years,home_diff_roll)
home_trend = home_slope*years+home_int

Now just for 1980 and after

In [174]:
occur_diff_roll1980 = rolling_mean(storm_eq_occur_diff['1980-01-01':'2015-01-01'], 3, center=True).dropna()
dmg_diff_roll1980 = rolling_mean(storm_eq_dmg_diff['1980-01-01':'2015-01-01'], 3, center=True).dropna()
kill_diff_roll1980 = rolling_mean(storm_eq_kill_diff['1980-01-01':'2015-01-01'], 3, center=True).dropna()
inj_diff_roll1980 = rolling_mean(storm_eq_inj_diff['1980-01-01':'2015-01-01'], 3, center=True).dropna()
home_diff_roll1980 = rolling_mean(storm_eq_home_diff['1980-01-01':'2015-01-01'], 3, center=True).dropna()
years1980 = occur_diff_roll1980.index.year

occur_slope1980, occur_int1980, occur_r1980, occur_p1980, occur_std1980 = stats.linregress(years1980,occur_diff_roll1980)
occur_trend1980 = occur_slope1980*years1980+occur_int1980

dmg_slope1980, dmg_int1980, dmg_r1980, dmg_p1980, dmg_std1980 = stats.linregress(years1980,dmg_diff_roll1980)
dmg_trend1980 = dmg_slope1980*years1980+dmg_int1980

kill_slope1980, kill_int1980, kill_r1980, kill_p1980, kill_std1980 = stats.linregress(years1980,kill_diff_roll1980)
kill_trend1980 = kill_slope1980*years1980+kill_int1980

inj_slope1980, inj_int1980, inj_r1980, inj_p1980, inj_std1980 = stats.linregress(years1980,inj_diff_roll1980)
inj_trend1980 = inj_slope1980*years1980+inj_int1980

home_slope1980, home_int1980, home_r1980, home_p1980, home_std1980 = stats.linregress(years1980,home_diff_roll1980)
home_trend1980 = home_slope1980*years1980+home_int1980

Plot differences between storm and earthquake for all indicators (rolling mean)

In [206]:
fig, ax = subplots(5, sharex=True)

ax[0].plot(storm_eq_occur_diff.index, storm_eq_occur_diff, label='Occurrence', color='k')
ax[0].plot(occur_diff_roll.index,occur_trend, linestyle='--', linewidth=1, color='k')
ax[0].set_title('Storm - Earthquake', size='24')
ax[0].set_ylabel('Occurrences')
#ax[0].locator_params(axis='y',nbins = 5)
ax[1].plot(storm_eq_dmg_diff.index, storm_eq_dmg_diff / 1000000000.0, label='Damage', color='k')
ax[1].plot(dmg_diff_roll.index,dmg_trend / 1000000000.0, linestyle='--', linewidth=1, color='k')
ax[1].set_ylabel('Losses \n  ($billion)')
ax[1].locator_params(axis='y',nbins = 5)
ax[2].plot(storm_eq_kill_diff.index, storm_eq_kill_diff / 1000.0, label='Killed', color='k')
ax[2].plot(kill_diff_roll.index,kill_trend / 1000.0, linestyle='--', linewidth=1, color='k')
ax[2].set_ylabel('Deaths \n (1,000s)')
ax[2].locator_params(axis='y',nbins = 5)
ax[3].plot(storm_eq_inj_diff.index, storm_eq_inj_diff / 1000.0, label='Injured', color='k')
ax[3].plot(inj_diff_roll.index,inj_trend / 1000.0, linestyle='--', linewidth=1, color='k')
ax[3].set_ylabel('Injuries \n (1,000s)')
ax[3].locator_params(axis='y',nbins = 5)
ax[4].plot(storm_eq_home_diff.index, storm_eq_home_diff / 1000.0, label='Homeless', color='k')
ax[4].plot(occur_diff_roll.index,home_trend / 1000.0, linestyle='--', linewidth=1, color='k')
ax[4].set_ylabel('Homeless \n (1,000s)')
ax[4].locator_params(axis='y',nbins = 5)
setp(ax[4].get_xticklabels(), fontsize=20)

#xlim(['1980-01-01', '2015-01-01'])
xlim(['1900-01-01', '2015-01-01'])
Out[206]:
(693596.0, 735599.0)