Context Navigation

source: trunk/anuga_core/anuga/utilities/quantity_setting_functions.py @ 9679

Last change on this file since 9679 was 9679, checked in by davies, 10 years ago
Bugfix in new nan_interpolation_region functionality of composite_quantity_setting_functions
File size: 26.1 KB

Line
1	"""
2
3	Function which can be useful when setting quantities
4
5	"""
6	import copy
7	import os
8	import anuga.utilities.spatialInputUtil as su
9
10	def make_nearestNeighbour_quantity_function(
11	quantity_xyValueIn,
12	domain,
13	threshold_distance = 9.0e+100,
14	background_value = 9.0e+100,
15	k_nearest_neighbours = 1,
16	):
17	"""
18	Function which makes another function, which can be used in set_quantity
19
20	Idea: For every point x,y in the domain, we want to set a quantity based on
21	the 'nearest-neighbours' from quantity_xyValue (a 3 column array with
22	x,y,quantity-value),
23	UNLESS the distance from x,y to the nearest-neighbour is >
24	threshold_distance.
25	In the latter case, we want to set the quantity value to
26	'background_value'
27
28	We need a function f(x,y) to do that. This routine makes the
29	function, with the desired quantity_xyValue points,
30	threshold_distance, and background_value
31
32	INPUTS:
33	@param quantity_xyValueIn -- A 3 column array with 'x,y, Value'
34	defining the points used to set the new quantity values in
35	georeferenced coordinates
36	@param domain -- The ANUGA domain
37	@param k_nearest_neighbors If > 1, then an inverse-distance-weighted average
38	of the k nearest neighbours is used
39	@param threshold_distance -- Points greater than this distance from
40	their nearest quantity_xyValue point are set to background_value
41	@param background_value -- see 'threshold_distance'
42
43	OUTPUTS:
44	A function f which can be passed to domain.set_quantity('myQuantity', f)
45
46	"""
47
48	import scipy
49	import scipy.interpolate
50	import scipy.spatial
51
52	if(len(quantity_xyValueIn.shape) > 1):
53	quantity_xyValue = quantity_xyValueIn
54	else:
55	# Treat the single-point case
56	quantity_xyValue = quantity_xyValueIn.reshape((1,3))
57
58	# Make a function which gives us the ROW-INDEX of the nearest xy point in
59	# quantity_xyValue
60	#quantity_xy_interpolator = scipy.interpolate.NearestNDInterpolator(
61	# quantity_xyValue[:,0:2],
62	# scipy.arange(len(quantity_xyValue[:,2])))
63
64	# Make a function which returns k-nearest-neighbour indices + distances
65	quantity_xy_interpolator = scipy.spatial.cKDTree(quantity_xyValue[:,0:2])
66
67	# Make a function of x,y which we can pass to domain.set_quantity
68	def quant_NN_fun(x,y):
69	"""
70	Function to assign quantity from the nearest point in quantity_xyValue,
71	UNLESS the point is more than 'threshold_distance' away from the
72	nearest point, in which case the background value is used
73
74	"""
75
76	import scipy
77	import scipy.interpolate
78	import scipy.spatial
79
80	# Since ANUGA stores x,y internally in non-georeferenced coordinates,
81	# we adjust them here
82	xll = domain.geo_reference.xllcorner
83	yll = domain.geo_reference.yllcorner
84	z = scipy.zeros(shape=(len(x), 2))
85	z[:,0] = x+xll
86	z[:,1] = y+yll
87
88	# This will hold the quantity values
89	quantity_output = x*0. + background_value
90	# Compute the index of the nearest-neighbour in quantity_xyValue
91	neighbour_data = quantity_xy_interpolator.query(z,
92	k=k_nearest_neighbours)
93
94	# Next find indices with distance < threshold_distance
95	if(k_nearest_neighbours==1):
96	dist_lt_thresh = neighbour_data[0] < threshold_distance
97	else:
98	dist_lt_thresh = neighbour_data[0][:,0] < threshold_distance
99
100	dist_lt_thresh = dist_lt_thresh.nonzero()[0]
101
102	# Initialise output
103	quantity_output = x*0 + background_value
104
105	# Interpolate
106	if len(dist_lt_thresh)>0:
107	numerator = 0
108	denominator = 0
109	for i in range(k_nearest_neighbours):
110	if(k_nearest_neighbours==1):
111	distances = neighbour_data[0][dist_lt_thresh]
112	indices = neighbour_data[1][dist_lt_thresh]
113	else:
114	distances = neighbour_data[0][dist_lt_thresh,i]
115	indices = neighbour_data[1][dist_lt_thresh,i]
116
117	inverse_distance = 1.0/(distances+1.0e-100)
118	values = quantity_xyValue[indices,2]
119	numerator += values*inverse_distance
120	denominator += inverse_distance
121
122	quantity_output[dist_lt_thresh] = numerator/denominator
123
124	return quantity_output
125
126	# Return the quantity function
127	return quant_NN_fun
128
129	###############################################################################
130
131	def composite_quantity_setting_function(poly_fun_pairs,
132	domain,
133	clip_range = None,
134	nan_treatment = 'exception',
135	nan_interpolation_region_polygon = None,
136	default_k_nearest_neighbours = 1,
137	default_raster_interpolation = 'pixel',
138	verbose=True):
139	""" Make a 'composite function' to set quantities -- applies different
140	functions inside different polygon regions.
141
142	poly_fun_pairs = [ [p0, f0], [p1, f1], ...]
143
144	Where:
145
146	fi is a function,
147	or a constant,
148	or a '.txt' or '.csv' file with comma separated xyz data
149	and an optional header row which contains letters,
150	or the name of a gdal-compatible rasterFile
151	(not ending in .txt or .csv),
152	or a numpy array with 3 columns
153
154	pi is a polygon (anuga polygon format),
155	or a polygon filename (shapefile or a csv format that
156	anuga.read_polygon will read),
157	or None ( equivalent to a polygon with zero area),
158	or 'All' (equivalent to a polygon covering everything)
159	or 'Extent' in the case that fi is a rasterFile name
160	(equivalent to a polygon with the same extent as the raster)
161
162	IMPORTANT: When polygons overlap, the first elements of the list are
163	given priority. The approach is:
164	First f0 is applied to all points in p0, and we record
165	that these points have been 'set'
166	Next f1 is applied to all points in p1 which have not
167	been 'set', and then we record those points as being 'set'
168	Next f2 is applied to all points in p2 which have not
169	been 'set', and then we record those points as being 'set'
170	... etc
171
172	INPUT:
173	@param poly_fun_pairs = [ [p0, f0], [p1, f1], ...]
174
175	where fi(x,y) is a function returning quantity values at points,
176	or any of the special cases below
177
178	SPECIAL fi CASES:
179	fi = a constant in which case points in the polygon are
180	set to that value,
181	fi = a .txt or .csv file name containing x, y, z data,
182	with comma separators and an optional header row
183	containing letters (nearest neighbour interpolation is used)
184	fi = a string rasterFile name (not ending in .txt or .csv)
185	which can be passed to quantityRasterFun to make a function
186	fi = a numpy array with 3 columns (x,y,Value) in which case
187	nearest-neighbour interpolation is used on the points
188
189	pi are polygons where we want to use fi inside
190	(anuga polygon format) or any of the special cases below
191	SPECIAL pi CASES:
192	If pi is a filename ending in .shp or a csv format that
193	anuga.read_polygon can read, we assume it contains a polygon
194	we have to read
195	If any pi = 'All', then we assume that ALL unset points are set
196	using the function. This CAN ONLY happen in the last [fi,pi]
197	pair where pi is not None (since fi will be applied to
198	all remaining points -- so anything else is probably an
199	input mistake)
200	If any pi = None, then that [fi,pi] pair is skipped
201	If pi = 'Extent' and fi is the name of a raster file, then the
202	extent of the raster file is used to define the polygon
203
204	@param domain = ANUGA domain object
205
206	@param clip_range = List with the same length as poly_fun_pairs,
207	of the form:
208	[ [min0, max0], [min1, max1], ...]
209	After f0 is applied in p0, its values will be 'clipped' to the
210	range
211	[min0, max0]
212	, and similarly for the other fi
213
214	@param nan_treatment = 'exception' or 'fall_through' -- string determining
215	what to do if F(x,y) is nan. The default 'exception' raises an exception.
216	The value 'fall_through' allows the function to try lower-priority
217	poly,fun pairs (in sequence) to set the value.
218
219	@param nan_interpolation_region_polygon = None, or 'All', or a list
220	of csv or shp filenames containing polygons, or a list of
221	anuga polygon objects.
222
223	If it is not None, then all x,y points which evaluate to nan
224	on their first preference dataset are recorded, and as a
225	final step, the values at these x,y points
226	which are inside the nan_interpolation_region_polygon
227	are interpolated from the other x,y,F(x,y) values.
228
229	Nearest neighbour interpolation is used, with
230	k_nearest_neighbours taken from default_k_nearest_neighbours.
231
232	Note that if nan_treatment = 'exception', then nan's will cause
233	exceptions earlier on in this routine, so you will need
234	nan_treatment = 'fall_through' to use this option.
235
236	Example of why you might want this:
237	Say you have 2 elevation datasets (one defining the
238	topography above MSL, and the other defining the topography
239	below MSL). There might be small nan gaps between them,
240	which you would like to fill with interpolation. That
241	can be done with this option, by including the nan regions
242	in one of the elevation-dataset-polygons pi.
243
244	@param default_k_nearest_neighbours = integer >=1 . The value of
245	k_nearest_neighbours passed to
246	make_nearestNeighbour_quantity_function when a 'special_case'
247	value of fi is passed in (either a point array or a .txt or
248	.csv point file), or when nan_interpolation_region_polygon is
249	not None
250
251	@param default_raster_interpolation = 'pixel' or 'bilinear'. The value of
252	'interpolation' passed to quantityRasterFun if a raster filename
253	is passed as one of the fi.
254
255	@param verbose TRUE/FALSE Print more information
256
257	OUTPUT: A function F(x,y) which can be used e.g. to set the quantity
258	domain.set_quantity('elevation', F)
259
260	"""
261	import os
262	import numpy
263	from anuga.geometry.polygon import inside_polygon
264
265
266	# Check that clip_range has the right form
267	if clip_range is not None:
268	if len(clip_range) != len(poly_fun_pairs):
269	msg = ' clip_range must be the same ' +\
270	'length as poly_fun_pairs, or None'
271	raise ValueError(msg)
272	# Check that min < = max
273	for i in range(len(clip_range)):
274	if clip_range[i][0] > clip_range[i][1]:
275	raise Exception('clip_range minima must be less than maxima')
276
277
278	def F(x,y):
279	"""This is the function returned by composite_quantity_setting_function
280	It can be passed to set_quantity
281	"""
282	isSet = numpy.zeros(len(x)) # 0/1 - record if each point has been set
283	quantityVal = x*0 + numpy.nan # Function return value
284
285	# Record points which evaluated to nan on their first preference
286	# dataset.
287	was_ever_nan = (x*0).astype(int)
288
289	lpf = len(poly_fun_pairs)
290	if(lpf <= 0):
291	raise Exception('Must have at least 1 fun-poly-pair')
292
293	# Make an array of 'transformed' spatial coordinates, for checking
294	# polygon inclusion
295	xll = domain.geo_reference.xllcorner
296	yll = domain.geo_reference.yllcorner
297	xy_array_trans = numpy.vstack([x+xll,y+yll]).transpose()
298
299	# Check that none of the pi polygons [except perhaps the last] is 'All'
300	for i in range(lpf-1):
301	if(poly_fun_pairs[i][0]=='All'):
302	# This is only ok if all the othe poly_fun_pairs are None
303	remaining_poly_fun_pairs_are_None = \
304	[poly_fun_pairs[j][0] is None for j in range(i+1,lpf)]
305	if(not all(remaining_poly_fun_pairs_are_None)):
306	raise Exception('Can only have the last polygon = All')
307
308	# Main Loop
309	# Apply the fi inside the pi
310	for i in range(lpf):
311	fi = poly_fun_pairs[i][1] # The function
312	pi = poly_fun_pairs[i][0] # The polygon
313
314	# Quick exit
315	if(pi is None):
316	continue
317
318	###################################################################
319	# Get indices fInds of points in polygon pi which are not already
320	# set
321	###################################################################
322	if(pi == 'All'):
323	# Get all unset points
324	fInside = (1-isSet)
325	fInds = (fInside==1).nonzero()[0]
326
327	else:
328
329	if(pi == 'Extent'):
330	# Here fi MUST be a gdal-compatible raster
331	if(not (type(fi) == str)):
332	msg = ' pi = "Extent" can only be used when fi is a' +\
333	' raster file name'
334	raise Exception(msg)
335
336	if(not os.path.exists(fi)):
337	msg = 'fi ' + str(fi) + ' is supposed to be a ' +\
338	' raster filename, but it could not be found'
339	raise Exception(msg)
340
341	# Then we get the extent from the raster itself
342	pi_path = su.getRasterExtent(fi,asPolygon=True)
343
344	elif( (type(pi) == str) and os.path.isfile(pi) ):
345	# pi is a file
346	pi_path = su.read_polygon(pi)
347
348	else:
349	# pi is the actual polygon data
350	pi_path = pi
351
352	# Get the insides of unset points inside pi_path
353	notSet = (isSet==0.).nonzero()[0]
354	fInds = inside_polygon(xy_array_trans[notSet,:], pi_path)
355	fInds = notSet[fInds]
356
357	if len(fInds) == 0:
358	# No points found, move on
359	continue
360
361	###################################################################
362	# Evaluate fi at the points inside pi
363	###################################################################
364
365	# We use various tricks to infer whether fi is a function,
366	# a constant, a file (raster or csv), or an array
367	if(hasattr(fi,'__call__')):
368	# fi is a function
369	quantityVal[fInds] = fi(x[fInds], y[fInds])
370
371	elif isinstance(fi, (int, long, float)):
372	# fi is a numerical constant
373	quantityVal[fInds] = fi*1.0
374
375	elif ( type(fi) is str and os.path.exists(fi)):
376	# fi is a file which is assumed to be
377	# a gdal-compatible raster OR an x,y,z elevation file
378	if os.path.splitext(fi)[1] in ['.txt', '.csv']:
379	fi_array = su.read_csv_optional_header(fi)
380	# Check the results
381	if fi_array.shape[1] is not 3:
382	print 'Treated input file ' + fi +\
383	' as xyz array with an optional header'
384	msg = 'Array should have 3 columns -- x,y,value'
385	raise Exception(msg)
386
387	newfi = make_nearestNeighbour_quantity_function(
388	fi_array, domain,
389	k_nearest_neighbours = default_k_nearest_neighbours)
390	quantityVal[fInds] = newfi(x[fInds], y[fInds])
391
392	else:
393	# Treating input file as a raster
394	newfi = quantityRasterFun(domain, fi,
395	interpolation = default_raster_interpolation)
396	quantityVal[fInds] = newfi(x[fInds], y[fInds])
397
398	elif(type(fi) is numpy.ndarray):
399	if fi.shape[1] is not 3:
400	msg = 'Array should have 3 columns -- x,y,value'
401	raise Exception(msg)
402	newfi = make_nearestNeighbour_quantity_function(fi, domain,
403	k_nearest_neighbours = default_k_nearest_neighbours)
404	quantityVal[fInds] = newfi(x[fInds], y[fInds])
405
406	else:
407	print 'Error with function from'
408	print fi
409	msg='Cannot make function from type ' + str(type(fi))
410	raise Exception, msg
411
412	###################################################################
413	# Check for nan values
414	###################################################################
415
416	nan_flag = (quantityVal[fInds] != quantityVal[fInds])
417	nan_inds = nan_flag.nonzero()[0]
418	was_ever_nan[fInds[nan_inds]] = 1
419
420	if len(nan_inds)>0:
421	if nan_treatment == 'exception':
422	msg = 'nan values generated by the poly_fun_pair at '\
423	'index ' + str(i) + ' '\
424	'in composite_quantity_setting_function. ' + \
425	'To allow these values to be set by later ' + \
426	'poly_fun pairs, pass the argument ' + \
427	'nan_treatment="fall_through" ' + \
428	'to composite_quantity_setting_function'
429	raise Exception(msg)
430
431	elif nan_treatment == 'fall_through':
432	msg = 'WARNING: nan values generated by the ' + \
433	'poly_fun_pair at index ' + str(i) + ' '\
434	'in composite_quantity_setting_function. ' + \
435	'They will be passed to later poly_fun_pairs'
436	if verbose: print msg
437	not_nan_inds = (1-nan_flag).nonzero()[0]
438
439	if len(not_nan_inds)>0:
440	fInds = fInds[not_nan_inds]
441	else:
442	# All values are nan
443	msg = '( Actually all the values were nan - ' + \
444	'Are you sure they should be? Possible error?)'
445	if verbose: print msg
446	continue
447
448	else:
449	msg = 'Found nan values in ' + \
450	'composite_quantity_setting_function but ' + \
451	'nan_treatment is not a recognized value'
452	raise Exception(msg)
453
454	# Record that the points have been set
455	isSet[fInds] = 1
456
457	# Enforce clip_range
458	if clip_range is not None:
459	lower_bound = clip_range[i][0]
460	upper_bound = clip_range[i][1]
461	quantityVal[fInds] = numpy.maximum(
462	quantityVal[fInds], lower_bound)
463	quantityVal[fInds] = numpy.minimum(
464	quantityVal[fInds], upper_bound)
465
466	# End of loop
467
468	# Find points which were nan on their first preference dataset + are
469	# inside nan_interpolation_region_polygon. Then reinterpolate their
470	# values from the other x,y, quantityVal points.
471	if (nan_interpolation_region_polygon is not None) &\
472	(was_ever_nan.sum() > 0):
473	if nan_interpolation_region_polygon == 'All':
474	points_to_reinterpolate = was_ever_nan.nonzero()[0]
475	else:
476	# nan_interpolation_region_polygon contains information on 1 or
477	# more polygons
478	# Inside those polygons, we need to re-interpolate points which
479	# first evaluted to na
480	possible_points_to_reint = was_ever_nan.nonzero()[0]
481	points_to_reinterpolate = numpy.array([]).astype(int)
482
483	for i in range(len(nan_interpolation_region_polygon)):
484	nan_pi = nan_interpolation_region_polygon[i]
485
486	# Ensure nan_pi = list of x,y points making a polygon
487	if(type(nan_pi) == str):
488	nan_pi = su.read_polygon(nan_pi)
489
490	points_in_nan_pi = inside_polygon(
491	xy_array_trans[possible_points_to_reint,:],
492	nan_pi)
493
494	if len(points_in_nan_pi)>0:
495	points_to_reinterpolate = numpy.hstack(
496	[points_to_reinterpolate,
497	possible_points_to_reint[points_in_nan_pi]])
498
499
500	if verbose:
501	print 'Re-interpolating ', len(points_to_reinterpolate),\
502	' points which were nan under their',\
503	' first-preference and are inside the',\
504	' nan_interpolation_region_polygon'
505
506	# Find the interpolation points = points not needing reinterpolation
507	ip = x*0 + 1
508	ip[points_to_reinterpolate] = 0
509	number_of_ip = ip.sum()
510	ip = ip.nonzero()[0]
511
512	if(number_of_ip < default_k_nearest_neighbours):
513	raise Exception('Too few non-nan points to interpolate from')
514
515	# Make function for re-interpolation. Note this requires
516	# x,y,z in georeferenced coordinates, whereas x,y are ANUGA
517	# coordinates
518	reinterp_F = make_nearestNeighbour_quantity_function(
519	numpy.vstack([xy_array_trans[ip,0], xy_array_trans[ip,1],
520	quantityVal[ip]]).transpose(),
521	domain,
522	k_nearest_neighbours = default_k_nearest_neighbours)
523
524	# re-interpolate
525	quantityVal[points_to_reinterpolate] = reinterp_F(
526	x[points_to_reinterpolate], y[points_to_reinterpolate])
527
528	isSet[points_to_reinterpolate] = 1
529
530	# Check there are no remaining nan values
531	if( min(isSet) != 1):
532	print 'Some points remain as nan, which is not allowed'
533	unset_inds = (isSet!=1).nonzero()[0]
534	lui = min(5, len(unset_inds))
535	print 'There are ', len(unset_inds), ' such points'
536	print 'Here are a few:'
537	for i in range(lui):
538	print x[unset_inds[i]]+xll, y[unset_inds[i]]+yll
539	raise Exception('It seems the input data needs to be fixed')
540
541	return quantityVal
542	# END OF FUNCTION F(x,y)
543
544	return F
545
546	##############################################################################
547
548	def quantityRasterFun(domain, rasterFile, interpolation='pixel'):
549	"""
550	Make a function whick takes x,y in ANUGA coordinates, and returns the values
551	on a raster rasterFile
552
553	This can be used to set a quantity, and takes care of the manual conversion
554	from ANUGA coordinates to spatial coordinates.
555
556	INPUTS: @param domain = ANUGA domain
557	@param rasterFile = Filename of the raster to extract point values
558	from
559	@param interpolation = 'pixel' (in which case the point value is
560	set from the pixel it is on) or 'bilinear' in which case
561	the point value is set from bilinear interpolation of
562	pixels.
563
564	OUTPUT: Function which takes x,y in ANUGA coordinates, and outputs their
565	corresponding raster values
566	"""
567	import scipy
568	from anuga.utilities.spatialInputUtil import rasterValuesAtPoints
569	def QFun(x,y):
570	xll=domain.geo_reference.xllcorner
571	yll=domain.geo_reference.yllcorner
572	inDat=scipy.vstack([x+xll,y+yll]).transpose()
573	return rasterValuesAtPoints(xy=inDat,rasterFile=rasterFile,
574	interpolation=interpolation)
575
576	return QFun
577
578	#################################################################################
579
580	def quantity_from_Pt_Pol_Data_and_Raster(Pt_Pol_Data, quantity_raster, domain):
581	"""
582	Function to make a function F(x,y) which returns the corresponding
583	values on quantity_raster, except if x,y is inside the polygon associated with
584	any element of Pt_Pol_Data, in which case a Pt_Pol_-specific nearest neighbour
585	interpolator is used.
586
587	This has been superceeded by composite_quantity_setting_function
588
589	INPUT:
590	@param Pt_Pol_Data = a list with [ [ Polygon_0, Pt_XYZ_0],
591	[ Polygon_1, Pt_XYZ_1],
592	...
593	]
594	Here Polygon_i is a polygon in ANUGA format,
595	and Pt_XYZ_i is a 3 column array of x,y,Value points
596	@param quantity_raster = A GDAL-compatible quantity raster
597	@param domain = ANUGA domain
598	"""
599
600	# Function to set quantity from raster
601	qFun1=quantityRasterFun(domain, rasterFile=quantity_raster)
602
603	# List of function/polygon pairs defining the Pt_Pol_ quantity data
604	qFunChanList=[]
605	for i in range(len(Pt_Pol_Data)):
606	qFunChanList.append([
607	Pt_Pol_Data[i][0],
608	make_nearestNeighbour_quantity_function(Pt_Pol_Data[i][1], domain)
609	])
610
611	#
612	qFun=composite_quantity_setting_function(qFunChanList+[['All', qFun1]], domain)
613
614	return qFun

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

Download in other formats: