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

Rev	Line
[3231]	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
[3303]	7	from math import sqrt
[3318]	8	from Scientific.Functions.Interpolation import InterpolatingFunction
[3334]	9	from Numeric import array, ravel, Float, zeros
	10	from random import choice
[3303]	11
	12	try:
	13	import kinds
	14	except ImportError:
[3318]	15	# Hand-built mockup of the things we need from the kinds package, since it
[3303]	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()
[3318]	27
[3303]	28
[3334]	29	from utilities.numerical_tools import ensure_numeric
[3318]	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
[3334]	35	depth_epsilon = epsilon
[3318]	36
[3303]	37	def inundation_damage(sww_file, exposure_file_in,
	38	exposure_file_out=None,
[3350]	39	ground_floor_height=0.3,
[3303]	40	overwrite=False, verbose=True,
	41	use_cache = True):
	42	"""
[3346]	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.
[3303]	49	"""
	50
	51	csv = Exposure_csv(exposure_file_in)
	52	geospatial = csv.get_location()
[3318]	53	max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
[3350]	54	geospatial,
	55	ground_floor_height=ground_floor_height,
	56	verbose=verbose,
	57	use_cache=use_cache)
[3337]	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	)
[3355]	64	results_dic = edm.calc_damage_and_costs(verbose_csv=True, verbose=verbose)
[3337]	65	for title, value in results_dic.iteritems():
[3350]	66	csv.set_column(title, value, overwrite=overwrite)
[3337]	67
[3318]	68	# Save info back to csv file
[3346]	69	if exposure_file_out == None:
	70	exposure_file_out = exposure_file_in
[3337]	71	csv.save(exposure_file_out)
[3303]	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()
[3318]	82	max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
[3303]	83	geospatial,
	84	verbose=verbose,
	85	use_cache=use_cache)
[3318]	86	csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
	87	csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
[3303]	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	"""
[3318]	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.
[3303]	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
[3352]	109	max_depths = [-ground_floor_height]*point_count
	110	max_momentums = [-ground_floor_height]*point_count
[3303]	111	for point_i, point in enumerate(points):
	112	for time in callable_sww.get_time():
	113	quantities = callable_sww(time,point_i)
[3310]	114	#print "quantities", quantities
[3309]	115
[3303]	116	w = quantities[0]
	117	z = quantities[1]
	118	uh = quantities[2]
	119	vh = quantities[3]
[3318]	120
[3334]	121	# -ground_floor_height is the minimum value.
[3303]	122	depth = w - z - ground_floor_height
[3309]	123
[3318]	124	if depth > max_depths[point_i]:
	125	max_depths[point_i] = depth
[3309]	126	if w == INF or z == INF or uh == INF or vh == INF:
	127	continue
	128	momentum = sqrt(uhuh + vhvh)
[3318]	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	"""
[3337]	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"
[3350]	143	MAX_DEPTH_TITLE = "Inundation height above ground floor (m)"
[3356]	144	STRUCT_COLLAPSED_TITLE = "collapsed structure if 1"
	145	STRUCT_INUNDATED_TITLE = "inundated structure if 1"
[3334]	146	double_brick_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
	147	[0.0-depth_epsilon, 0.0],
[3337]	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],
[3334]	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],
[3337]	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],
[3334]	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],
[3337]	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],
[3334]	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:]))
[3337]	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]]
[3334]	205
	206	def __init__(self,max_depths, shore_distances, walls,
	207	struct_costs, content_costs):
[3350]	208	"""
	209	max depth is Inundation height above ground floor (m), so
	210	the ground floor has been taken into account.
	211	"""
[3346]	212	self.max_depths = [float(x) for x in max_depths]
	213	self.shore_distances = [float(x) for x in shore_distances]
[3318]	214	self.walls = walls
[3346]	215	self.struct_costs = [float(x) for x in struct_costs]
	216	self.content_costs = [float(x) for x in content_costs]
[3334]	217
[3318]	218	self.structure_count = len(self.max_depths)
[3334]	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
[3355]	226	def calc_damage_and_costs(self, verbose_csv=False, verbose=False):
[3337]	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
[3350]	236	,self.MAX_DEPTH_TITLE:self.max_depths
[3355]	237	,self.STRUCT_COLLAPSED_TITLE:self.struct_collapsed
	238	,self.STRUCT_INUNDATED_TITLE:self.struct_inundated
[3337]	239	}
	240	if verbose_csv:
	241	results_dict[self.COLLAPSE_CSV_INFO_TITLE] = self.collapse_csv_info
	242	return results_dict
	243
[3318]	244	def calc_damage_percentages(self):
	245
	246	# the data being created
[3334]	247	struct_damage = zeros(self.structure_count,Float)
	248	contents_damage = zeros(self.structure_count,Float)
[3355]	249	self.struct_inundated = ['']* self.structure_count
[3318]	250
[3334]	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
[3355]	263
	264	# The definition of inundated is if the max_depth is > 0.0
	265	self.struct_inundated[i] = 1.0
	266
[3334]	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)
[3337]	272
[3334]	273	self.struct_damage = struct_damage
	274	self.contents_damage = contents_damage
[3318]	275
[3303]	276
[3334]	277	def calc_cost(self):
[3346]	278	# ensure_numeric does not cut it.
[3334]	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)
[3337]	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
[3334]	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	"""
[3355]	338
	339	self.struct_collapsed = ['']* self.structure_count
[3334]	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
[3337]	349	self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
[3334]	350	+ str(int(collapse_count)) + ' collapsed out of ' \
[3337]	351	+ str(len(house_indexes)) + ')'
[3334]	352	for _ in range(int(collapse_count)):
	353	house_index = choice(house_indexes)
[3337]	354	self.struct_damage[house_index] = 1.0
	355	self.contents_damage[house_index] = 1.0
[3334]	356	house_indexes.remove(house_index)
[3355]	357	self.struct_collapsed[house_index] = 1
	358
[3334]	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
[3318]	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
[3334]	383	dict = {1:10, 2:20}
	384	print dict.get(100,dict[2])
	385
[3337]	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: