Context Navigation

source: inundation/damage/inundation_damage.py @ 3380

Last change on this file since 3380 was 3356, checked in by duncan, 19 years ago
slight change
File size: 17.0 KB

Line
1	"""Classes for implementing damage curves and economic damage
2
3	Duncan Gray, Ole Nielsen, Jane Sexton, Nick Bartzis
4	Geoscience Australia, 2006
5	"""
6
7	from math import sqrt
8	from Scientific.Functions.Interpolation import InterpolatingFunction
9	from Numeric import array, ravel, Float, zeros
10	from random import choice
11
12	try:
13	import kinds
14	except ImportError:
15	# Hand-built mockup of the things we need from the kinds package, since it
16	# was recently removed from the standard Numeric distro. Some users may
17	# not have it by default.
18	class _bunch:
19	pass
20
21	class _kinds(_bunch):
22	default_float_kind = _bunch()
23	default_float_kind.MIN = 2.2250738585072014e-308 #smallest +ve number
24	default_float_kind.MAX = 1.7976931348623157e+308
25
26	kinds = _kinds()
27
28
29	from utilities.numerical_tools import ensure_numeric
30	from pyvolution.data_manager import Exposure_csv
31	from pyvolution.util import file_function
32	from geospatial_data.geospatial_data import ensure_absolute
33	from utilities.numerical_tools import INF
34	from anuga_config import epsilon
35	depth_epsilon = epsilon
36
37	def inundation_damage(sww_file, exposure_file_in,
38	exposure_file_out=None,
39	ground_floor_height=0.3,
40	overwrite=False, verbose=True,
41	use_cache = True):
42	"""
43	This is the main function for calculating tsunami damage due to
44	inundation. It gets the location of structures from the exposure
45	file and gets the inundation of these structures from the
46	sww file.
47
48	It then calculates the damage loss.
49	"""
50
51	csv = Exposure_csv(exposure_file_in)
52	geospatial = csv.get_location()
53	max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
54	geospatial,
55	ground_floor_height=ground_floor_height,
56	verbose=verbose,
57	use_cache=use_cache)
58	edm = EventDamageModel(max_depths,
59	csv.get_column('SHORE_DIST'),
60	csv.get_column('WALLS'),
61	csv.get_column('STR_VALUE'),
62	csv.get_column('C_VALUE')
63	)
64	results_dic = edm.calc_damage_and_costs(verbose_csv=True, verbose=verbose)
65	for title, value in results_dic.iteritems():
66	csv.set_column(title, value, overwrite=overwrite)
67
68	# Save info back to csv file
69	if exposure_file_out == None:
70	exposure_file_out = exposure_file_in
71	csv.save(exposure_file_out)
72
73	def add_depth_and_momentum2csv(sww_file, exposure_file_in,
74	exposure_file_out=None,
75	overwrite=False, verbose=True,
76	use_cache = True):
77	"""
78	"""
79
80	csv = Exposure_csv(exposure_file_in)
81	geospatial = csv.get_location()
82	max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
83	geospatial,
84	verbose=verbose,
85	use_cache=use_cache)
86	csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
87	csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
88	csv.save(exposure_file_out)
89
90	def calc_max_depth_and_momentum(sww_file, points,
91	ground_floor_height=0.0,
92	verbose=True,
93	use_cache = True):
94	"""
95	Calculate the maximum inundation height above ground floor (mihagf) for a list
96	of locations.
97
98	The mihagf value is in the range -ground_floor_height to overflow errors.
99	"""
100	quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
101	points = ensure_absolute(points)
102	point_count = len(points)
103	callable_sww = file_function(sww_file,
104	quantities=quantities,
105	interpolation_points=points,
106	verbose=verbose,
107	use_cache=use_cache)
108	# initialise the max lists
109	max_depths = [-ground_floor_height]*point_count
110	max_momentums = [-ground_floor_height]*point_count
111	for point_i, point in enumerate(points):
112	for time in callable_sww.get_time():
113	quantities = callable_sww(time,point_i)
114	#print "quantities", quantities
115
116	w = quantities[0]
117	z = quantities[1]
118	uh = quantities[2]
119	vh = quantities[3]
120
121	# -ground_floor_height is the minimum value.
122	depth = w - z - ground_floor_height
123
124	if depth > max_depths[point_i]:
125	max_depths[point_i] = depth
126	if w == INF or z == INF or uh == INF or vh == INF:
127	continue
128	momentum = sqrt(uhuh + vhvh)
129	if momentum > max_momentums[point_i]:
130	max_momentums[point_i] = momentum
131	return max_depths, max_momentums
132
133	class EventDamageModel:
134	"""
135
136
137	"""
138	STRUCT_LOSS_TITLE = "Structure Loss ($)"
139	CONTENTS_LOSS_TITLE = "Contents Loss ($)"
140	CONTENTS_DAMAGE_TITLE = "Contents damaged (fraction)"
141	STRUCT_DAMAGE_TITLE = "Structure damaged (fraction)"
142	COLLAPSE_CSV_INFO_TITLE = "Calculation notes"
143	MAX_DEPTH_TITLE = "Inundation height above ground floor (m)"
144	STRUCT_COLLAPSED_TITLE = "collapsed structure if 1"
145	STRUCT_INUNDATED_TITLE = "inundated structure if 1"
146	double_brick_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
147	[0.0-depth_epsilon, 0.0],
148	[0.0,0.016],
149	[0.1,0.150],
150	[0.3,0.425],
151	[0.5,0.449],
152	[1.0,0.572],
153	[1.5,0.582],
154	[2.0,0.587],
155	[2.5,0.647],
156	[kinds.default_float_kind.MAX,64.7]])
157	double_brick_damage_curve = InterpolatingFunction( \
158	(ravel(double_brick_damage_array[:,0:1]),),
159	ravel(double_brick_damage_array[:,1:]))
160
161	brick_veeer_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
162	[0.0-depth_epsilon, 0.0],
163	[0.0,0.016],
164	[0.1,0.169],
165	[0.3,0.445],
166	[0.5,0.472],
167	[1.0,0.618],
168	[1.5,0.629],
169	[2.0,0.633],
170	[2.5,0.694],
171	[kinds.default_float_kind.MAX,69.4]])
172	brick_veeer_damage_curve = InterpolatingFunction( \
173	(ravel(brick_veeer_damage_array[:,0:1]),),
174	ravel(brick_veeer_damage_array[:,1:]))
175	struct_damage_curve = {'Double Brick':double_brick_damage_curve,
176	'Brick Veneer':brick_veeer_damage_curve}
177	default_struct_damage_curve = brick_veeer_damage_curve
178
179	contents_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
180	[0.0-depth_epsilon, 0.0],
181	[0.0,0.013],
182	[0.1,0.102],
183	[0.3,0.381],
184	[0.5,0.500],
185	[1.0,0.970],
186	[1.5,0.976],
187	[2.0,0.986],
188	[kinds.default_float_kind.MAX,98.6]])
189	contents_damage_curve = InterpolatingFunction( \
190	(ravel(contents_damage_array[:,0:1]),),
191	ravel(contents_damage_array[:,1:]))
192
193	#building collapse probability
194	# inundation depth above ground floor, m
195	depth_upper_limits = [depth_epsilon, 1.0, 2.0, 3.0, 5.0, kinds.default_float_kind.MAX]
196	# shore mistance, m
197	shore_upper_limits = [125,200,250, kinds.default_float_kind.MAX]
198	# Building collapse probability
199	collapse_probability = [[0.0, 0.0, 0.0, 0.0], #Code below assumes 0.0
200	[0.05, 0.02, 0.01, 0.0],
201	[0.6, 0.3, 0.1, 0.05],
202	[0.8, 0.4, 0.25, 0.15],
203	[0.95, 0.7, 0.5, 0.3],
204	[0.99, 0.9, 0.65, 0.45]]
205
206	def __init__(self,max_depths, shore_distances, walls,
207	struct_costs, content_costs):
208	"""
209	max depth is Inundation height above ground floor (m), so
210	the ground floor has been taken into account.
211	"""
212	self.max_depths = [float(x) for x in max_depths]
213	self.shore_distances = [float(x) for x in shore_distances]
214	self.walls = walls
215	self.struct_costs = [float(x) for x in struct_costs]
216	self.content_costs = [float(x) for x in content_costs]
217
218	self.structure_count = len(self.max_depths)
219	#Fixme expand
220	assert self.structure_count == len(self.shore_distances)
221	assert self.structure_count == len(self.walls)
222	assert self.structure_count == len(self.struct_costs)
223	assert self.structure_count == len(self.content_costs)
224	#assert self.structure_count == len(self.)
225
226	def calc_damage_and_costs(self, verbose_csv=False, verbose=False):
227	self.calc_damage_percentages()
228	collapse_probability = self.calc_collapse_probability()
229	self._calc_collapse_structures(collapse_probability,
230	verbose_csv=verbose_csv)
231	self.calc_cost()
232	results_dict = {self.STRUCT_LOSS_TITLE:self.struct_loss
233	,self.STRUCT_DAMAGE_TITLE:self.struct_damage
234	,self.CONTENTS_LOSS_TITLE:self.contents_loss
235	,self.CONTENTS_DAMAGE_TITLE:self.contents_damage
236	,self.MAX_DEPTH_TITLE:self.max_depths
237	,self.STRUCT_COLLAPSED_TITLE:self.struct_collapsed
238	,self.STRUCT_INUNDATED_TITLE:self.struct_inundated
239	}
240	if verbose_csv:
241	results_dict[self.COLLAPSE_CSV_INFO_TITLE] = self.collapse_csv_info
242	return results_dict
243
244	def calc_damage_percentages(self):
245
246	# the data being created
247	struct_damage = zeros(self.structure_count,Float)
248	contents_damage = zeros(self.structure_count,Float)
249	self.struct_inundated = ['']* self.structure_count
250
251	for i,max_depth,shore_distance,wall in map(None,
252	range(self.structure_count),
253	self.max_depths,
254	self.shore_distances,
255	self.walls):
256	#print "i",i
257	#print "max_depth",max_depth
258	#print "shore_distance",shore_distance
259	#print "wall",wall
260	## WARNING SKIP IF DEPTH < 0.0
261	if 0.0 > max_depth:
262	continue
263
264	# The definition of inundated is if the max_depth is > 0.0
265	self.struct_inundated[i] = 1.0
266
267	#calc structural damage %
268	damage_curve = self.struct_damage_curve.get(wall,
269	self.default_struct_damage_curve)
270	struct_damage[i] = damage_curve(max_depth)
271	contents_damage[i] = self.contents_damage_curve(max_depth)
272
273	self.struct_damage = struct_damage
274	self.contents_damage = contents_damage
275
276
277	def calc_cost(self):
278	# ensure_numeric does not cut it.
279	self.struct_loss = self.struct_damage * \
280	ensure_numeric(self.struct_costs)
281	self.contents_loss = self.contents_damage * \
282	ensure_numeric(self.content_costs)
283
284	def calc_collapse_probability(self):
285	"""
286	return a dict of which structures have x probability of collapse.
287	key is collapse probability
288	value is list of struct indexes with key probability of collapse
289	"""
290	# I could've done this is the calc_damage_percentages and
291	# Just had one loop.
292	# But for ease of testing and bug finding I'm seperating the loops.
293	# I'm make the outer loop for both of them the same though,
294	# so this loop can easily be folded into the other loop.
295
296	# dict of which structures have x probability of collapse.
297	# key of collapse probability
298	# value of list of struct indexes
299	struct_coll_prob = {}
300
301	for i,max_depth,shore_distance,wall in map(None,
302	range(self.structure_count),
303	self.max_depths,
304	self.shore_distances,
305	self.walls):
306	#print "i",i
307	#print "max_depth",max_depth
308	#print "shore_distance",shore_distance
309	#print "wall",wall
310	# WARNING ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
311	depth_upper_limits = self.depth_upper_limits
312	shore_upper_limits = self.shore_upper_limits
313	collapse_probability = self.collapse_probability
314	if max_depth <= depth_upper_limits[0]:
315	continue
316	start = 1
317	for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
318	#Have to change i_depth so it indexes into the lists correctly
319	i_depth += start
320	if max_depth <= depth_limit:
321	for i_shore, shore_limit in enumerate(shore_upper_limits):
322	if shore_distance <= shore_limit:
323	coll_prob = collapse_probability[i_depth][i_shore]
324	if 0.0 == collapse_probability[i_depth][i_shore]:
325	break
326	struct_coll_prob.setdefault(coll_prob,[]).append(i)
327	break
328	break
329
330	return struct_coll_prob
331
332	def _calc_collapse_structures(self, collapse_probability,
333	verbose_csv=False):
334	"""
335	Given the collapse probabilities, throw the dice
336	and collapse some houses
337	"""
338
339	self.struct_collapsed = ['']* self.structure_count
340	if verbose_csv:
341	self.collapse_csv_info = ['']* self.structure_count
342	#for a given 'bin', work out how many houses will collapse
343	for probability, house_indexes in collapse_probability.iteritems():
344	collapse_count = round(len(house_indexes) *probability)
345
346	if verbose_csv:
347	for i in house_indexes:
348	# This could be sped up I think
349	self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
350	+ str(int(collapse_count)) + ' collapsed out of ' \
351	+ str(len(house_indexes)) + ')'
352	for _ in range(int(collapse_count)):
353	house_index = choice(house_indexes)
354	self.struct_damage[house_index] = 1.0
355	self.contents_damage[house_index] = 1.0
356	house_indexes.remove(house_index)
357	self.struct_collapsed[house_index] = 1
358
359	# Warning, the collapse_probability list now lists
360	# houses that did not collapse, (though not all of them)
361	#print "",self.collapse_csv_info
362	#############################################################################
363	if __name__ == "__main__":
364	from Scientific.Functions.Interpolation import InterpolatingFunction
365	from Numeric import array, ravel
366	a = array([[0,0],[1,10]])
367	c = array([0,1])
368
369	axis = ravel(a[:,0:1])
370	values = ravel(a[:,1:])
371
372	#axis = array([0,1])
373	#values = array([0,10])
374	print "axis",axis
375	print "values",values
376	i = InterpolatingFunction((axis,), values)
377	print "value",i(0.5)
378
379
380	i = InterpolatingFunction((array([0,1]),),array([0,10]))
381	print "value",i(0.5)
382
383	dict = {1:10, 2:20}
384	print dict.get(100,dict[2])
385
386	inundation_damage('source.sww', 'augmented_buildings_high.csv',
387	'augmented_buildings_high_dsg.csv')

Note: See TracBrowser for help on using the repository browser.

Download in other formats: