Context Navigation

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

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

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

Download in other formats: