Context Navigation

source: anuga_core/source/anuga/damage_modelling/inundation_damage.py @ 4605

Last change on this file since 4605 was 4585, checked in by nick, 18 years ago
add print statement
File size: 18.6 KB

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

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

Download in other formats: