1 | #12345678901234567890123456789012345678901234567890123456789012345678901234567890 |
---|
2 | """This module contains various auxiliary function used by pyvolution. |
---|
3 | |
---|
4 | It is also a clearing house for functions that may later earn a module |
---|
5 | of their own. |
---|
6 | """ |
---|
7 | |
---|
8 | import anuga.utilities.polygon |
---|
9 | import sys |
---|
10 | import os |
---|
11 | |
---|
12 | from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,getcwd |
---|
13 | from os.path import exists, basename, split,join |
---|
14 | from warnings import warn |
---|
15 | from shutil import copy |
---|
16 | |
---|
17 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
18 | from Scientific.IO.NetCDF import NetCDFFile |
---|
19 | from Numeric import arange, choose, zeros, Float, array, allclose, take, compress |
---|
20 | |
---|
21 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
22 | from math import sqrt, atan, degrees |
---|
23 | |
---|
24 | # FIXME (Ole): Temporary short cuts - |
---|
25 | # FIXME (Ole): remove and update scripts where they are used |
---|
26 | from anuga.utilities.system_tools import get_revision_number |
---|
27 | from anuga.utilities.system_tools import store_version_info |
---|
28 | |
---|
29 | |
---|
30 | ## |
---|
31 | # @brief Read time history of data from NetCDF file, return callable object. |
---|
32 | # @param filename Name of .sww or .tms file. |
---|
33 | # @param domain Associated domain object. |
---|
34 | # @param quantities Name of quantity to be interpolated or a list of names. |
---|
35 | # @param interpolation_points List of absolute UTM coordinates for points |
---|
36 | # (N x 2) or geospatial object or |
---|
37 | # points file name at which values are sought. |
---|
38 | # @param time_thinning |
---|
39 | # @param verbose True if this function is to be verbose. |
---|
40 | # @param use_cache True means that caching of intermediate result is attempted. |
---|
41 | # @param boundary_polygon |
---|
42 | # @return A callable object. |
---|
43 | def file_function(filename, |
---|
44 | domain=None, |
---|
45 | quantities=None, |
---|
46 | interpolation_points=None, |
---|
47 | time_thinning=1, |
---|
48 | verbose=False, |
---|
49 | use_cache=False, |
---|
50 | boundary_polygon=None): |
---|
51 | """Read time history of spatial data from NetCDF file and return |
---|
52 | a callable object. |
---|
53 | |
---|
54 | Input variables: |
---|
55 | |
---|
56 | filename - Name of sww or tms file |
---|
57 | |
---|
58 | If the file has extension 'sww' then it is assumed to be spatio-temporal |
---|
59 | or temporal and the callable object will have the form f(t,x,y) or f(t) |
---|
60 | depending on whether the file contains spatial data |
---|
61 | |
---|
62 | If the file has extension 'tms' then it is assumed to be temporal only |
---|
63 | and the callable object will have the form f(t) |
---|
64 | |
---|
65 | Either form will return interpolated values based on the input file |
---|
66 | using the underlying interpolation_function. |
---|
67 | |
---|
68 | domain - Associated domain object |
---|
69 | If domain is specified, model time (domain.starttime) |
---|
70 | will be checked and possibly modified. |
---|
71 | |
---|
72 | All times are assumed to be in UTC |
---|
73 | |
---|
74 | All spatial information is assumed to be in absolute UTM coordinates. |
---|
75 | |
---|
76 | quantities - the name of the quantity to be interpolated or a |
---|
77 | list of quantity names. The resulting function will return |
---|
78 | a tuple of values - one for each quantity |
---|
79 | If quantities are None, the default quantities are |
---|
80 | ['stage', 'xmomentum', 'ymomentum'] |
---|
81 | |
---|
82 | |
---|
83 | interpolation_points - list of absolute UTM coordinates for points (N x 2) |
---|
84 | or geospatial object or points file name at which values are sought |
---|
85 | |
---|
86 | time_thinning - |
---|
87 | |
---|
88 | verbose - |
---|
89 | |
---|
90 | use_cache: True means that caching of intermediate result of |
---|
91 | Interpolation_function is attempted |
---|
92 | |
---|
93 | boundary_polygon - |
---|
94 | |
---|
95 | |
---|
96 | See Interpolation function in anuga.fit_interpolate.interpolation for |
---|
97 | further documentation |
---|
98 | """ |
---|
99 | |
---|
100 | # FIXME (OLE): Should check origin of domain against that of file |
---|
101 | # In fact, this is where origin should be converted to that of domain |
---|
102 | # Also, check that file covers domain fully. |
---|
103 | |
---|
104 | # Take into account: |
---|
105 | # - domain's georef |
---|
106 | # - sww file's georef |
---|
107 | # - interpolation points as absolute UTM coordinates |
---|
108 | |
---|
109 | if quantities is None: |
---|
110 | if verbose: |
---|
111 | msg = 'Quantities specified in file_function are None,' |
---|
112 | msg += ' so I will use stage, xmomentum, and ymomentum in that order' |
---|
113 | print msg |
---|
114 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
115 | |
---|
116 | # Use domain's startime if available |
---|
117 | if domain is not None: |
---|
118 | domain_starttime = domain.get_starttime() |
---|
119 | else: |
---|
120 | domain_starttime = None |
---|
121 | |
---|
122 | # Build arguments and keyword arguments for use with caching or apply. |
---|
123 | args = (filename,) |
---|
124 | |
---|
125 | # FIXME (Ole): Caching this function will not work well |
---|
126 | # if domain is passed in as instances change hash code. |
---|
127 | # Instead we pass in those attributes that are needed (and return them |
---|
128 | # if modified) |
---|
129 | kwargs = {'quantities': quantities, |
---|
130 | 'interpolation_points': interpolation_points, |
---|
131 | 'domain_starttime': domain_starttime, |
---|
132 | 'time_thinning': time_thinning, |
---|
133 | 'verbose': verbose, |
---|
134 | 'boundary_polygon': boundary_polygon} |
---|
135 | |
---|
136 | # Call underlying engine with or without caching |
---|
137 | if use_cache is True: |
---|
138 | try: |
---|
139 | from caching import cache |
---|
140 | except: |
---|
141 | msg = 'Caching was requested, but caching module'+\ |
---|
142 | 'could not be imported' |
---|
143 | raise msg |
---|
144 | |
---|
145 | f, starttime = cache(_file_function, |
---|
146 | args, kwargs, |
---|
147 | dependencies=[filename], |
---|
148 | compression=False, |
---|
149 | verbose=verbose) |
---|
150 | else: |
---|
151 | f, starttime = apply(_file_function, |
---|
152 | args, kwargs) |
---|
153 | |
---|
154 | #FIXME (Ole): Pass cache arguments, such as compression, in some sort of |
---|
155 | #structure |
---|
156 | |
---|
157 | f.starttime = starttime |
---|
158 | f.filename = filename |
---|
159 | |
---|
160 | if domain is not None: |
---|
161 | #Update domain.startime if it is *earlier* than starttime from file |
---|
162 | if starttime > domain.starttime: |
---|
163 | msg = 'WARNING: Start time as specified in domain (%f)'\ |
---|
164 | %domain.starttime |
---|
165 | msg += ' is earlier than the starttime of file %s (%f).'\ |
---|
166 | %(filename, starttime) |
---|
167 | msg += ' Modifying domain starttime accordingly.' |
---|
168 | |
---|
169 | if verbose: print msg |
---|
170 | |
---|
171 | domain.set_starttime(starttime) #Modifying model time |
---|
172 | |
---|
173 | if verbose: print 'Domain starttime is now set to %f'\ |
---|
174 | %domain.starttime |
---|
175 | return f |
---|
176 | |
---|
177 | |
---|
178 | ## |
---|
179 | # @brief ?? |
---|
180 | # @param filename Name of .sww or .tms file. |
---|
181 | # @param domain Associated domain object. |
---|
182 | # @param quantities Name of quantity to be interpolated or a list of names. |
---|
183 | # @param interpolation_points List of absolute UTM coordinates for points |
---|
184 | # (N x 2) or geospatial object or |
---|
185 | # points file name at which values are sought. |
---|
186 | # @param time_thinning |
---|
187 | # @param verbose True if this function is to be verbose. |
---|
188 | # @param use_cache True means that caching of intermediate result is attempted. |
---|
189 | # @param boundary_polygon |
---|
190 | def _file_function(filename, |
---|
191 | quantities=None, |
---|
192 | interpolation_points=None, |
---|
193 | domain_starttime=None, |
---|
194 | time_thinning=1, |
---|
195 | verbose=False, |
---|
196 | boundary_polygon=None): |
---|
197 | """Internal function |
---|
198 | |
---|
199 | See file_function for documentatiton |
---|
200 | """ |
---|
201 | |
---|
202 | assert type(filename) == type(''),\ |
---|
203 | 'First argument to File_function must be a string' |
---|
204 | |
---|
205 | try: |
---|
206 | fid = open(filename) |
---|
207 | except Exception, e: |
---|
208 | msg = 'File "%s" could not be opened: Error="%s"' % (filename, e) |
---|
209 | raise msg |
---|
210 | |
---|
211 | # read first line of file, guess file type |
---|
212 | line = fid.readline() |
---|
213 | fid.close() |
---|
214 | |
---|
215 | if line[:3] == 'CDF': |
---|
216 | return get_netcdf_file_function(filename, |
---|
217 | quantities, |
---|
218 | interpolation_points, |
---|
219 | domain_starttime, |
---|
220 | time_thinning=time_thinning, |
---|
221 | verbose=verbose, |
---|
222 | boundary_polygon=boundary_polygon) |
---|
223 | else: |
---|
224 | # FIXME (Ole): Could add csv file here to address Ted Rigby's |
---|
225 | # suggestion about reading hydrographs. |
---|
226 | # This may also deal with the gist of ticket:289 |
---|
227 | raise 'Must be a NetCDF File' |
---|
228 | |
---|
229 | |
---|
230 | ## |
---|
231 | # @brief ?? |
---|
232 | # @param filename Name of .sww or .tms file. |
---|
233 | # @param quantity_names Name of quantity to be interpolated or a list of names. |
---|
234 | # @param interpolation_points List of absolute UTM coordinates for points |
---|
235 | # (N x 2) or geospatial object or |
---|
236 | # points file name at which values are sought. |
---|
237 | # @param domain_starttime Start time from domain object. |
---|
238 | # @param time_thinning ?? |
---|
239 | # @param verbose True if this function is to be verbose. |
---|
240 | # @param boundary_polygon ?? |
---|
241 | # @return A callable object. |
---|
242 | def get_netcdf_file_function(filename, |
---|
243 | quantity_names=None, |
---|
244 | interpolation_points=None, |
---|
245 | domain_starttime=None, |
---|
246 | time_thinning=1, |
---|
247 | verbose=False, |
---|
248 | boundary_polygon=None): |
---|
249 | """Read time history of spatial data from NetCDF sww file and |
---|
250 | return a callable object f(t,x,y) |
---|
251 | which will return interpolated values based on the input file. |
---|
252 | |
---|
253 | Model time (domain_starttime) |
---|
254 | will be checked, possibly modified and returned |
---|
255 | |
---|
256 | All times are assumed to be in UTC |
---|
257 | |
---|
258 | See Interpolation function for further documetation |
---|
259 | """ |
---|
260 | |
---|
261 | # FIXME: Check that model origin is the same as file's origin |
---|
262 | # (both in UTM coordinates) |
---|
263 | # If not - modify those from file to match domain |
---|
264 | # (origin should be passed in) |
---|
265 | # Take this code from e.g. dem2pts in data_manager.py |
---|
266 | # FIXME: Use geo_reference to read and write xllcorner... |
---|
267 | |
---|
268 | import time, calendar, types |
---|
269 | from anuga.config import time_format |
---|
270 | from Numeric import array, zeros, Float, alltrue, concatenate, reshape |
---|
271 | |
---|
272 | # Open NetCDF file |
---|
273 | if verbose: print 'Reading', filename |
---|
274 | |
---|
275 | fid = NetCDFFile(filename, 'r') |
---|
276 | |
---|
277 | if type(quantity_names) == types.StringType: |
---|
278 | quantity_names = [quantity_names] |
---|
279 | |
---|
280 | if quantity_names is None or len(quantity_names) < 1: |
---|
281 | msg = 'No quantities are specified in file_function' |
---|
282 | raise Exception, msg |
---|
283 | |
---|
284 | if interpolation_points is not None: |
---|
285 | interpolation_points = ensure_absolute(interpolation_points) |
---|
286 | msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1] |
---|
287 | assert interpolation_points.shape[1] == 2, msg |
---|
288 | |
---|
289 | # Now assert that requested quantitites (and the independent ones) |
---|
290 | # are present in file |
---|
291 | missing = [] |
---|
292 | for quantity in ['time'] + quantity_names: |
---|
293 | if not fid.variables.has_key(quantity): |
---|
294 | missing.append(quantity) |
---|
295 | |
---|
296 | if len(missing) > 0: |
---|
297 | msg = 'Quantities %s could not be found in file %s'\ |
---|
298 | % (str(missing), filename) |
---|
299 | fid.close() |
---|
300 | raise Exception, msg |
---|
301 | |
---|
302 | # Decide whether this data has a spatial dimension |
---|
303 | spatial = True |
---|
304 | for quantity in ['x', 'y']: |
---|
305 | if not fid.variables.has_key(quantity): |
---|
306 | spatial = False |
---|
307 | |
---|
308 | if filename[-3:] == 'tms' and spatial is True: |
---|
309 | msg = 'Files of type tms must not contain spatial information' |
---|
310 | raise msg |
---|
311 | |
---|
312 | if filename[-3:] == 'sww' and spatial is False: |
---|
313 | msg = 'Files of type sww must contain spatial information' |
---|
314 | raise msg |
---|
315 | |
---|
316 | if filename[-3:] == 'sts' and spatial is False: |
---|
317 | #What if mux file only contains one point |
---|
318 | msg = 'Files of type sts must contain spatial information' |
---|
319 | raise msg |
---|
320 | |
---|
321 | if filename[-3:] == 'sts' and boundary_polygon is None: |
---|
322 | #What if mux file only contains one point |
---|
323 | msg = 'Files of type sts require boundary polygon' |
---|
324 | raise msg |
---|
325 | |
---|
326 | # Get first timestep |
---|
327 | try: |
---|
328 | starttime = fid.starttime[0] |
---|
329 | except ValueError: |
---|
330 | msg = 'Could not read starttime from file %s' %filename |
---|
331 | raise msg |
---|
332 | |
---|
333 | # Get variables |
---|
334 | # if verbose: print 'Get variables' |
---|
335 | time = fid.variables['time'][:] |
---|
336 | |
---|
337 | # Get time independent stuff |
---|
338 | if spatial: |
---|
339 | # Get origin |
---|
340 | xllcorner = fid.xllcorner[0] |
---|
341 | yllcorner = fid.yllcorner[0] |
---|
342 | zone = fid.zone[0] |
---|
343 | |
---|
344 | x = fid.variables['x'][:] |
---|
345 | y = fid.variables['y'][:] |
---|
346 | if filename[-3:] == 'sww': |
---|
347 | triangles = fid.variables['volumes'][:] |
---|
348 | |
---|
349 | x = reshape(x, (len(x),1)) |
---|
350 | y = reshape(y, (len(y),1)) |
---|
351 | vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array |
---|
352 | |
---|
353 | if boundary_polygon is not None: |
---|
354 | #removes sts points that do not lie on boundary |
---|
355 | boundary_polygon=ensure_numeric(boundary_polygon) |
---|
356 | boundary_polygon[:,0] -= xllcorner |
---|
357 | boundary_polygon[:,1] -= yllcorner |
---|
358 | temp=[] |
---|
359 | boundary_id=[] |
---|
360 | gauge_id=[] |
---|
361 | for i in range(len(boundary_polygon)): |
---|
362 | for j in range(len(x)): |
---|
363 | if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4): |
---|
364 | #FIX ME: |
---|
365 | #currently gauges lat and long is stored as float and |
---|
366 | #then cast to double. This cuases slight repositioning |
---|
367 | #of vertex_coordinates. |
---|
368 | temp.append(boundary_polygon[i]) |
---|
369 | gauge_id.append(j) |
---|
370 | boundary_id.append(i) |
---|
371 | break |
---|
372 | gauge_neighbour_id=[] |
---|
373 | for i in range(len(boundary_id)-1): |
---|
374 | if boundary_id[i]+1==boundary_id[i+1]: |
---|
375 | gauge_neighbour_id.append(i+1) |
---|
376 | else: |
---|
377 | gauge_neighbour_id.append(-1) |
---|
378 | if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \ |
---|
379 | and boundary_id[0]==0: |
---|
380 | gauge_neighbour_id.append(0) |
---|
381 | else: |
---|
382 | gauge_neighbour_id.append(-1) |
---|
383 | gauge_neighbour_id=ensure_numeric(gauge_neighbour_id) |
---|
384 | |
---|
385 | if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \ |
---|
386 | != len(temp)-1: |
---|
387 | msg='incorrect number of segments' |
---|
388 | raise msg |
---|
389 | vertex_coordinates=ensure_numeric(temp) |
---|
390 | if len(vertex_coordinates)==0: |
---|
391 | msg = 'None of the sts gauges fall on the boundary' |
---|
392 | raise msg |
---|
393 | else: |
---|
394 | gauge_neighbour_id=None |
---|
395 | |
---|
396 | if interpolation_points is not None: |
---|
397 | # Adjust for georef |
---|
398 | interpolation_points[:,0] -= xllcorner |
---|
399 | interpolation_points[:,1] -= yllcorner |
---|
400 | else: |
---|
401 | gauge_neighbour_id=None |
---|
402 | |
---|
403 | if domain_starttime is not None: |
---|
404 | # If domain_startime is *later* than starttime, |
---|
405 | # move time back - relative to domain's time |
---|
406 | if domain_starttime > starttime: |
---|
407 | time = time - domain_starttime + starttime |
---|
408 | |
---|
409 | # FIXME Use method in geo to reconcile |
---|
410 | # if spatial: |
---|
411 | # assert domain.geo_reference.xllcorner == xllcorner |
---|
412 | # assert domain.geo_reference.yllcorner == yllcorner |
---|
413 | # assert domain.geo_reference.zone == zone |
---|
414 | |
---|
415 | if verbose: |
---|
416 | print 'File_function data obtained from: %s' %filename |
---|
417 | print ' References:' |
---|
418 | #print ' Datum ....' #FIXME |
---|
419 | if spatial: |
---|
420 | print ' Lower left corner: [%f, %f]'\ |
---|
421 | %(xllcorner, yllcorner) |
---|
422 | print ' Start time: %f' %starttime |
---|
423 | |
---|
424 | |
---|
425 | # Produce values for desired data points at |
---|
426 | # each timestep for each quantity |
---|
427 | quantities = {} |
---|
428 | for i, name in enumerate(quantity_names): |
---|
429 | quantities[name] = fid.variables[name][:] |
---|
430 | if boundary_polygon is not None: |
---|
431 | #removes sts points that do not lie on boundary |
---|
432 | quantities[name] = take(quantities[name],gauge_id,1) |
---|
433 | |
---|
434 | # Close sww, tms or sts netcdf file |
---|
435 | fid.close() |
---|
436 | |
---|
437 | from anuga.fit_interpolate.interpolate import Interpolation_function |
---|
438 | |
---|
439 | if not spatial: |
---|
440 | vertex_coordinates = triangles = interpolation_points = None |
---|
441 | if filename[-3:] == 'sts':#added |
---|
442 | triangles = None |
---|
443 | #vertex coordinates is position of urs gauges |
---|
444 | |
---|
445 | # Return Interpolation_function instance as well as |
---|
446 | # starttime for use to possible modify that of domain |
---|
447 | return (Interpolation_function(time, |
---|
448 | quantities, |
---|
449 | quantity_names, |
---|
450 | vertex_coordinates, |
---|
451 | triangles, |
---|
452 | interpolation_points, |
---|
453 | time_thinning=time_thinning, |
---|
454 | verbose=verbose, |
---|
455 | gauge_neighbour_id=gauge_neighbour_id), |
---|
456 | starttime) |
---|
457 | |
---|
458 | # NOTE (Ole): Caching Interpolation function is too slow as |
---|
459 | # the very long parameters need to be hashed. |
---|
460 | |
---|
461 | |
---|
462 | ## |
---|
463 | # @brief Replace multiple substrings in a string. |
---|
464 | # @param text The string to operate on. |
---|
465 | # @param dictionary A dict containing replacements, key->value. |
---|
466 | # @return The new string. |
---|
467 | def multiple_replace(text, dictionary): |
---|
468 | """Multiple replace of words in text |
---|
469 | |
---|
470 | text: String to be modified |
---|
471 | dictionary: Mapping of words that are to be substituted |
---|
472 | |
---|
473 | Python Cookbook 3.14 page 88 and page 90 |
---|
474 | http://code.activestate.com/recipes/81330/ |
---|
475 | """ |
---|
476 | |
---|
477 | import re |
---|
478 | |
---|
479 | #Create a regular expression from all of the dictionary keys |
---|
480 | #matching only entire words |
---|
481 | regex = re.compile(r'\b'+ \ |
---|
482 | r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \ |
---|
483 | r'\b' ) |
---|
484 | |
---|
485 | #For each match, lookup the corresponding value in the dictionary |
---|
486 | return regex.sub(lambda match: dictionary[match.group(0)], text) |
---|
487 | |
---|
488 | |
---|
489 | ## |
---|
490 | # @brief Apply arbitrary expressions to the values of a dict. |
---|
491 | # @param expression A string expression to apply. |
---|
492 | # @param dictionary The dictionary to apply the expression to. |
---|
493 | def apply_expression_to_dictionary(expression, dictionary): |
---|
494 | """Apply arbitrary expression to values of dictionary |
---|
495 | |
---|
496 | Given an expression in terms of the keys, replace key by the |
---|
497 | corresponding values and evaluate. |
---|
498 | |
---|
499 | expression: Arbitrary, e.g. arithmetric, expression relating keys |
---|
500 | from dictionary. |
---|
501 | |
---|
502 | dictionary: Mapping between symbols in expression and objects that |
---|
503 | will be evaluated by expression. |
---|
504 | Values in dictionary must support operators given in |
---|
505 | expression e.g. by overloading |
---|
506 | |
---|
507 | Due to a limitation with Numeric, this can not evaluate 0/0 |
---|
508 | In general, the user can fix by adding 1e-30 to the numerator. |
---|
509 | SciPy core can handle this situation. |
---|
510 | """ |
---|
511 | |
---|
512 | import types |
---|
513 | import re |
---|
514 | |
---|
515 | assert isinstance(expression, basestring) |
---|
516 | assert type(dictionary) == types.DictType |
---|
517 | |
---|
518 | #Convert dictionary values to textual representations suitable for eval |
---|
519 | D = {} |
---|
520 | for key in dictionary: |
---|
521 | D[key] = 'dictionary["%s"]' % key |
---|
522 | |
---|
523 | #Perform substitution of variables |
---|
524 | expression = multiple_replace(expression, D) |
---|
525 | |
---|
526 | #Evaluate and return |
---|
527 | try: |
---|
528 | return eval(expression) |
---|
529 | except NameError, e: |
---|
530 | msg = 'Expression "%s" could not be evaluated: %s' % (expression, e) |
---|
531 | raise NameError, msg |
---|
532 | except ValueError, e: |
---|
533 | msg = 'Expression "%s" could not be evaluated: %s' % (expression, e) |
---|
534 | raise ValueError, msg |
---|
535 | |
---|
536 | |
---|
537 | ## |
---|
538 | # @brief Format a float into a string. |
---|
539 | # @param value Float value to format. |
---|
540 | # @param format The format to use (%.2f is default). |
---|
541 | # @return The formatted float as a string. |
---|
542 | def get_textual_float(value, format = '%.2f'): |
---|
543 | """Get textual representation of floating point numbers |
---|
544 | and accept None as valid entry |
---|
545 | |
---|
546 | format is a string - default = '%.2f' |
---|
547 | """ |
---|
548 | |
---|
549 | if value is None: |
---|
550 | return 'None' |
---|
551 | else: |
---|
552 | try: |
---|
553 | float(value) |
---|
554 | except: |
---|
555 | # May this is a vector |
---|
556 | if len(value) > 1: |
---|
557 | s = '(' |
---|
558 | for v in value: |
---|
559 | s += get_textual_float(v, format) + ', ' |
---|
560 | |
---|
561 | s = s[:-2] + ')' # Strip trailing comma and close |
---|
562 | return s |
---|
563 | else: |
---|
564 | raise 'Illegal input to get_textual_float:', value |
---|
565 | else: |
---|
566 | return format % float(value) |
---|
567 | |
---|
568 | |
---|
569 | ################################################################################# |
---|
570 | # OBSOLETE STUFF |
---|
571 | ################################################################################# |
---|
572 | |
---|
573 | # @note TEMP |
---|
574 | def angle(v1, v2): |
---|
575 | """Temporary Interface to new location""" |
---|
576 | |
---|
577 | import anuga.utilities.numerical_tools as NT |
---|
578 | |
---|
579 | msg = 'angle has moved from util.py. ' |
---|
580 | msg += 'Please use "from anuga.utilities.numerical_tools import angle"' |
---|
581 | warn(msg, DeprecationWarning) |
---|
582 | |
---|
583 | return NT.angle(v1, v2) |
---|
584 | |
---|
585 | |
---|
586 | # @note TEMP |
---|
587 | def anglediff(v0, v1): |
---|
588 | """Temporary Interface to new location""" |
---|
589 | |
---|
590 | import anuga.utilities.numerical_tools as NT |
---|
591 | |
---|
592 | msg = 'anglediff has moved from util.py. ' |
---|
593 | msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"' |
---|
594 | warn(msg, DeprecationWarning) |
---|
595 | |
---|
596 | return NT.anglediff(v0, v1) |
---|
597 | |
---|
598 | |
---|
599 | # @note TEMP |
---|
600 | def mean(x): |
---|
601 | """Temporary Interface to new location""" |
---|
602 | |
---|
603 | import anuga.utilities.numerical_tools as NT |
---|
604 | |
---|
605 | msg = 'mean has moved from util.py. ' |
---|
606 | msg += 'Please use "from anuga.utilities.numerical_tools import mean"' |
---|
607 | warn(msg, DeprecationWarning) |
---|
608 | |
---|
609 | return NT.mean(x) |
---|
610 | |
---|
611 | |
---|
612 | # @note TEMP |
---|
613 | def point_on_line(*args, **kwargs): |
---|
614 | """Temporary Interface to new location""" |
---|
615 | |
---|
616 | msg = 'point_on_line has moved from util.py. ' |
---|
617 | msg += 'Please use "from anuga.utilities.polygon import point_on_line"' |
---|
618 | warn(msg, DeprecationWarning) |
---|
619 | |
---|
620 | return utilities.polygon.point_on_line(*args, **kwargs) |
---|
621 | |
---|
622 | |
---|
623 | # @note TEMP |
---|
624 | def inside_polygon(*args, **kwargs): |
---|
625 | """Temporary Interface to new location""" |
---|
626 | |
---|
627 | print 'inside_polygon has moved from util.py. ', |
---|
628 | print 'Please use "from anuga.utilities.polygon import inside_polygon"' |
---|
629 | |
---|
630 | return utilities.polygon.inside_polygon(*args, **kwargs) |
---|
631 | |
---|
632 | |
---|
633 | # @note TEMP |
---|
634 | def outside_polygon(*args, **kwargs): |
---|
635 | """Temporary Interface to new location""" |
---|
636 | |
---|
637 | print 'outside_polygon has moved from util.py. ', |
---|
638 | print 'Please use "from anuga.utilities.polygon import outside_polygon"' |
---|
639 | |
---|
640 | return utilities.polygon.outside_polygon(*args, **kwargs) |
---|
641 | |
---|
642 | |
---|
643 | # @note TEMP |
---|
644 | def separate_points_by_polygon(*args, **kwargs): |
---|
645 | """Temporary Interface to new location""" |
---|
646 | |
---|
647 | print 'separate_points_by_polygon has moved from util.py. ', |
---|
648 | print 'Please use "from anuga.utilities.polygon import ' \ |
---|
649 | 'separate_points_by_polygon"' |
---|
650 | |
---|
651 | return utilities.polygon.separate_points_by_polygon(*args, **kwargs) |
---|
652 | |
---|
653 | |
---|
654 | # @note TEMP |
---|
655 | def read_polygon(*args, **kwargs): |
---|
656 | """Temporary Interface to new location""" |
---|
657 | |
---|
658 | print 'read_polygon has moved from util.py. ', |
---|
659 | print 'Please use "from anuga.utilities.polygon import read_polygon"' |
---|
660 | |
---|
661 | return utilities.polygon.read_polygon(*args, **kwargs) |
---|
662 | |
---|
663 | |
---|
664 | # @note TEMP |
---|
665 | def populate_polygon(*args, **kwargs): |
---|
666 | """Temporary Interface to new location""" |
---|
667 | |
---|
668 | print 'populate_polygon has moved from util.py. ', |
---|
669 | print 'Please use "from anuga.utilities.polygon import populate_polygon"' |
---|
670 | |
---|
671 | return utilities.polygon.populate_polygon(*args, **kwargs) |
---|
672 | |
---|
673 | |
---|
674 | ################################################################################# |
---|
675 | # End of obsolete stuff ? |
---|
676 | ################################################################################# |
---|
677 | |
---|
678 | # @note TEMP |
---|
679 | def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', |
---|
680 | verbose=False): |
---|
681 | """Temporary Interface to new location""" |
---|
682 | from anuga.shallow_water.data_manager import start_screen_catcher \ |
---|
683 | as dm_start_screen_catcher |
---|
684 | |
---|
685 | print 'start_screen_catcher has moved from util.py. ', |
---|
686 | print 'Please use "from anuga.shallow_water.data_manager import ' \ |
---|
687 | 'start_screen_catcher"' |
---|
688 | |
---|
689 | return dm_start_screen_catcher(dir_name, myid='', numprocs='', |
---|
690 | extra_info='', verbose=False) |
---|
691 | |
---|
692 | |
---|
693 | ## |
---|
694 | # @brief Read a .sww file and plot the time series. |
---|
695 | # @param swwfiles Dictionary of .sww files. |
---|
696 | # @param gauge_filename Name of gauge file. |
---|
697 | # @param production_dirs ?? |
---|
698 | # @param report If True, write figures to report directory. |
---|
699 | # @param reportname Name of generated report (if any). |
---|
700 | # @param plot_quantity List containing quantities to plot. |
---|
701 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
702 | # @param surface If True, then generate solution surface with 3d plot. |
---|
703 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
704 | # @param time_max End of user defined time range for plotting purposes. |
---|
705 | # @param time_thinning ?? |
---|
706 | # @param time_unit ?? |
---|
707 | # @param title_on If True, export standard graphics with title. |
---|
708 | # @param use_cache If True, use caching. |
---|
709 | # @param verbose If True, this function is verbose. |
---|
710 | def sww2timeseries(swwfiles, |
---|
711 | gauge_filename, |
---|
712 | production_dirs, |
---|
713 | report=None, |
---|
714 | reportname=None, |
---|
715 | plot_quantity=None, |
---|
716 | generate_fig=False, |
---|
717 | surface=None, |
---|
718 | time_min=None, |
---|
719 | time_max=None, |
---|
720 | time_thinning=1, |
---|
721 | time_unit=None, |
---|
722 | title_on=None, |
---|
723 | use_cache=False, |
---|
724 | verbose=False): |
---|
725 | """ Read sww file and plot the time series for the |
---|
726 | prescribed quantities at defined gauge locations and |
---|
727 | prescribed time range. |
---|
728 | |
---|
729 | Input variables: |
---|
730 | |
---|
731 | swwfiles - dictionary of sww files with label_ids (used in |
---|
732 | generating latex output. It will be part of |
---|
733 | the directory name of file_loc (typically the timestamp). |
---|
734 | Helps to differentiate latex files for different |
---|
735 | simulations for a particular scenario. |
---|
736 | - assume that all conserved quantities have been stored |
---|
737 | - assume each sww file has been simulated with same timestep |
---|
738 | |
---|
739 | gauge_filename - name of file containing gauge data |
---|
740 | - easting, northing, name , elevation? |
---|
741 | - OR (this is not yet done) |
---|
742 | - structure which can be converted to a Numeric array, |
---|
743 | such as a geospatial data object |
---|
744 | |
---|
745 | production_dirs - A list of list, example {20061101_121212: '1 in 10000', |
---|
746 | 'boundaries': 'urs boundary'} |
---|
747 | this will use the second part as the label and the |
---|
748 | first part as the ? |
---|
749 | #FIXME: Is it a list or a dictionary |
---|
750 | # This is probably obsolete by now |
---|
751 | |
---|
752 | report - if True, then write figures to report_figures directory in |
---|
753 | relevant production directory |
---|
754 | - if False, figures are already stored with sww file |
---|
755 | - default to False |
---|
756 | |
---|
757 | reportname - name for report if wishing to generate report |
---|
758 | |
---|
759 | plot_quantity - list containing quantities to plot, they must |
---|
760 | be the name of an existing quantity or one of |
---|
761 | the following possibilities |
---|
762 | - possibilities: |
---|
763 | - stage; 'stage' |
---|
764 | - depth; 'depth' |
---|
765 | - speed; calculated as absolute momentum |
---|
766 | (pointwise) divided by depth; 'speed' |
---|
767 | - bearing; calculated as the angle of the momentum |
---|
768 | vector (xmomentum, ymomentum) from the North; 'bearing' |
---|
769 | - absolute momentum; calculated as |
---|
770 | sqrt(xmomentum^2 + ymomentum^2); 'momentum' |
---|
771 | - x momentum; 'xmomentum' |
---|
772 | - y momentum; 'ymomentum' |
---|
773 | - default will be ['stage', 'speed', 'bearing'] |
---|
774 | |
---|
775 | generate_fig - if True, generate figures as well as csv file |
---|
776 | - if False, csv files created only |
---|
777 | |
---|
778 | surface - if True, then generate solution surface with 3d plot |
---|
779 | and save to current working directory |
---|
780 | - default = False |
---|
781 | |
---|
782 | time_min - beginning of user defined time range for plotting purposes |
---|
783 | - default will be first available time found in swwfile |
---|
784 | |
---|
785 | time_max - end of user defined time range for plotting purposes |
---|
786 | - default will be last available time found in swwfile |
---|
787 | |
---|
788 | title_on - if True, export standard graphics with title |
---|
789 | - if False, export standard graphics without title |
---|
790 | |
---|
791 | |
---|
792 | Output: |
---|
793 | |
---|
794 | - time series data stored in .csv for later use if required. |
---|
795 | Name = gauges_timeseries followed by gauge name |
---|
796 | - latex file will be generated in same directory as where script is |
---|
797 | run (usually production scenario directory. |
---|
798 | Name = latexoutputlabel_id.tex |
---|
799 | |
---|
800 | Other important information: |
---|
801 | |
---|
802 | It is assumed that the used has stored all the conserved quantities |
---|
803 | and elevation during the scenario run, i.e. |
---|
804 | ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
805 | If this has not occurred then sww2timeseries will not work. |
---|
806 | |
---|
807 | |
---|
808 | Usage example |
---|
809 | texname = sww2timeseries({project.boundary_name + '.sww': ''}, |
---|
810 | project.polygons_dir + sep + 'boundary_extent.csv', |
---|
811 | project.anuga_dir, |
---|
812 | report = False, |
---|
813 | plot_quantity = ['stage', 'speed', 'bearing'], |
---|
814 | time_min = None, |
---|
815 | time_max = None, |
---|
816 | title_on = True, |
---|
817 | verbose = True) |
---|
818 | |
---|
819 | """ |
---|
820 | |
---|
821 | msg = 'NOTE: A new function is available to create csv files from sww ' |
---|
822 | msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util' |
---|
823 | msg += ' PLUS another new function to create graphs from csv files called ' |
---|
824 | msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util' |
---|
825 | print msg |
---|
826 | |
---|
827 | k = _sww2timeseries(swwfiles, |
---|
828 | gauge_filename, |
---|
829 | production_dirs, |
---|
830 | report, |
---|
831 | reportname, |
---|
832 | plot_quantity, |
---|
833 | generate_fig, |
---|
834 | surface, |
---|
835 | time_min, |
---|
836 | time_max, |
---|
837 | time_thinning, |
---|
838 | time_unit, |
---|
839 | title_on, |
---|
840 | use_cache, |
---|
841 | verbose) |
---|
842 | return k |
---|
843 | |
---|
844 | |
---|
845 | ## |
---|
846 | # @brief Read a .sww file and plot the time series. |
---|
847 | # @param swwfiles Dictionary of .sww files. |
---|
848 | # @param gauge_filename Name of gauge file. |
---|
849 | # @param production_dirs ?? |
---|
850 | # @param report If True, write figures to report directory. |
---|
851 | # @param reportname Name of generated report (if any). |
---|
852 | # @param plot_quantity List containing quantities to plot. |
---|
853 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
854 | # @param surface If True, then generate solution surface with 3d plot. |
---|
855 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
856 | # @param time_max End of user defined time range for plotting purposes. |
---|
857 | # @param time_thinning ?? |
---|
858 | # @param time_unit ?? |
---|
859 | # @param title_on If True, export standard graphics with title. |
---|
860 | # @param use_cache If True, use caching. |
---|
861 | # @param verbose If True, this function is verbose. |
---|
862 | def _sww2timeseries(swwfiles, |
---|
863 | gauge_filename, |
---|
864 | production_dirs, |
---|
865 | report = None, |
---|
866 | reportname = None, |
---|
867 | plot_quantity = None, |
---|
868 | generate_fig = False, |
---|
869 | surface = None, |
---|
870 | time_min = None, |
---|
871 | time_max = None, |
---|
872 | time_thinning = 1, |
---|
873 | time_unit = None, |
---|
874 | title_on = None, |
---|
875 | use_cache = False, |
---|
876 | verbose = False): |
---|
877 | |
---|
878 | assert type(gauge_filename) == type(''), 'Gauge filename must be a string' |
---|
879 | |
---|
880 | try: |
---|
881 | fid = open(gauge_filename) |
---|
882 | except Exception, e: |
---|
883 | msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e) |
---|
884 | raise msg |
---|
885 | |
---|
886 | if report is None: |
---|
887 | report = False |
---|
888 | |
---|
889 | if plot_quantity is None: |
---|
890 | plot_quantity = ['depth', 'speed'] |
---|
891 | else: |
---|
892 | assert type(plot_quantity) == list, 'plot_quantity must be a list' |
---|
893 | check_list(plot_quantity) |
---|
894 | |
---|
895 | if surface is None: |
---|
896 | surface = False |
---|
897 | |
---|
898 | if time_unit is None: |
---|
899 | time_unit = 'hours' |
---|
900 | |
---|
901 | if title_on is None: |
---|
902 | title_on = True |
---|
903 | |
---|
904 | if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename |
---|
905 | |
---|
906 | gauges, locations, elev = get_gauges_from_file(gauge_filename) |
---|
907 | |
---|
908 | sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
909 | |
---|
910 | file_loc = [] |
---|
911 | f_list = [] |
---|
912 | label_id = [] |
---|
913 | leg_label = [] |
---|
914 | themaxT = 0.0 |
---|
915 | theminT = 0.0 |
---|
916 | |
---|
917 | for swwfile in swwfiles.keys(): |
---|
918 | try: |
---|
919 | fid = open(swwfile) |
---|
920 | except Exception, e: |
---|
921 | msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e) |
---|
922 | raise msg |
---|
923 | |
---|
924 | print 'swwfile', swwfile |
---|
925 | |
---|
926 | # Extract parent dir name and use as label |
---|
927 | path, _ = os.path.split(swwfile) |
---|
928 | _, label = os.path.split(path) |
---|
929 | |
---|
930 | #print 'label', label |
---|
931 | leg_label.append(label) |
---|
932 | |
---|
933 | f = file_function(swwfile, |
---|
934 | quantities = sww_quantity, |
---|
935 | interpolation_points = gauges, |
---|
936 | time_thinning = time_thinning, |
---|
937 | verbose = verbose, |
---|
938 | use_cache = use_cache) |
---|
939 | |
---|
940 | # determine which gauges are contained in sww file |
---|
941 | count = 0 |
---|
942 | gauge_index = [] |
---|
943 | for k, g in enumerate(gauges): |
---|
944 | if f(0.0, point_id = k)[2] > 1.0e6: |
---|
945 | count += 1 |
---|
946 | if count == 1: print 'Gauges not contained here:' |
---|
947 | print locations[k] |
---|
948 | else: |
---|
949 | gauge_index.append(k) |
---|
950 | |
---|
951 | if len(gauge_index) > 0: |
---|
952 | print 'Gauges contained here: \n', |
---|
953 | else: |
---|
954 | print 'No gauges contained here. \n' |
---|
955 | for i in range(len(gauge_index)): |
---|
956 | print locations[gauge_index[i]] |
---|
957 | |
---|
958 | index = swwfile.rfind(sep) |
---|
959 | file_loc.append(swwfile[:index+1]) |
---|
960 | label_id.append(swwfiles[swwfile]) |
---|
961 | |
---|
962 | f_list.append(f) |
---|
963 | maxT = max(f.get_time()) |
---|
964 | minT = min(f.get_time()) |
---|
965 | if maxT > themaxT: themaxT = maxT |
---|
966 | if minT > theminT: theminT = minT |
---|
967 | |
---|
968 | if time_min is None: |
---|
969 | time_min = theminT # min(T) |
---|
970 | else: |
---|
971 | if time_min < theminT: # min(T): |
---|
972 | msg = 'Minimum time entered not correct - please try again' |
---|
973 | raise Exception, msg |
---|
974 | |
---|
975 | if time_max is None: |
---|
976 | time_max = themaxT # max(T) |
---|
977 | else: |
---|
978 | if time_max > themaxT: # max(T): |
---|
979 | msg = 'Maximum time entered not correct - please try again' |
---|
980 | raise Exception, msg |
---|
981 | |
---|
982 | if verbose and len(gauge_index) > 0: |
---|
983 | print 'Inputs OK - going to generate figures' |
---|
984 | |
---|
985 | if len(gauge_index) <> 0: |
---|
986 | texfile, elev_output = \ |
---|
987 | generate_figures(plot_quantity, file_loc, report, reportname, |
---|
988 | surface, leg_label, f_list, gauges, locations, |
---|
989 | elev, gauge_index, production_dirs, time_min, |
---|
990 | time_max, time_unit, title_on, label_id, |
---|
991 | generate_fig, verbose) |
---|
992 | else: |
---|
993 | texfile = '' |
---|
994 | elev_output = [] |
---|
995 | |
---|
996 | return texfile, elev_output |
---|
997 | |
---|
998 | |
---|
999 | ## |
---|
1000 | # @brief Read gauge info from a file. |
---|
1001 | # @param filename The name of the file to read. |
---|
1002 | # @return A (gauges, gaugelocation, elev) tuple. |
---|
1003 | def get_gauges_from_file(filename): |
---|
1004 | """ Read in gauge information from file |
---|
1005 | """ |
---|
1006 | |
---|
1007 | from os import sep, getcwd, access, F_OK, mkdir |
---|
1008 | |
---|
1009 | # Get data from the gauge file |
---|
1010 | fid = open(filename) |
---|
1011 | lines = fid.readlines() |
---|
1012 | fid.close() |
---|
1013 | |
---|
1014 | gauges = [] |
---|
1015 | gaugelocation = [] |
---|
1016 | elev = [] |
---|
1017 | |
---|
1018 | # Check header information |
---|
1019 | line1 = lines[0] |
---|
1020 | line11 = line1.split(',') |
---|
1021 | |
---|
1022 | if isinstance(line11[0], str) is True: |
---|
1023 | # We have found text in the first line |
---|
1024 | east_index = None |
---|
1025 | north_index = None |
---|
1026 | name_index = None |
---|
1027 | elev_index = None |
---|
1028 | |
---|
1029 | for i in range(len(line11)): |
---|
1030 | if line11[i].strip().lower() == 'easting': east_index = i |
---|
1031 | if line11[i].strip().lower() == 'northing': north_index = i |
---|
1032 | if line11[i].strip().lower() == 'name': name_index = i |
---|
1033 | if line11[i].strip().lower() == 'elevation': elev_index = i |
---|
1034 | |
---|
1035 | if east_index < len(line11) and north_index < len(line11): |
---|
1036 | pass |
---|
1037 | else: |
---|
1038 | msg = 'WARNING: %s does not contain correct header information' \ |
---|
1039 | % filename |
---|
1040 | msg += 'The header must be: easting, northing, name, elevation' |
---|
1041 | raise Exception, msg |
---|
1042 | |
---|
1043 | if elev_index is None: |
---|
1044 | raise Exception |
---|
1045 | |
---|
1046 | if name_index is None: |
---|
1047 | raise Exception |
---|
1048 | |
---|
1049 | lines = lines[1:] # Remove header from data |
---|
1050 | else: |
---|
1051 | # No header, assume that this is a simple easting, northing file |
---|
1052 | |
---|
1053 | msg = 'There was no header in file %s and the number of columns is %d' \ |
---|
1054 | % (filename, len(line11)) |
---|
1055 | msg += '- was assuming two columns corresponding to Easting and Northing' |
---|
1056 | assert len(line11) == 2, msg |
---|
1057 | |
---|
1058 | east_index = 0 |
---|
1059 | north_index = 1 |
---|
1060 | |
---|
1061 | N = len(lines) |
---|
1062 | elev = [-9999]*N |
---|
1063 | gaugelocation = range(N) |
---|
1064 | |
---|
1065 | # Read in gauge data |
---|
1066 | for line in lines: |
---|
1067 | fields = line.split(',') |
---|
1068 | |
---|
1069 | gauges.append([float(fields[east_index]), float(fields[north_index])]) |
---|
1070 | |
---|
1071 | if len(fields) > 2: |
---|
1072 | elev.append(float(fields[elev_index])) |
---|
1073 | loc = fields[name_index] |
---|
1074 | gaugelocation.append(loc.strip('\n')) |
---|
1075 | |
---|
1076 | return gauges, gaugelocation, elev |
---|
1077 | |
---|
1078 | |
---|
1079 | ## |
---|
1080 | # @brief Check that input quantities in quantity list are legal. |
---|
1081 | # @param quantity Quantity list to check. |
---|
1082 | # @note Raises an exception of list is not legal. |
---|
1083 | def check_list(quantity): |
---|
1084 | """ Check that input quantities in quantity list are possible |
---|
1085 | """ |
---|
1086 | import sys |
---|
1087 | from sets import Set as set |
---|
1088 | |
---|
1089 | all_quantity = ['stage', 'depth', 'momentum', 'xmomentum', |
---|
1090 | 'ymomentum', 'speed', 'bearing', 'elevation'] |
---|
1091 | |
---|
1092 | # convert all quanitiy names to lowercase |
---|
1093 | for i,j in enumerate(quantity): |
---|
1094 | quantity[i] = quantity[i].lower() |
---|
1095 | |
---|
1096 | # check that all names in 'quantity' appear in 'all_quantity' |
---|
1097 | p = list(set(quantity).difference(set(all_quantity))) |
---|
1098 | if len(p) != 0: |
---|
1099 | msg = 'Quantities %s do not exist - please try again' %p |
---|
1100 | raise Exception, msg |
---|
1101 | |
---|
1102 | |
---|
1103 | ## |
---|
1104 | # @brief Calculate velocity bearing from North. |
---|
1105 | # @param uh ?? |
---|
1106 | # @param vh ?? |
---|
1107 | # @return The calculated bearing. |
---|
1108 | def calc_bearing(uh, vh): |
---|
1109 | """ Calculate velocity bearing from North |
---|
1110 | """ |
---|
1111 | #FIXME (Ole): I reckon we should refactor this one to use |
---|
1112 | # the function angle() in utilities/numerical_tools |
---|
1113 | # |
---|
1114 | # It will be a simple matter of |
---|
1115 | # * converting from radians to degrees |
---|
1116 | # * moving the reference direction from [1,0] to North |
---|
1117 | # * changing from counter clockwise to clocwise. |
---|
1118 | |
---|
1119 | angle = degrees(atan(vh/(uh+1.e-15))) |
---|
1120 | |
---|
1121 | if (0 < angle < 90.0): |
---|
1122 | if vh > 0: |
---|
1123 | bearing = 90.0 - abs(angle) |
---|
1124 | if vh < 0: |
---|
1125 | bearing = 270.0 - abs(angle) |
---|
1126 | |
---|
1127 | if (-90 < angle < 0): |
---|
1128 | if vh < 0: |
---|
1129 | bearing = 90.0 - (angle) |
---|
1130 | if vh > 0: |
---|
1131 | bearing = 270.0 - (angle) |
---|
1132 | if angle == 0: bearing = 0.0 |
---|
1133 | |
---|
1134 | return bearing |
---|
1135 | |
---|
1136 | |
---|
1137 | ## |
---|
1138 | # @brief Generate figures from quantities and gauges for each sww file. |
---|
1139 | # @param plot_quantity ?? |
---|
1140 | # @param file_loc ?? |
---|
1141 | # @param report ?? |
---|
1142 | # @param reportname ?? |
---|
1143 | # @param surface ?? |
---|
1144 | # @param leg_label ?? |
---|
1145 | # @param f_list ?? |
---|
1146 | # @param gauges ?? |
---|
1147 | # @param locations ?? |
---|
1148 | # @param elev ?? |
---|
1149 | # @param gauge_index ?? |
---|
1150 | # @param production_dirs ?? |
---|
1151 | # @param time_min ?? |
---|
1152 | # @param time_max ?? |
---|
1153 | # @param time_unit ?? |
---|
1154 | # @param title_on ?? |
---|
1155 | # @param label_id ?? |
---|
1156 | # @param generate_fig ?? |
---|
1157 | # @param verbose?? |
---|
1158 | # @return (texfile2, elev_output) |
---|
1159 | def generate_figures(plot_quantity, file_loc, report, reportname, surface, |
---|
1160 | leg_label, f_list, gauges, locations, elev, gauge_index, |
---|
1161 | production_dirs, time_min, time_max, time_unit, |
---|
1162 | title_on, label_id, generate_fig, verbose): |
---|
1163 | """ Generate figures based on required quantities and gauges for |
---|
1164 | each sww file |
---|
1165 | """ |
---|
1166 | from Numeric import ones, allclose, zeros, Float, ravel |
---|
1167 | from os import sep, altsep, getcwd, mkdir, access, F_OK, environ |
---|
1168 | |
---|
1169 | if generate_fig is True: |
---|
1170 | from pylab import ion, hold, plot, axis, figure, legend, savefig, \ |
---|
1171 | xlabel, ylabel, title, close, subplot |
---|
1172 | |
---|
1173 | if surface is True: |
---|
1174 | import pylab as p1 |
---|
1175 | import mpl3d.mplot3d as p3 |
---|
1176 | |
---|
1177 | if report == True: |
---|
1178 | texdir = getcwd()+sep+'report'+sep |
---|
1179 | if access(texdir,F_OK) == 0: |
---|
1180 | mkdir (texdir) |
---|
1181 | if len(label_id) == 1: |
---|
1182 | label_id1 = label_id[0].replace(sep,'') |
---|
1183 | label_id2 = label_id1.replace('_','') |
---|
1184 | texfile = texdir + reportname + '%s' % label_id2 |
---|
1185 | texfile2 = reportname + '%s' % label_id2 |
---|
1186 | texfilename = texfile + '.tex' |
---|
1187 | fid = open(texfilename, 'w') |
---|
1188 | |
---|
1189 | if verbose: print '\n Latex output printed to %s \n' %texfilename |
---|
1190 | else: |
---|
1191 | texfile = texdir+reportname |
---|
1192 | texfile2 = reportname |
---|
1193 | texfilename = texfile + '.tex' |
---|
1194 | fid = open(texfilename, 'w') |
---|
1195 | |
---|
1196 | if verbose: print '\n Latex output printed to %s \n' %texfilename |
---|
1197 | else: |
---|
1198 | texfile = '' |
---|
1199 | texfile2 = '' |
---|
1200 | |
---|
1201 | p = len(f_list) |
---|
1202 | n = [] |
---|
1203 | n0 = 0 |
---|
1204 | for i in range(len(f_list)): |
---|
1205 | n.append(len(f_list[i].get_time())) |
---|
1206 | if n[i] > n0: n0 = n[i] |
---|
1207 | n0 = int(n0) |
---|
1208 | m = len(locations) |
---|
1209 | model_time = zeros((n0, m, p), Float) |
---|
1210 | stages = zeros((n0, m, p), Float) |
---|
1211 | elevations = zeros((n0, m, p), Float) |
---|
1212 | momenta = zeros((n0, m, p), Float) |
---|
1213 | xmom = zeros((n0, m, p), Float) |
---|
1214 | ymom = zeros((n0, m, p), Float) |
---|
1215 | speed = zeros((n0, m, p), Float) |
---|
1216 | bearings = zeros((n0, m, p), Float) |
---|
1217 | due_east = 90.0*ones((n0, 1), Float) |
---|
1218 | due_west = 270.0*ones((n0, 1), Float) |
---|
1219 | depths = zeros((n0, m, p), Float) |
---|
1220 | eastings = zeros((n0, m, p), Float) |
---|
1221 | min_stages = [] |
---|
1222 | max_stages = [] |
---|
1223 | min_momentums = [] |
---|
1224 | max_momentums = [] |
---|
1225 | max_xmomentums = [] |
---|
1226 | max_ymomentums = [] |
---|
1227 | min_xmomentums = [] |
---|
1228 | min_ymomentums = [] |
---|
1229 | max_speeds = [] |
---|
1230 | min_speeds = [] |
---|
1231 | max_depths = [] |
---|
1232 | model_time_plot3d = zeros((n0, m), Float) |
---|
1233 | stages_plot3d = zeros((n0, m), Float) |
---|
1234 | eastings_plot3d = zeros((n0, m),Float) |
---|
1235 | if time_unit is 'mins': scale = 60.0 |
---|
1236 | if time_unit is 'hours': scale = 3600.0 |
---|
1237 | |
---|
1238 | ##### loop over each swwfile ##### |
---|
1239 | for j, f in enumerate(f_list): |
---|
1240 | if verbose: print 'swwfile %d of %d' % (j, len(f_list)) |
---|
1241 | |
---|
1242 | starttime = f.starttime |
---|
1243 | comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv' |
---|
1244 | fid_compare = open(comparefile, 'w') |
---|
1245 | file0 = file_loc[j] + 'gauges_t0.csv' |
---|
1246 | fid_0 = open(file0, 'w') |
---|
1247 | |
---|
1248 | ##### loop over each gauge ##### |
---|
1249 | for k in gauge_index: |
---|
1250 | if verbose: print 'Gauge %d of %d' % (k, len(gauges)) |
---|
1251 | |
---|
1252 | g = gauges[k] |
---|
1253 | min_stage = 10 |
---|
1254 | max_stage = 0 |
---|
1255 | max_momentum = max_xmomentum = max_ymomentum = 0 |
---|
1256 | min_momentum = min_xmomentum = min_ymomentum = 100 |
---|
1257 | max_speed = 0 |
---|
1258 | min_speed = 0 |
---|
1259 | max_depth = 0 |
---|
1260 | gaugeloc = str(locations[k]) |
---|
1261 | thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \ |
---|
1262 | + gaugeloc + '.csv' |
---|
1263 | fid_out = open(thisfile, 'w') |
---|
1264 | s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' |
---|
1265 | fid_out.write(s) |
---|
1266 | |
---|
1267 | #### generate quantities ####### |
---|
1268 | for i, t in enumerate(f.get_time()): |
---|
1269 | if time_min <= t <= time_max: |
---|
1270 | w = f(t, point_id = k)[0] |
---|
1271 | z = f(t, point_id = k)[1] |
---|
1272 | uh = f(t, point_id = k)[2] |
---|
1273 | vh = f(t, point_id = k)[3] |
---|
1274 | depth = w-z |
---|
1275 | m = sqrt(uh*uh + vh*vh) |
---|
1276 | if depth < 0.001: |
---|
1277 | vel = 0.0 |
---|
1278 | else: |
---|
1279 | vel = m / (depth + 1.e-6/depth) |
---|
1280 | bearing = calc_bearing(uh, vh) |
---|
1281 | model_time[i,k,j] = (t + starttime)/scale #t/60.0 |
---|
1282 | stages[i,k,j] = w |
---|
1283 | elevations[i,k,j] = z |
---|
1284 | xmom[i,k,j] = uh |
---|
1285 | ymom[i,k,j] = vh |
---|
1286 | momenta[i,k,j] = m |
---|
1287 | speed[i,k,j] = vel |
---|
1288 | bearings[i,k,j] = bearing |
---|
1289 | depths[i,k,j] = depth |
---|
1290 | thisgauge = gauges[k] |
---|
1291 | eastings[i,k,j] = thisgauge[0] |
---|
1292 | s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \ |
---|
1293 | % (t, w, m, vel, z, uh, vh, bearing) |
---|
1294 | fid_out.write(s) |
---|
1295 | if t == 0: |
---|
1296 | s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w) |
---|
1297 | fid_0.write(s) |
---|
1298 | if t/60.0 <= 13920: tindex = i |
---|
1299 | if w > max_stage: max_stage = w |
---|
1300 | if w < min_stage: min_stage = w |
---|
1301 | if m > max_momentum: max_momentum = m |
---|
1302 | if m < min_momentum: min_momentum = m |
---|
1303 | if uh > max_xmomentum: max_xmomentum = uh |
---|
1304 | if vh > max_ymomentum: max_ymomentum = vh |
---|
1305 | if uh < min_xmomentum: min_xmomentum = uh |
---|
1306 | if vh < min_ymomentum: min_ymomentum = vh |
---|
1307 | if vel > max_speed: max_speed = vel |
---|
1308 | if vel < min_speed: min_speed = vel |
---|
1309 | if z > 0 and depth > max_depth: max_depth = depth |
---|
1310 | |
---|
1311 | |
---|
1312 | s = '%.2f, %.2f, %.2f, %.2f, %s\n' \ |
---|
1313 | % (max_stage, min_stage, z, thisgauge[0], leg_label[j]) |
---|
1314 | fid_compare.write(s) |
---|
1315 | max_stages.append(max_stage) |
---|
1316 | min_stages.append(min_stage) |
---|
1317 | max_momentums.append(max_momentum) |
---|
1318 | max_xmomentums.append(max_xmomentum) |
---|
1319 | max_ymomentums.append(max_ymomentum) |
---|
1320 | min_xmomentums.append(min_xmomentum) |
---|
1321 | min_ymomentums.append(min_ymomentum) |
---|
1322 | min_momentums.append(min_momentum) |
---|
1323 | max_depths.append(max_depth) |
---|
1324 | max_speeds.append(max_speed) |
---|
1325 | min_speeds.append(min_speed) |
---|
1326 | #### finished generating quantities for each swwfile ##### |
---|
1327 | |
---|
1328 | model_time_plot3d[:,:] = model_time[:,:,j] |
---|
1329 | stages_plot3d[:,:] = stages[:,:,j] |
---|
1330 | eastings_plot3d[:,] = eastings[:,:,j] |
---|
1331 | |
---|
1332 | if surface is True: |
---|
1333 | print 'Printing surface figure' |
---|
1334 | for i in range(2): |
---|
1335 | fig = p1.figure(10) |
---|
1336 | ax = p3.Axes3D(fig) |
---|
1337 | if len(gauges) > 80: |
---|
1338 | ax.plot_surface(model_time[:,:,j], |
---|
1339 | eastings[:,:,j], |
---|
1340 | stages[:,:,j]) |
---|
1341 | else: |
---|
1342 | ax.plot3D(ravel(eastings[:,:,j]), |
---|
1343 | ravel(model_time[:,:,j]), |
---|
1344 | ravel(stages[:,:,j])) |
---|
1345 | ax.set_xlabel('time') |
---|
1346 | ax.set_ylabel('x') |
---|
1347 | ax.set_zlabel('stage') |
---|
1348 | fig.add_axes(ax) |
---|
1349 | p1.show() |
---|
1350 | surfacefig = 'solution_surface%s' % leg_label[j] |
---|
1351 | p1.savefig(surfacefig) |
---|
1352 | p1.close() |
---|
1353 | |
---|
1354 | #### finished generating quantities for all swwfiles ##### |
---|
1355 | |
---|
1356 | # x profile for given time |
---|
1357 | if surface is True: |
---|
1358 | figure(11) |
---|
1359 | plot(eastings[tindex,:,j], stages[tindex,:,j]) |
---|
1360 | xlabel('x') |
---|
1361 | ylabel('stage') |
---|
1362 | profilefig = 'solution_xprofile' |
---|
1363 | savefig('profilefig') |
---|
1364 | |
---|
1365 | elev_output = [] |
---|
1366 | if generate_fig is True: |
---|
1367 | depth_axis = axis([starttime/scale, time_max/scale, -0.1, |
---|
1368 | max(max_depths)*1.1]) |
---|
1369 | stage_axis = axis([starttime/scale, time_max/scale, |
---|
1370 | min(min_stages), max(max_stages)*1.1]) |
---|
1371 | vel_axis = axis([starttime/scale, time_max/scale, |
---|
1372 | min(min_speeds), max(max_speeds)*1.1]) |
---|
1373 | mom_axis = axis([starttime/scale, time_max/scale, |
---|
1374 | min(min_momentums), max(max_momentums)*1.1]) |
---|
1375 | xmom_axis = axis([starttime/scale, time_max/scale, |
---|
1376 | min(min_xmomentums), max(max_xmomentums)*1.1]) |
---|
1377 | ymom_axis = axis([starttime/scale, time_max/scale, |
---|
1378 | min(min_ymomentums), max(max_ymomentums)*1.1]) |
---|
1379 | cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] |
---|
1380 | nn = len(plot_quantity) |
---|
1381 | no_cols = 2 |
---|
1382 | |
---|
1383 | if len(label_id) > 1: graphname_report = [] |
---|
1384 | pp = 1 |
---|
1385 | div = 11. |
---|
1386 | cc = 0 |
---|
1387 | for k in gauge_index: |
---|
1388 | g = gauges[k] |
---|
1389 | count1 = 0 |
---|
1390 | if report == True and len(label_id) > 1: |
---|
1391 | s = '\\begin{figure}[ht] \n' \ |
---|
1392 | '\\centering \n' \ |
---|
1393 | '\\begin{tabular}{cc} \n' |
---|
1394 | fid.write(s) |
---|
1395 | if len(label_id) > 1: graphname_report = [] |
---|
1396 | |
---|
1397 | #### generate figures for each gauge #### |
---|
1398 | for j, f in enumerate(f_list): |
---|
1399 | ion() |
---|
1400 | hold(True) |
---|
1401 | count = 0 |
---|
1402 | where1 = 0 |
---|
1403 | where2 = 0 |
---|
1404 | word_quantity = '' |
---|
1405 | if report == True and len(label_id) == 1: |
---|
1406 | s = '\\begin{figure}[hbt] \n' \ |
---|
1407 | '\\centering \n' \ |
---|
1408 | '\\begin{tabular}{cc} \n' |
---|
1409 | fid.write(s) |
---|
1410 | |
---|
1411 | for which_quantity in plot_quantity: |
---|
1412 | count += 1 |
---|
1413 | where1 += 1 |
---|
1414 | figure(count, frameon = False) |
---|
1415 | if which_quantity == 'depth': |
---|
1416 | plot(model_time[0:n[j]-1,k,j], |
---|
1417 | depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1418 | units = 'm' |
---|
1419 | axis(depth_axis) |
---|
1420 | if which_quantity == 'stage': |
---|
1421 | if elevations[0,k,j] <= 0: |
---|
1422 | plot(model_time[0:n[j]-1,k,j], |
---|
1423 | stages[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1424 | axis(stage_axis) |
---|
1425 | else: |
---|
1426 | plot(model_time[0:n[j]-1,k,j], |
---|
1427 | depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1428 | #axis(depth_axis) |
---|
1429 | units = 'm' |
---|
1430 | if which_quantity == 'momentum': |
---|
1431 | plot(model_time[0:n[j]-1,k,j], |
---|
1432 | momenta[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1433 | axis(mom_axis) |
---|
1434 | units = 'm^2 / sec' |
---|
1435 | if which_quantity == 'xmomentum': |
---|
1436 | plot(model_time[0:n[j]-1,k,j], |
---|
1437 | xmom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1438 | axis(xmom_axis) |
---|
1439 | units = 'm^2 / sec' |
---|
1440 | if which_quantity == 'ymomentum': |
---|
1441 | plot(model_time[0:n[j]-1,k,j], |
---|
1442 | ymom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1443 | axis(ymom_axis) |
---|
1444 | units = 'm^2 / sec' |
---|
1445 | if which_quantity == 'speed': |
---|
1446 | plot(model_time[0:n[j]-1,k,j], |
---|
1447 | speed[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1448 | axis(vel_axis) |
---|
1449 | units = 'm / sec' |
---|
1450 | if which_quantity == 'bearing': |
---|
1451 | plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-', |
---|
1452 | model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', |
---|
1453 | model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.') |
---|
1454 | units = 'degrees from North' |
---|
1455 | #ax = axis([time_min, time_max, 0.0, 360.0]) |
---|
1456 | legend(('Bearing','West','East')) |
---|
1457 | |
---|
1458 | if time_unit is 'mins': xlabel('time (mins)') |
---|
1459 | if time_unit is 'hours': xlabel('time (hours)') |
---|
1460 | #if which_quantity == 'stage' \ |
---|
1461 | # and elevations[0:n[j]-1,k,j] > 0: |
---|
1462 | # ylabel('%s (%s)' %('depth', units)) |
---|
1463 | #else: |
---|
1464 | # ylabel('%s (%s)' %(which_quantity, units)) |
---|
1465 | #ylabel('%s (%s)' %('wave height', units)) |
---|
1466 | ylabel('%s (%s)' %(which_quantity, units)) |
---|
1467 | if len(label_id) > 1: legend((leg_label),loc='upper right') |
---|
1468 | |
---|
1469 | #gaugeloc1 = gaugeloc.replace(' ','') |
---|
1470 | #gaugeloc2 = gaugeloc1.replace('_','') |
---|
1471 | gaugeloc2 = str(locations[k]).replace(' ','') |
---|
1472 | graphname = '%sgauge%s_%s' %(file_loc[j], |
---|
1473 | gaugeloc2, |
---|
1474 | which_quantity) |
---|
1475 | |
---|
1476 | if report == True and len(label_id) > 1: |
---|
1477 | figdir = getcwd()+sep+'report_figures'+sep |
---|
1478 | if access(figdir,F_OK) == 0 : |
---|
1479 | mkdir (figdir) |
---|
1480 | latex_file_loc = figdir.replace(sep,altsep) |
---|
1481 | # storing files in production directory |
---|
1482 | graphname_latex = '%sgauge%s%s' \ |
---|
1483 | % (latex_file_loc, gaugeloc2, |
---|
1484 | which_quantity) |
---|
1485 | # giving location in latex output file |
---|
1486 | graphname_report_input = '%sgauge%s%s' % \ |
---|
1487 | ('..' + altsep + |
---|
1488 | 'report_figures' + altsep, |
---|
1489 | gaugeloc2, which_quantity) |
---|
1490 | graphname_report.append(graphname_report_input) |
---|
1491 | |
---|
1492 | # save figures in production directory for report |
---|
1493 | savefig(graphname_latex) |
---|
1494 | |
---|
1495 | if report == True: |
---|
1496 | figdir = getcwd() + sep + 'report_figures' + sep |
---|
1497 | if access(figdir,F_OK) == 0: |
---|
1498 | mkdir(figdir) |
---|
1499 | latex_file_loc = figdir.replace(sep,altsep) |
---|
1500 | |
---|
1501 | if len(label_id) == 1: |
---|
1502 | # storing files in production directory |
---|
1503 | graphname_latex = '%sgauge%s%s%s' % \ |
---|
1504 | (latex_file_loc, gaugeloc2, |
---|
1505 | which_quantity, label_id2) |
---|
1506 | # giving location in latex output file |
---|
1507 | graphname_report = '%sgauge%s%s%s' % \ |
---|
1508 | ('..' + altsep + |
---|
1509 | 'report_figures' + altsep, |
---|
1510 | gaugeloc2, which_quantity, |
---|
1511 | label_id2) |
---|
1512 | s = '\includegraphics' \ |
---|
1513 | '[width=0.49\linewidth, height=50mm]{%s%s}' % \ |
---|
1514 | (graphname_report, '.png') |
---|
1515 | fid.write(s) |
---|
1516 | if where1 % 2 == 0: |
---|
1517 | s = '\\\\ \n' |
---|
1518 | where1 = 0 |
---|
1519 | else: |
---|
1520 | s = '& \n' |
---|
1521 | fid.write(s) |
---|
1522 | savefig(graphname_latex) |
---|
1523 | |
---|
1524 | if title_on == True: |
---|
1525 | title('%s scenario: %s at %s gauge' % \ |
---|
1526 | (label_id, which_quantity, gaugeloc2)) |
---|
1527 | #title('Gauge %s (MOST elevation %.2f, ' \ |
---|
1528 | # 'ANUGA elevation %.2f)' % \ |
---|
1529 | # (gaugeloc2, elevations[10,k,0], |
---|
1530 | # elevations[10,k,1])) |
---|
1531 | |
---|
1532 | savefig(graphname) # save figures with sww file |
---|
1533 | |
---|
1534 | if report == True and len(label_id) == 1: |
---|
1535 | for i in range(nn-1): |
---|
1536 | if nn > 2: |
---|
1537 | if plot_quantity[i] == 'stage' \ |
---|
1538 | and elevations[0,k,j] > 0: |
---|
1539 | word_quantity += 'depth' + ', ' |
---|
1540 | else: |
---|
1541 | word_quantity += plot_quantity[i] + ', ' |
---|
1542 | else: |
---|
1543 | if plot_quantity[i] == 'stage' \ |
---|
1544 | and elevations[0,k,j] > 0: |
---|
1545 | word_quantity += 'depth' + ', ' |
---|
1546 | else: |
---|
1547 | word_quantity += plot_quantity[i] |
---|
1548 | |
---|
1549 | if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0: |
---|
1550 | word_quantity += ' and ' + 'depth' |
---|
1551 | else: |
---|
1552 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
1553 | caption = 'Time series for %s at %s location ' \ |
---|
1554 | '(elevation %.2fm)' % \ |
---|
1555 | (word_quantity, locations[k], elev[k]) |
---|
1556 | if elev[k] == 0.0: |
---|
1557 | caption = 'Time series for %s at %s location ' \ |
---|
1558 | '(elevation %.2fm)' % \ |
---|
1559 | (word_quantity, locations[k], |
---|
1560 | elevations[0,k,j]) |
---|
1561 | east = gauges[0] |
---|
1562 | north = gauges[1] |
---|
1563 | elev_output.append([locations[k], east, north, |
---|
1564 | elevations[0,k,j]]) |
---|
1565 | label = '%sgauge%s' % (label_id2, gaugeloc2) |
---|
1566 | s = '\end{tabular} \n' \ |
---|
1567 | '\\caption{%s} \n' \ |
---|
1568 | '\label{fig:%s} \n' \ |
---|
1569 | '\end{figure} \n \n' % (caption, label) |
---|
1570 | fid.write(s) |
---|
1571 | cc += 1 |
---|
1572 | if cc % 6 == 0: fid.write('\\clearpage \n') |
---|
1573 | savefig(graphname_latex) |
---|
1574 | |
---|
1575 | if report == True and len(label_id) > 1: |
---|
1576 | for i in range(nn-1): |
---|
1577 | if nn > 2: |
---|
1578 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1579 | word_quantity += 'depth' + ',' |
---|
1580 | else: |
---|
1581 | word_quantity += plot_quantity[i] + ', ' |
---|
1582 | else: |
---|
1583 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1584 | word_quantity += 'depth' |
---|
1585 | else: |
---|
1586 | word_quantity += plot_quantity[i] |
---|
1587 | where1 = 0 |
---|
1588 | count1 += 1 |
---|
1589 | index = j*len(plot_quantity) |
---|
1590 | for which_quantity in plot_quantity: |
---|
1591 | where1 += 1 |
---|
1592 | s = '\includegraphics' \ |
---|
1593 | '[width=0.49\linewidth, height=50mm]{%s%s}' % \ |
---|
1594 | (graphname_report[index], '.png') |
---|
1595 | index += 1 |
---|
1596 | fid.write(s) |
---|
1597 | if where1 % 2 == 0: |
---|
1598 | s = '\\\\ \n' |
---|
1599 | where1 = 0 |
---|
1600 | else: |
---|
1601 | s = '& \n' |
---|
1602 | fid.write(s) |
---|
1603 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
1604 | label = 'gauge%s' %(gaugeloc2) |
---|
1605 | caption = 'Time series for %s at %s location ' \ |
---|
1606 | '(elevation %.2fm)' % \ |
---|
1607 | (word_quantity, locations[k], elev[k]) |
---|
1608 | if elev[k] == 0.0: |
---|
1609 | caption = 'Time series for %s at %s location ' \ |
---|
1610 | '(elevation %.2fm)' % \ |
---|
1611 | (word_quantity, locations[k], |
---|
1612 | elevations[0,k,j]) |
---|
1613 | thisgauge = gauges[k] |
---|
1614 | east = thisgauge[0] |
---|
1615 | north = thisgauge[1] |
---|
1616 | elev_output.append([locations[k], east, north, |
---|
1617 | elevations[0,k,j]]) |
---|
1618 | |
---|
1619 | s = '\end{tabular} \n' \ |
---|
1620 | '\\caption{%s} \n' \ |
---|
1621 | '\label{fig:%s} \n' \ |
---|
1622 | '\end{figure} \n \n' % (caption, label) |
---|
1623 | fid.write(s) |
---|
1624 | if float((k+1)/div - pp) == 0.: |
---|
1625 | fid.write('\\clearpage \n') |
---|
1626 | pp += 1 |
---|
1627 | #### finished generating figures ### |
---|
1628 | |
---|
1629 | close('all') |
---|
1630 | |
---|
1631 | return texfile2, elev_output |
---|
1632 | |
---|
1633 | |
---|
1634 | # FIXME (DSG): Add unit test, make general, not just 2 files, |
---|
1635 | # but any number of files. |
---|
1636 | ## |
---|
1637 | # @brief ?? |
---|
1638 | # @param dir_name ?? |
---|
1639 | # @param filename1 ?? |
---|
1640 | # @param filename2 ?? |
---|
1641 | # @return ?? |
---|
1642 | # @note TEMP |
---|
1643 | def copy_code_files(dir_name, filename1, filename2): |
---|
1644 | """Temporary Interface to new location""" |
---|
1645 | |
---|
1646 | from anuga.shallow_water.data_manager import \ |
---|
1647 | copy_code_files as dm_copy_code_files |
---|
1648 | print 'copy_code_files has moved from util.py. ', |
---|
1649 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1650 | import copy_code_files"' |
---|
1651 | |
---|
1652 | return dm_copy_code_files(dir_name, filename1, filename2) |
---|
1653 | |
---|
1654 | |
---|
1655 | ## |
---|
1656 | # @brief Create a nested sub-directory path. |
---|
1657 | # @param root_directory The base diretory path. |
---|
1658 | # @param directories An iterable of sub-directory names. |
---|
1659 | # @return The final joined directory path. |
---|
1660 | # @note If each sub-directory doesn't exist, it will be created. |
---|
1661 | def add_directories(root_directory, directories): |
---|
1662 | """ |
---|
1663 | Add the first sub-directory in 'directories' to root_directory. |
---|
1664 | Then add the second sub-directory to the accumulating path and so on. |
---|
1665 | |
---|
1666 | Return the path of the final directory. |
---|
1667 | |
---|
1668 | This is handy for specifying and creating a directory where data will go. |
---|
1669 | """ |
---|
1670 | dir = root_directory |
---|
1671 | for new_dir in directories: |
---|
1672 | dir = os.path.join(dir, new_dir) |
---|
1673 | if not access(dir,F_OK): |
---|
1674 | mkdir(dir) |
---|
1675 | return dir |
---|
1676 | |
---|
1677 | |
---|
1678 | ## |
---|
1679 | # @brief |
---|
1680 | # @param filename |
---|
1681 | # @param separator_value |
---|
1682 | # @return |
---|
1683 | # @note TEMP |
---|
1684 | def get_data_from_file(filename, separator_value=','): |
---|
1685 | """Temporary Interface to new location""" |
---|
1686 | from anuga.shallow_water.data_manager import \ |
---|
1687 | get_data_from_file as dm_get_data_from_file |
---|
1688 | print 'get_data_from_file has moved from util.py' |
---|
1689 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1690 | import get_data_from_file"' |
---|
1691 | |
---|
1692 | return dm_get_data_from_file(filename,separator_value = ',') |
---|
1693 | |
---|
1694 | |
---|
1695 | ## |
---|
1696 | # @brief |
---|
1697 | # @param verbose |
---|
1698 | # @param kwargs |
---|
1699 | # @return |
---|
1700 | # @note TEMP |
---|
1701 | def store_parameters(verbose=False,**kwargs): |
---|
1702 | """Temporary Interface to new location""" |
---|
1703 | |
---|
1704 | from anuga.shallow_water.data_manager \ |
---|
1705 | import store_parameters as dm_store_parameters |
---|
1706 | print 'store_parameters has moved from util.py.' |
---|
1707 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1708 | import store_parameters"' |
---|
1709 | |
---|
1710 | return dm_store_parameters(verbose=False,**kwargs) |
---|
1711 | |
---|
1712 | |
---|
1713 | ## |
---|
1714 | # @brief Remove vertices that are not associated with any triangle. |
---|
1715 | # @param verts An iterable (or array) of points. |
---|
1716 | # @param triangles An iterable of 3 element tuples. |
---|
1717 | # @param number_of_full_nodes ?? |
---|
1718 | # @return (verts, triangles) where 'verts' has been updated. |
---|
1719 | def remove_lone_verts(verts, triangles, number_of_full_nodes=None): |
---|
1720 | """Removes vertices that are not associated with any triangles. |
---|
1721 | |
---|
1722 | verts is a list/array of points. |
---|
1723 | triangles is a list of 3 element tuples. Each tuple represents a triangle. |
---|
1724 | number_of_full_nodes relate to parallelism when a mesh has an |
---|
1725 | extra layer of ghost points. |
---|
1726 | """ |
---|
1727 | |
---|
1728 | verts = ensure_numeric(verts) |
---|
1729 | triangles = ensure_numeric(triangles) |
---|
1730 | |
---|
1731 | N = len(verts) |
---|
1732 | |
---|
1733 | # initialise the array to easily find the index of the first loner |
---|
1734 | # ie, if N=3 -> [6,5,4] |
---|
1735 | loners=arange(2*N, N, -1) |
---|
1736 | for t in triangles: |
---|
1737 | for vert in t: |
---|
1738 | loners[vert]= vert # all non-loners will have loners[i]=i |
---|
1739 | |
---|
1740 | lone_start = 2*N - max(loners) # The index of the first loner |
---|
1741 | |
---|
1742 | if lone_start-1 == N: |
---|
1743 | # no loners |
---|
1744 | pass |
---|
1745 | elif min(loners[lone_start:N]) > N: |
---|
1746 | # All the loners are at the end of the vert array |
---|
1747 | verts = verts[0:lone_start] |
---|
1748 | else: |
---|
1749 | # change the loners list so it can be used to modify triangles |
---|
1750 | # Remove the loners from verts |
---|
1751 | # Could've used X=compress(less(loners,N),loners) |
---|
1752 | # verts=take(verts,X) to Remove the loners from verts |
---|
1753 | # but I think it would use more memory |
---|
1754 | new_i = lone_start # point at first loner - first 'shuffle down' target |
---|
1755 | for i in range(lone_start, N): |
---|
1756 | if loners[i] >= N: # [i] is a loner, leave alone |
---|
1757 | pass |
---|
1758 | else: # a non-loner, move down |
---|
1759 | loners[i] = new_i |
---|
1760 | verts[new_i] = verts[i] |
---|
1761 | new_i += 1 |
---|
1762 | verts = verts[0:new_i] |
---|
1763 | |
---|
1764 | # Modify the triangles |
---|
1765 | #print "loners", loners |
---|
1766 | #print "triangles before", triangles |
---|
1767 | triangles = choose(triangles,loners) |
---|
1768 | #print "triangles after", triangles |
---|
1769 | return verts, triangles |
---|
1770 | |
---|
1771 | |
---|
1772 | def get_centroid_values(x, triangles): |
---|
1773 | """Compute centroid values from vertex values |
---|
1774 | |
---|
1775 | x: Values at vertices of triangular mesh |
---|
1776 | triangles: Nx3 integer array pointing to vertex information |
---|
1777 | for each of the N triangels. Elements of triangles are |
---|
1778 | indices into x |
---|
1779 | """ |
---|
1780 | |
---|
1781 | |
---|
1782 | xc = zeros(triangles.shape[0], Float) # Space for centroid info |
---|
1783 | |
---|
1784 | for k in range(triangles.shape[0]): |
---|
1785 | # Indices of vertices |
---|
1786 | i0 = triangles[k][0] |
---|
1787 | i1 = triangles[k][1] |
---|
1788 | i2 = triangles[k][2] |
---|
1789 | |
---|
1790 | xc[k] = (x[i0] + x[i1] + x[i2])/3 |
---|
1791 | |
---|
1792 | |
---|
1793 | return xc |
---|
1794 | |
---|
1795 | # @note TEMP |
---|
1796 | def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]}, |
---|
1797 | output_dir='', |
---|
1798 | base_name='', |
---|
1799 | plot_numbers=['3-5'], |
---|
1800 | quantities=['speed','stage','momentum'], |
---|
1801 | assess_all_csv_files=True, |
---|
1802 | extra_plot_name='test' ): |
---|
1803 | |
---|
1804 | msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs ', |
---|
1805 | msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import csv2timeseries_graphs"' |
---|
1806 | |
---|
1807 | raise Exception, msg |
---|
1808 | return csv2timeseries_graphs(directories_dic, |
---|
1809 | output_dir, |
---|
1810 | base_name, |
---|
1811 | plot_numbers, |
---|
1812 | quantities, |
---|
1813 | extra_plot_name, |
---|
1814 | assess_all_csv_files |
---|
1815 | ) |
---|
1816 | |
---|
1817 | def csv2timeseries_graphs(directories_dic={}, |
---|
1818 | output_dir='', |
---|
1819 | base_name=None, |
---|
1820 | plot_numbers='', |
---|
1821 | quantities=['stage'], |
---|
1822 | extra_plot_name='', |
---|
1823 | assess_all_csv_files=True, |
---|
1824 | create_latex=False, |
---|
1825 | verbose=False): |
---|
1826 | |
---|
1827 | """ |
---|
1828 | Read in csv files that have the right header information and |
---|
1829 | plot time series such as Stage, Speed, etc. Will also plot several |
---|
1830 | time series on one plot. Filenames must follow this convention, |
---|
1831 | <base_name><plot_number>.csv eg gauge_timeseries3.csv |
---|
1832 | |
---|
1833 | NOTE: relies that 'elevation' is in the csv file! |
---|
1834 | |
---|
1835 | Each file represents a location and within each file there are |
---|
1836 | time, quantity columns. |
---|
1837 | |
---|
1838 | For example: |
---|
1839 | if "directories_dic" defines 4 directories and in each directories |
---|
1840 | there is a csv files corresponding to the right "plot_numbers", |
---|
1841 | this will create a plot with 4 lines one for each directory AND |
---|
1842 | one plot for each "quantities". ??? FIXME: unclear. |
---|
1843 | |
---|
1844 | Usage: |
---|
1845 | csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0], |
---|
1846 | 'fixed_wave'+sep:['Fixed Wave',0,0]}, |
---|
1847 | output_dir='fixed_wave'+sep, |
---|
1848 | base_name='gauge_timeseries_', |
---|
1849 | plot_numbers='', |
---|
1850 | quantities=['stage','speed'], |
---|
1851 | extra_plot_name='', |
---|
1852 | assess_all_csv_files=True, |
---|
1853 | create_latex=False, |
---|
1854 | verbose=True) |
---|
1855 | this will create one plot for stage with both 'slide' and |
---|
1856 | 'fixed_wave' lines on it for stage and speed for each csv |
---|
1857 | file with 'gauge_timeseries_' as the prefix. The graghs |
---|
1858 | will be in the output directory 'fixed_wave' and the graph |
---|
1859 | axis will be determined by assessing all the |
---|
1860 | |
---|
1861 | ANOTHER EXAMPLE |
---|
1862 | new_csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0], |
---|
1863 | 'fixed_wave'+sep:['Fixed Wave',0,0]}, |
---|
1864 | output_dir='fixed_wave'+sep, |
---|
1865 | base_name='gauge_timeseries_', |
---|
1866 | plot_numbers=['1-3'], |
---|
1867 | quantities=['stage','speed'], |
---|
1868 | extra_plot_name='', |
---|
1869 | assess_all_csv_files=False, |
---|
1870 | create_latex=False, |
---|
1871 | verbose=True) |
---|
1872 | This will plot csv files called gauge_timeseries_1.csv and |
---|
1873 | gauge_timeseries3.csv from both 'slide' and 'fixed_wave' directories |
---|
1874 | to 'fixed_wave'. There will be 4 plots created two speed and two stage |
---|
1875 | one for each csv file. There will be two lines on each of these plots. |
---|
1876 | And the axis will have been determined from only these files, had |
---|
1877 | assess_all_csv_files = True all csv file with 'gauges_timeseries_' prefix |
---|
1878 | would of been assessed. |
---|
1879 | |
---|
1880 | ANOTHER EXAMPLE |
---|
1881 | csv2timeseries_graphs({'J:'+sep+'anuga_validation'+sep:['new',20,-.1], |
---|
1882 | 'J:'+sep+'conical_island'+sep:['test',0,0]}, |
---|
1883 | output_dir='', |
---|
1884 | plot_numbers=['1','3'], |
---|
1885 | quantities=['stage','depth','bearing'], |
---|
1886 | base_name='gauge_b', |
---|
1887 | assess_all_csv_files=True, |
---|
1888 | verbose=True) |
---|
1889 | |
---|
1890 | This will produce one plot for each quantity (therefore 3) in the current directory, |
---|
1891 | each plot will have 2 lines on them. The first plot named 'new' will have the time |
---|
1892 | offseted by 20secs and the stage height adjusted by -0.1m |
---|
1893 | |
---|
1894 | Inputs: |
---|
1895 | directories_dic: dictionary of directory with values (plot |
---|
1896 | legend name for directory), (start time of |
---|
1897 | the time series) and the (value to add to |
---|
1898 | stage if needed). For example |
---|
1899 | {dir1:['Anuga_ons',5000, 0], |
---|
1900 | dir2:['b_emoth',5000,1.5], |
---|
1901 | dir3:['b_ons',5000,1.5]} |
---|
1902 | Having multiple directories defined will plot them on |
---|
1903 | one plot, therefore there will be 3 lines on each of these |
---|
1904 | plot. If you only want one line per plot call csv2timeseries_graph |
---|
1905 | separately for each directory, eg only have one directory in the |
---|
1906 | 'directories_dic' in each call. |
---|
1907 | |
---|
1908 | output_dir: directory for the plot outputs. Only important to define |
---|
1909 | when you have more than one directory in your directories_dic, |
---|
1910 | if you have not defined it and you have multiple directories in |
---|
1911 | 'directories_dic' there will be plots in each directory however only |
---|
1912 | one directory will contain the complete plot/graphs. |
---|
1913 | |
---|
1914 | base_name: Is used a couple of times. 1) to find the csv files to be plotted |
---|
1915 | if there is no 'plot_numbers' then csv files with 'base_name' are |
---|
1916 | plotted and 2) in the title of the plots, the lenght of base_name is |
---|
1917 | removed from the front of the filename to be used in the title. |
---|
1918 | This could be changed if needed. |
---|
1919 | Note is ignored if assess_all_csv_files=True |
---|
1920 | |
---|
1921 | plot_numbers: a String list of numbers to plot. For example |
---|
1922 | [0-4,10,15-17] will read and attempt to plot |
---|
1923 | the follow 0,1,2,3,4,10,15,16,17 |
---|
1924 | NOTE: if no plot numbers this will create |
---|
1925 | one plot per quantity, per gauge |
---|
1926 | quantities: Will get available quantities from the header in the csv file. |
---|
1927 | quantities must be one of these. |
---|
1928 | NOTE: ALL QUANTITY NAMES MUST BE lower case! |
---|
1929 | |
---|
1930 | extra_plot_name: A string that is appended to the end of the |
---|
1931 | output filename. |
---|
1932 | |
---|
1933 | assess_all_csv_files: if true it will read ALL csv file with |
---|
1934 | "base_name", regardless of 'plot_numbers' |
---|
1935 | and determine a uniform set of axes for |
---|
1936 | Stage, Speed and Momentum. IF FALSE it |
---|
1937 | will only read the csv file within the |
---|
1938 | 'plot_numbers' |
---|
1939 | |
---|
1940 | create_latex: NOT IMPLEMENTED YET!! sorry Jane.... |
---|
1941 | |
---|
1942 | OUTPUTS: saves the plots to |
---|
1943 | <output_dir><base_name><plot_number><extra_plot_name>.png |
---|
1944 | |
---|
1945 | |
---|
1946 | """ |
---|
1947 | try: |
---|
1948 | import pylab |
---|
1949 | except ImportError: |
---|
1950 | msg='csv2timeseries_graphs needs pylab to be installed correctly' |
---|
1951 | raise msg |
---|
1952 | #ANUGA don't need pylab to work so the system doesn't |
---|
1953 | #rely on pylab being installed |
---|
1954 | return |
---|
1955 | |
---|
1956 | from os import sep |
---|
1957 | from anuga.shallow_water.data_manager import \ |
---|
1958 | get_all_files_with_extension, csv2dict |
---|
1959 | |
---|
1960 | seconds_in_hour = 3600 |
---|
1961 | seconds_in_minutes = 60 |
---|
1962 | |
---|
1963 | quantities_label={} |
---|
1964 | # quantities_label['time'] = 'time (hours)' |
---|
1965 | quantities_label['time'] = 'time (minutes)' |
---|
1966 | quantities_label['stage'] = 'wave height (m)' |
---|
1967 | quantities_label['speed'] = 'speed (m/s)' |
---|
1968 | quantities_label['momentum'] = 'momentum (m^2/sec)' |
---|
1969 | quantities_label['depth'] = 'water depth (m)' |
---|
1970 | quantities_label['xmomentum'] = 'momentum (m^2/sec)' |
---|
1971 | quantities_label['ymomentum'] = 'momentum (m^2/sec)' |
---|
1972 | quantities_label['bearing'] = 'degrees (o)' |
---|
1973 | quantities_label['elevation'] = 'elevation (m)' |
---|
1974 | |
---|
1975 | |
---|
1976 | if extra_plot_name != '': |
---|
1977 | extra_plot_name='_'+extra_plot_name |
---|
1978 | |
---|
1979 | new_plot_numbers=[] |
---|
1980 | #change plot_numbers to list, eg ['0-4','10'] |
---|
1981 | #to ['0','1','2','3','4','10'] |
---|
1982 | for i, num_string in enumerate(plot_numbers): |
---|
1983 | if '-' in num_string: |
---|
1984 | start = int(num_string[:num_string.rfind('-')]) |
---|
1985 | end = int(num_string[num_string.rfind('-')+1:])+1 |
---|
1986 | for x in range(start, end): |
---|
1987 | new_plot_numbers.append(str(x)) |
---|
1988 | else: |
---|
1989 | new_plot_numbers.append(num_string) |
---|
1990 | |
---|
1991 | #finds all the files that fit the specs provided and return a list of them |
---|
1992 | #so to help find a uniform max and min for the plots... |
---|
1993 | list_filenames=[] |
---|
1994 | all_csv_filenames=[] |
---|
1995 | if verbose: print 'Determining files to access for axes ranges \n' |
---|
1996 | |
---|
1997 | #print directories_dic.keys(), base_name |
---|
1998 | |
---|
1999 | for i,directory in enumerate(directories_dic.keys()): |
---|
2000 | all_csv_filenames.append(get_all_files_with_extension(directory, |
---|
2001 | base_name,'.csv')) |
---|
2002 | |
---|
2003 | filenames=[] |
---|
2004 | if plot_numbers == '': |
---|
2005 | list_filenames.append(get_all_files_with_extension(directory, |
---|
2006 | base_name,'.csv')) |
---|
2007 | else: |
---|
2008 | for number in new_plot_numbers: |
---|
2009 | # print 'number!!!', base_name, number |
---|
2010 | filenames.append(base_name+number) |
---|
2011 | # print filenames |
---|
2012 | list_filenames.append(filenames) |
---|
2013 | |
---|
2014 | |
---|
2015 | |
---|
2016 | #print "list_filenames", list_filenames |
---|
2017 | |
---|
2018 | #use all the files to get the values for the plot axis |
---|
2019 | max_start_time= -1000. |
---|
2020 | min_start_time = 100000 |
---|
2021 | |
---|
2022 | |
---|
2023 | if verbose: print 'Determining uniform axes \n' |
---|
2024 | #this entire loop is to determine the min and max range for the |
---|
2025 | #axes of the plots |
---|
2026 | |
---|
2027 | # quantities.insert(0,'elevation') |
---|
2028 | quantities.insert(0,'time') |
---|
2029 | |
---|
2030 | directory_quantity_value={} |
---|
2031 | # quantity_value={} |
---|
2032 | min_quantity_value={} |
---|
2033 | max_quantity_value={} |
---|
2034 | |
---|
2035 | |
---|
2036 | for i, directory in enumerate(directories_dic.keys()): |
---|
2037 | filename_quantity_value={} |
---|
2038 | if assess_all_csv_files==False: |
---|
2039 | which_csv_to_assess = list_filenames[i] |
---|
2040 | else: |
---|
2041 | #gets list of filenames for directory "i" |
---|
2042 | which_csv_to_assess = all_csv_filenames[i] |
---|
2043 | # print'IN DIR', list_filenames[i] |
---|
2044 | |
---|
2045 | |
---|
2046 | |
---|
2047 | |
---|
2048 | for j, filename in enumerate(which_csv_to_assess): |
---|
2049 | quantity_value={} |
---|
2050 | |
---|
2051 | dir_filename=join(directory,filename) |
---|
2052 | attribute_dic, title_index_dic = csv2dict(dir_filename+ |
---|
2053 | '.csv') |
---|
2054 | directory_start_time = directories_dic[directory][1] |
---|
2055 | directory_add_tide = directories_dic[directory][2] |
---|
2056 | |
---|
2057 | if verbose: print 'reading: %s.csv' %dir_filename |
---|
2058 | # print 'keys',attribute_dic.keys() |
---|
2059 | #add time to get values |
---|
2060 | for k, quantity in enumerate(quantities): |
---|
2061 | quantity_value[quantity] = [float(x) for x in attribute_dic[quantity]] |
---|
2062 | |
---|
2063 | #add tide to stage if provided |
---|
2064 | if quantity == 'stage': |
---|
2065 | quantity_value[quantity]=array(quantity_value[quantity])+directory_add_tide |
---|
2066 | |
---|
2067 | #condition to find max and mins for all the plots |
---|
2068 | # populate the list with something when i=0 and j=0 and |
---|
2069 | # then compare to the other values to determine abs max and min |
---|
2070 | if i==0 and j==0: |
---|
2071 | |
---|
2072 | min_quantity_value[quantity], \ |
---|
2073 | max_quantity_value[quantity] = get_min_max_values(quantity_value[quantity]) |
---|
2074 | |
---|
2075 | if quantity != 'time': |
---|
2076 | min_quantity_value[quantity] = min_quantity_value[quantity] *1.1 |
---|
2077 | max_quantity_value[quantity] = max_quantity_value[quantity] *1.1 |
---|
2078 | |
---|
2079 | |
---|
2080 | # print '1 min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename |
---|
2081 | else: |
---|
2082 | # print 'min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename |
---|
2083 | min, max = get_min_max_values(quantity_value[quantity]) |
---|
2084 | # print "MIN",min, max |
---|
2085 | |
---|
2086 | #min and max are multipled by "1+increase_axis" to get axes that are slighty bigger |
---|
2087 | # than the max and mins so the plots look good. |
---|
2088 | |
---|
2089 | increase_axis = (max-min)*0.05 |
---|
2090 | # print quantity, "MIN MAX", max, min |
---|
2091 | if min<=min_quantity_value[quantity]: |
---|
2092 | if quantity == 'time': |
---|
2093 | min_quantity_value[quantity]=min |
---|
2094 | else: |
---|
2095 | if round(min,2) == 0.00: |
---|
2096 | min_quantity_value[quantity]=-increase_axis |
---|
2097 | # min_quantity_value[quantity]=-2. |
---|
2098 | #min_quantity_value[quantity]= -max_quantity_value[quantity]*increase_axis |
---|
2099 | else: |
---|
2100 | # min_quantity_value[quantity]=min*(1+increase_axis) |
---|
2101 | min_quantity_value[quantity]=min-increase_axis |
---|
2102 | # print quantity, min_quantity_value[quantity] |
---|
2103 | |
---|
2104 | if max>max_quantity_value[quantity]: |
---|
2105 | if quantity == 'time': |
---|
2106 | max_quantity_value[quantity]=max |
---|
2107 | else: |
---|
2108 | max_quantity_value[quantity]=max+increase_axis |
---|
2109 | # max_quantity_value[quantity]=max*(1+increase_axis) |
---|
2110 | # print quantity, max_quantity_value[quantity],increase_axis |
---|
2111 | |
---|
2112 | # print 'min,maj',quantity, min_quantity_value[quantity],max_quantity_value[quantity] |
---|
2113 | |
---|
2114 | |
---|
2115 | |
---|
2116 | |
---|
2117 | #set the time... ??? |
---|
2118 | if min_start_time > directory_start_time: |
---|
2119 | min_start_time = directory_start_time |
---|
2120 | if max_start_time < directory_start_time: |
---|
2121 | max_start_time = directory_start_time |
---|
2122 | #print 'start_time' , max_start_time, min_start_time |
---|
2123 | |
---|
2124 | filename_quantity_value[filename]=quantity_value |
---|
2125 | |
---|
2126 | directory_quantity_value[directory]=filename_quantity_value |
---|
2127 | |
---|
2128 | #final step to unifrom axis for the graphs |
---|
2129 | quantities_axis={} |
---|
2130 | |
---|
2131 | for i, quantity in enumerate(quantities): |
---|
2132 | quantities_axis[quantity] = (float(min_start_time)/float(seconds_in_minutes), |
---|
2133 | (float(max_quantity_value['time'])+float(max_start_time))\ |
---|
2134 | /float(seconds_in_minutes), |
---|
2135 | min_quantity_value[quantity], |
---|
2136 | max_quantity_value[quantity]) |
---|
2137 | if verbose and (quantity != 'time' and quantity != 'elevation'): |
---|
2138 | print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' %(quantity, |
---|
2139 | quantities_axis[quantity][0], |
---|
2140 | quantities_axis[quantity][1], |
---|
2141 | quantities_label['time'], |
---|
2142 | quantities_axis[quantity][2], |
---|
2143 | quantities_axis[quantity][3], |
---|
2144 | quantities_label[quantity]) |
---|
2145 | |
---|
2146 | #print quantities_axis[quantity] |
---|
2147 | |
---|
2148 | cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k'] |
---|
2149 | |
---|
2150 | if verbose: print 'Now start to plot \n' |
---|
2151 | |
---|
2152 | i_max = len(directories_dic.keys()) |
---|
2153 | legend_list_dic={} |
---|
2154 | legend_list =[] |
---|
2155 | for i, directory in enumerate(directories_dic.keys()): |
---|
2156 | |
---|
2157 | if verbose: print'Plotting in %s %s' %(directory, new_plot_numbers) |
---|
2158 | # print 'LIST',list_filenames |
---|
2159 | #FIXME THIS SORT IS VERY IMPORTANT, without it the assigned plot numbers may not work correctly |
---|
2160 | #there must be a better way |
---|
2161 | list_filenames[i].sort() |
---|
2162 | for j, filename in enumerate(list_filenames[i]): |
---|
2163 | # print'IN plot', list_filenames[i] |
---|
2164 | |
---|
2165 | if verbose: print'Starting %s' %filename |
---|
2166 | directory_name = directories_dic[directory][0] |
---|
2167 | directory_start_time = directories_dic[directory][1] |
---|
2168 | directory_add_tide = directories_dic[directory][2] |
---|
2169 | |
---|
2170 | #create an if about the start time and tide hieght |
---|
2171 | #if they don't exist |
---|
2172 | #print 'i %s,j %s, number %s, file %s' %(i,j,number,file) |
---|
2173 | attribute_dic, title_index_dic = csv2dict(directory+sep+filename+'.csv') |
---|
2174 | #get data from dict in to list |
---|
2175 | #do maths to list by changing to array |
---|
2176 | t=(array(directory_quantity_value[directory][filename]['time'])+directory_start_time)/seconds_in_minutes |
---|
2177 | |
---|
2178 | #finds the maximum elevation, used only as a test |
---|
2179 | # and as info in the graphs |
---|
2180 | max_ele=-100000 |
---|
2181 | min_ele=100000 |
---|
2182 | elevation = [float(x) for x in attribute_dic["elevation"]] |
---|
2183 | |
---|
2184 | min_ele, max_ele = get_min_max_values(elevation) |
---|
2185 | |
---|
2186 | if min_ele != max_ele: |
---|
2187 | print "Note! Elevation changes in %s" %dir_filename |
---|
2188 | |
---|
2189 | # creates a dictionary with keys that is the filename and attributes are a list of |
---|
2190 | # lists containing 'directory_name' and 'elevation'. This is used to make the contents |
---|
2191 | # for the legends in the graphs, this is the name of the model and the elevation. |
---|
2192 | # All in this great one liner from DG. If the key 'filename' doesn't exist it creates the |
---|
2193 | # entry if the entry exist it appends to the key. |
---|
2194 | |
---|
2195 | legend_list_dic.setdefault(filename,[]).append([directory_name,round(max_ele,3)]) |
---|
2196 | |
---|
2197 | # creates a LIST for the legend on the last iteration of the directories |
---|
2198 | # which is when "legend_list_dic" has been fully populated. Creates a list of strings |
---|
2199 | # which is used in the legend |
---|
2200 | # only runs on the last iteration for all the gauges(csv) files |
---|
2201 | # empties the list before creating it |
---|
2202 | if i==i_max-1: |
---|
2203 | legend_list=[] |
---|
2204 | |
---|
2205 | #print 'DIC',legend_list_dic |
---|
2206 | for name_and_elevation in legend_list_dic[filename]: |
---|
2207 | legend_list.append('%s (elevation = %sm)'%(name_and_elevation[0],name_and_elevation[1])) |
---|
2208 | |
---|
2209 | #print 'filename',filename, quantities |
---|
2210 | #skip time and elevation so it is not plotted! |
---|
2211 | for k, quantity in enumerate(quantities): |
---|
2212 | if quantity != 'time' and quantity != 'elevation': |
---|
2213 | |
---|
2214 | num=int(k*100+j) |
---|
2215 | pylab.figure(num) |
---|
2216 | pylab.ylabel(quantities_label[quantity]) |
---|
2217 | pylab.plot(t, directory_quantity_value[directory][filename][quantity], c = cstr[i], linewidth=1) |
---|
2218 | pylab.xlabel(quantities_label['time']) |
---|
2219 | pylab.axis(quantities_axis[quantity]) |
---|
2220 | pylab.legend(legend_list,loc='upper right') |
---|
2221 | |
---|
2222 | pylab.title('%s at %s gauge' %(quantity,filename[len(base_name):])) |
---|
2223 | if output_dir == '': |
---|
2224 | figname = '%s%s%s_%s%s.png' %(directory,sep, |
---|
2225 | filename, |
---|
2226 | quantity, |
---|
2227 | extra_plot_name) |
---|
2228 | else: |
---|
2229 | figname = '%s%s%s_%s%s.png' %(output_dir,sep, |
---|
2230 | filename, |
---|
2231 | quantity, |
---|
2232 | extra_plot_name) |
---|
2233 | if verbose: print 'saving figure here %s' %figname |
---|
2234 | pylab.savefig(figname) |
---|
2235 | |
---|
2236 | if verbose: print 'Closing all plots' |
---|
2237 | pylab.close('all') |
---|
2238 | del pylab |
---|
2239 | if verbose: print 'Finished closing plots' |
---|
2240 | |
---|
2241 | def get_min_max_values(list=None): |
---|
2242 | """ |
---|
2243 | Returns the min and max of the list it was provided. |
---|
2244 | """ |
---|
2245 | if list == None: print 'List must be provided' |
---|
2246 | |
---|
2247 | return min(list), max(list) |
---|
2248 | |
---|
2249 | |
---|
2250 | |
---|
2251 | def get_runup_data_for_locations_from_file(gauge_filename, |
---|
2252 | sww_filename, |
---|
2253 | runup_filename, |
---|
2254 | size=10, |
---|
2255 | verbose=False): |
---|
2256 | |
---|
2257 | """this will read a csv file with the header x,y. Then look in a square |
---|
2258 | 'size'x2 around this position for the 'max_inundaiton_height' in the |
---|
2259 | 'sww_filename' and report the findings in the 'runup_filename |
---|
2260 | |
---|
2261 | WARNING: NO TESTS! |
---|
2262 | """ |
---|
2263 | |
---|
2264 | from anuga.shallow_water.data_manager import get_all_directories_with_name,\ |
---|
2265 | get_maximum_inundation_data,\ |
---|
2266 | csv2dict |
---|
2267 | |
---|
2268 | file = open(runup_filename,"w") |
---|
2269 | file.write("easting,northing,runup \n ") |
---|
2270 | file.close() |
---|
2271 | |
---|
2272 | #read gauge csv file to dictionary |
---|
2273 | attribute_dic, title_index_dic = csv2dict(gauge_filename) |
---|
2274 | northing = [float(x) for x in attribute_dic["y"]] |
---|
2275 | easting = [float(x) for x in attribute_dic["x"]] |
---|
2276 | |
---|
2277 | print 'Reading %s' %sww_filename |
---|
2278 | |
---|
2279 | runup_locations=[] |
---|
2280 | for i, x in enumerate(northing): |
---|
2281 | # print 'easting,northing',i,easting[i],northing[i] |
---|
2282 | poly = [[int(easting[i]+size),int(northing[i]+size)], |
---|
2283 | [int(easting[i]+size),int(northing[i]-size)], |
---|
2284 | [int(easting[i]-size),int(northing[i]-size)], |
---|
2285 | [int(easting[i]-size),int(northing[i]+size)]] |
---|
2286 | |
---|
2287 | run_up, x_y = get_maximum_inundation_data(filename=sww_filename, |
---|
2288 | polygon=poly, |
---|
2289 | verbose=False) |
---|
2290 | #if no runup will return 0 instead of NONE |
---|
2291 | if run_up==None: run_up=0 |
---|
2292 | if x_y==None: x_y=[0,0] |
---|
2293 | |
---|
2294 | if verbose: |
---|
2295 | print 'maximum inundation runup near %s is %s meters' %(x_y,run_up) |
---|
2296 | |
---|
2297 | #writes to file |
---|
2298 | file = open(runup_filename,"a") |
---|
2299 | temp = '%s,%s,%s \n' %(x_y[0], x_y[1], run_up) |
---|
2300 | file.write(temp) |
---|
2301 | file.close() |
---|
2302 | |
---|
2303 | |
---|
2304 | def sww2csv_gauges(sww_file, |
---|
2305 | gauge_file, |
---|
2306 | out_name='gauge_', |
---|
2307 | quantities = ['stage', 'depth', 'elevation', |
---|
2308 | 'xmomentum', 'ymomentum'], |
---|
2309 | verbose=False, |
---|
2310 | use_cache = True): |
---|
2311 | """ |
---|
2312 | |
---|
2313 | Inputs: |
---|
2314 | |
---|
2315 | NOTE: if using csv2timeseries_graphs after creating csv file, |
---|
2316 | it is essential to export quantities 'depth' and 'elevation'. |
---|
2317 | 'depth' is good to analyse gauges on land and elevation is used |
---|
2318 | automatically by csv2timeseries_graphs in the legend. |
---|
2319 | |
---|
2320 | sww_file: path to any sww file |
---|
2321 | |
---|
2322 | gauge_file: Assumes that it follows this format |
---|
2323 | name, easting, northing, elevation |
---|
2324 | point1, 100.3, 50.2, 10.0 |
---|
2325 | point2, 10.3, 70.3, 78.0 |
---|
2326 | |
---|
2327 | NOTE: order of column can change but names eg 'easting', 'elevation' |
---|
2328 | must be the same! ALL lowercaps! |
---|
2329 | |
---|
2330 | out_name: prefix for output file name (default is 'gauge_') |
---|
2331 | |
---|
2332 | Outputs: |
---|
2333 | one file for each gauge/point location in the points file. They |
---|
2334 | will be named with this format in the same directory as the 'sww_file' |
---|
2335 | <out_name><name>.csv |
---|
2336 | eg gauge_point1.csv if <out_name> not supplied |
---|
2337 | myfile_2_point1.csv if <out_name> ='myfile_2_' |
---|
2338 | |
---|
2339 | |
---|
2340 | They will all have a header |
---|
2341 | |
---|
2342 | Usage: sww2csv_gauges(sww_file='test1.sww', |
---|
2343 | quantities = ['stage', 'elevation','depth','bearing'], |
---|
2344 | gauge_file='gauge.txt') |
---|
2345 | |
---|
2346 | Interpolate the quantities at a given set of locations, given |
---|
2347 | an sww file. |
---|
2348 | The results are written to a csv file. |
---|
2349 | |
---|
2350 | In the future let points be a points file. |
---|
2351 | And the user choose the quantities. |
---|
2352 | |
---|
2353 | This is currently quite specific. |
---|
2354 | If it needs to be more general, change things. |
---|
2355 | |
---|
2356 | This is really returning speed, not velocity. |
---|
2357 | |
---|
2358 | |
---|
2359 | """ |
---|
2360 | |
---|
2361 | from csv import reader,writer |
---|
2362 | from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN |
---|
2363 | from Numeric import array, resize, shape, Float, zeros, take, argsort, argmin |
---|
2364 | import string |
---|
2365 | from anuga.shallow_water.data_manager import get_all_swwfiles |
---|
2366 | # quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
2367 | #print "points",points |
---|
2368 | |
---|
2369 | assert type(gauge_file) == type(''),\ |
---|
2370 | 'Gauge filename must be a string' |
---|
2371 | |
---|
2372 | assert type(out_name) == type(''),\ |
---|
2373 | 'Output filename prefix must be a string' |
---|
2374 | |
---|
2375 | try: |
---|
2376 | # fid = open(gauge_file) |
---|
2377 | point_reader = reader(file(gauge_file)) |
---|
2378 | |
---|
2379 | except Exception, e: |
---|
2380 | msg = 'File "%s" could not be opened: Error="%s"'\ |
---|
2381 | %(gauge_file, e) |
---|
2382 | raise msg |
---|
2383 | if verbose: print '\n Gauges obtained from: %s \n' %gauge_file |
---|
2384 | |
---|
2385 | |
---|
2386 | point_reader = reader(file(gauge_file)) |
---|
2387 | points = [] |
---|
2388 | point_name = [] |
---|
2389 | |
---|
2390 | #read point info from file |
---|
2391 | for i,row in enumerate(point_reader): |
---|
2392 | # print 'i',i,'row',row |
---|
2393 | #read header and determine the column numbers to read correcty. |
---|
2394 | if i==0: |
---|
2395 | for j,value in enumerate(row): |
---|
2396 | # print 'j',j,value, row |
---|
2397 | if value.strip()=='easting':easting=j |
---|
2398 | if value.strip()=='northing':northing=j |
---|
2399 | if value.strip()=='name':name=j |
---|
2400 | if value.strip()=='elevation':elevation=j |
---|
2401 | else: |
---|
2402 | # print i,'easting',easting,'northing',northing, row[easting] |
---|
2403 | points.append([float(row[easting]),float(row[northing])]) |
---|
2404 | point_name.append(row[name]) |
---|
2405 | |
---|
2406 | #convert to array for file_function |
---|
2407 | points_array = array(points,Float) |
---|
2408 | |
---|
2409 | points_array = ensure_absolute(points_array) |
---|
2410 | |
---|
2411 | dir_name, base = os.path.split(sww_file) |
---|
2412 | #print 'dirname',dir_name, base |
---|
2413 | #need to get current directory so when path and file |
---|
2414 | #are "joined" below the directory is correct |
---|
2415 | if dir_name == '': |
---|
2416 | dir_name =getcwd() |
---|
2417 | |
---|
2418 | if access(sww_file,R_OK): |
---|
2419 | if verbose: print 'File %s exists' %(sww_file) |
---|
2420 | else: |
---|
2421 | msg = 'File "%s" could not be opened: no read permission'\ |
---|
2422 | %(sww_file) |
---|
2423 | raise msg |
---|
2424 | |
---|
2425 | sww_files = get_all_swwfiles(look_in_dir=dir_name, |
---|
2426 | base_name=base, |
---|
2427 | verbose=verbose) |
---|
2428 | |
---|
2429 | #to make all the quantities lower case for file_function |
---|
2430 | quantities = [quantity.lower() for quantity in quantities] |
---|
2431 | |
---|
2432 | # what is quantities are needed from sww file to calculate output quantities |
---|
2433 | # also |
---|
2434 | |
---|
2435 | core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
2436 | |
---|
2437 | for sww_file in sww_files: |
---|
2438 | |
---|
2439 | # print 'sww_file',sww_file |
---|
2440 | sww_file = join(dir_name, sww_file+'.sww') |
---|
2441 | # print 'sww_file',sww_file, core_quantities |
---|
2442 | |
---|
2443 | callable_sww = file_function(sww_file, |
---|
2444 | quantities=core_quantities, |
---|
2445 | interpolation_points=points_array, |
---|
2446 | verbose=verbose, |
---|
2447 | use_cache=use_cache) |
---|
2448 | |
---|
2449 | gauge_file = out_name |
---|
2450 | |
---|
2451 | heading = [quantity for quantity in quantities] |
---|
2452 | heading.insert(0,'time') |
---|
2453 | |
---|
2454 | # print 'start time', callable_sww.starttime, heading, quantities |
---|
2455 | |
---|
2456 | #create a list of csv writers for all the points and write header |
---|
2457 | points_writer = [] |
---|
2458 | for i,point in enumerate(points): |
---|
2459 | #print 'gauge file:',dir_name+sep+'gauge_'+point_name[i]+'.csv' |
---|
2460 | points_writer.append(writer(file(dir_name+sep+gauge_file+point_name[i]+'.csv', "wb"))) |
---|
2461 | points_writer[i].writerow(heading) |
---|
2462 | |
---|
2463 | if verbose: print 'Writing csv files' |
---|
2464 | |
---|
2465 | for time in callable_sww.get_time(): |
---|
2466 | |
---|
2467 | for point_i, point in enumerate(points_array): |
---|
2468 | #add domain starttime to relative time. |
---|
2469 | points_list = [time+callable_sww.starttime] |
---|
2470 | # print'time',time,'point_i',point_i,point, points_array |
---|
2471 | point_quantities = callable_sww(time,point_i) |
---|
2472 | # print "quantities", point_quantities |
---|
2473 | |
---|
2474 | for quantity in quantities: |
---|
2475 | if quantity==NAN: |
---|
2476 | print 'quantity does not exist in' %callable_sww.get_name |
---|
2477 | else: |
---|
2478 | if quantity == 'stage': |
---|
2479 | points_list.append(point_quantities[0]) |
---|
2480 | |
---|
2481 | if quantity == 'elevation': |
---|
2482 | points_list.append(point_quantities[1]) |
---|
2483 | |
---|
2484 | if quantity == 'xmomentum': |
---|
2485 | points_list.append(point_quantities[2]) |
---|
2486 | |
---|
2487 | if quantity == 'ymomentum': |
---|
2488 | points_list.append(point_quantities[3]) |
---|
2489 | |
---|
2490 | if quantity == 'depth': |
---|
2491 | points_list.append(point_quantities[0] - point_quantities[1]) |
---|
2492 | |
---|
2493 | |
---|
2494 | if quantity == 'momentum': |
---|
2495 | momentum = sqrt(point_quantities[2]**2 +\ |
---|
2496 | point_quantities[3]**2) |
---|
2497 | points_list.append(momentum) |
---|
2498 | |
---|
2499 | if quantity == 'speed': |
---|
2500 | #if depth is less than 0.001 then speed = 0.0 |
---|
2501 | if point_quantities[0] - point_quantities[1] < 0.001: |
---|
2502 | vel = 0.0 |
---|
2503 | else: |
---|
2504 | if point_quantities[2] < 1.0e6: |
---|
2505 | momentum = sqrt(point_quantities[2]**2 +\ |
---|
2506 | point_quantities[3]**2) |
---|
2507 | # vel = momentum/depth |
---|
2508 | vel = momentum/(point_quantities[0] - point_quantities[1]) |
---|
2509 | # vel = momentum/(depth + 1.e-6/depth) |
---|
2510 | else: |
---|
2511 | momentum = 0 |
---|
2512 | vel = 0 |
---|
2513 | |
---|
2514 | points_list.append(vel) |
---|
2515 | |
---|
2516 | if quantity == 'bearing': |
---|
2517 | points_list.append(calc_bearing(point_quantities[2], |
---|
2518 | point_quantities[3])) |
---|
2519 | |
---|
2520 | #print 'list',points_list |
---|
2521 | |
---|
2522 | points_writer[point_i].writerow(points_list) |
---|
2523 | |
---|
2524 | |
---|
2525 | def greens_law(d1,d2,h1,verbose=False): |
---|
2526 | """ |
---|
2527 | |
---|
2528 | Green's Law allows an approximation of wave amplitude at |
---|
2529 | a given depth based on the fourh root of the ratio of two depths |
---|
2530 | and the amplitude at another given depth. |
---|
2531 | |
---|
2532 | Note, wave amplitude is equal to stage. |
---|
2533 | |
---|
2534 | Inputs: |
---|
2535 | |
---|
2536 | d1, d2 - the two depths |
---|
2537 | h1 - the wave amplitude at d1 |
---|
2538 | h2 - the derived amplitude at d2 |
---|
2539 | |
---|
2540 | h2 = h1 (d1/d2)^(1/4), where d2 cannot equal 0. |
---|
2541 | |
---|
2542 | """ |
---|
2543 | |
---|
2544 | d1 = ensure_numeric(d1) |
---|
2545 | d2 = ensure_numeric(d2) |
---|
2546 | h1 = ensure_numeric(h1) |
---|
2547 | |
---|
2548 | if d1 <= 0.0: |
---|
2549 | msg = 'the first depth, d1 (%f), must be strictly positive' %(d1) |
---|
2550 | raise Exception(msg) |
---|
2551 | |
---|
2552 | if d2 <= 0.0: |
---|
2553 | msg = 'the second depth, d2 (%f), must be strictly positive' %(d2) |
---|
2554 | raise Exception(msg) |
---|
2555 | |
---|
2556 | if h1 <= 0.0: |
---|
2557 | msg = 'the wave amplitude, h1 (%f), must be strictly positive' %(h1) |
---|
2558 | raise Exception(msg) |
---|
2559 | |
---|
2560 | h2 = h1*(d1/d2)**0.25 |
---|
2561 | |
---|
2562 | assert h2 > 0 |
---|
2563 | |
---|
2564 | return h2 |
---|
2565 | |
---|
2566 | |
---|
2567 | def square_root(s): |
---|
2568 | return sqrt(s) |
---|