1 | """This module contains various auxiliary function used by pyvolution. |
---|
2 | |
---|
3 | It is also a clearing house for functions that may later earn a module |
---|
4 | of their own. |
---|
5 | """ |
---|
6 | |
---|
7 | import anuga.utilities.polygon |
---|
8 | import sys |
---|
9 | import os |
---|
10 | |
---|
11 | from os import remove, mkdir, access, F_OK, W_OK, sep,mkdir |
---|
12 | from os.path import exists, basename, split,join |
---|
13 | from warnings import warn |
---|
14 | from shutil import copy |
---|
15 | |
---|
16 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
17 | from Numeric import arange, choose, zeros, Float, array |
---|
18 | |
---|
19 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
20 | |
---|
21 | def file_function(filename, |
---|
22 | domain=None, |
---|
23 | quantities=None, |
---|
24 | interpolation_points=None, |
---|
25 | time_thinning=1, |
---|
26 | verbose=False, |
---|
27 | use_cache=False): |
---|
28 | """Read time history of spatial data from NetCDF file and return |
---|
29 | a callable object. |
---|
30 | |
---|
31 | Input variables: |
---|
32 | |
---|
33 | filename - Name of sww or tms file |
---|
34 | |
---|
35 | If the file has extension 'sww' then it is assumed to be spatio-temporal |
---|
36 | or temporal and the callable object will have the form f(t,x,y) or f(t) |
---|
37 | depending on whether the file contains spatial data |
---|
38 | |
---|
39 | If the file has extension 'tms' then it is assumed to be temporal only |
---|
40 | and the callable object will have the form f(t) |
---|
41 | |
---|
42 | Either form will return interpolated values based on the input file |
---|
43 | using the underlying interpolation_function. |
---|
44 | |
---|
45 | domain - Associated domain object |
---|
46 | If domain is specified, model time (domain.starttime) |
---|
47 | will be checked and possibly modified. |
---|
48 | |
---|
49 | All times are assumed to be in UTC |
---|
50 | |
---|
51 | All spatial information is assumed to be in absolute UTM coordinates. |
---|
52 | |
---|
53 | quantities - the name of the quantity to be interpolated or a |
---|
54 | list of quantity names. The resulting function will return |
---|
55 | a tuple of values - one for each quantity |
---|
56 | If quantities are None, domain's conserved quantities |
---|
57 | are used. |
---|
58 | |
---|
59 | interpolation_points - list of absolute UTM coordinates for points (N x 2) |
---|
60 | or geospatial object or points file name at which values are sought |
---|
61 | |
---|
62 | use_cache: True means that caching of intermediate result of |
---|
63 | Interpolation_function is attempted |
---|
64 | |
---|
65 | |
---|
66 | See Interpolation function for further documentation |
---|
67 | """ |
---|
68 | |
---|
69 | |
---|
70 | #FIXME (OLE): Should check origin of domain against that of file |
---|
71 | #In fact, this is where origin should be converted to that of domain |
---|
72 | #Also, check that file covers domain fully. |
---|
73 | |
---|
74 | #Take into account: |
---|
75 | #- domain's georef |
---|
76 | #- sww file's georef |
---|
77 | #- interpolation points as absolute UTM coordinates |
---|
78 | |
---|
79 | |
---|
80 | |
---|
81 | |
---|
82 | # Use domain's conserved_quantity names as defaults |
---|
83 | if domain is not None: |
---|
84 | if quantities is None: |
---|
85 | quantities = domain.conserved_quantities |
---|
86 | |
---|
87 | domain_starttime = domain.starttime |
---|
88 | else: |
---|
89 | domain_starttime = None |
---|
90 | |
---|
91 | |
---|
92 | # Build arguments and keyword arguments for use with caching or apply. |
---|
93 | args = (filename,) |
---|
94 | |
---|
95 | |
---|
96 | # FIXME (Ole): Caching this function will not work well |
---|
97 | # if domain is passed in as instances change hash code. |
---|
98 | # Instead we pass in those attributes that are needed (and return them |
---|
99 | # if modified) |
---|
100 | kwargs = {'quantities': quantities, |
---|
101 | 'interpolation_points': interpolation_points, |
---|
102 | 'domain_starttime': domain_starttime, |
---|
103 | 'time_thinning': time_thinning, |
---|
104 | 'verbose': verbose} |
---|
105 | |
---|
106 | |
---|
107 | # Call underlying engine with or without caching |
---|
108 | if use_cache is True: |
---|
109 | try: |
---|
110 | from caching import cache |
---|
111 | except: |
---|
112 | msg = 'Caching was requested, but caching module'+\ |
---|
113 | 'could not be imported' |
---|
114 | raise msg |
---|
115 | |
---|
116 | f, starttime = cache(_file_function, |
---|
117 | args, kwargs, |
---|
118 | dependencies=[filename], |
---|
119 | compression=False, |
---|
120 | verbose=verbose) |
---|
121 | |
---|
122 | else: |
---|
123 | f, starttime = apply(_file_function, |
---|
124 | args, kwargs) |
---|
125 | |
---|
126 | |
---|
127 | #FIXME (Ole): Pass cache arguments, such as compression, in some sort of |
---|
128 | #structure |
---|
129 | |
---|
130 | f.starttime = starttime |
---|
131 | |
---|
132 | if domain is not None: |
---|
133 | #Update domain.startime if it is *earlier* than starttime from file |
---|
134 | if starttime > domain.starttime: |
---|
135 | msg = 'WARNING: Start time as specified in domain (%f)'\ |
---|
136 | %domain.starttime |
---|
137 | msg += ' is earlier than the starttime of file %s (%f).'\ |
---|
138 | %(filename, starttime) |
---|
139 | msg += ' Modifying domain starttime accordingly.' |
---|
140 | |
---|
141 | if verbose: print msg |
---|
142 | domain.starttime = starttime #Modifying model time |
---|
143 | if verbose: print 'Domain starttime is now set to %f'\ |
---|
144 | %domain.starttime |
---|
145 | |
---|
146 | return f |
---|
147 | |
---|
148 | |
---|
149 | |
---|
150 | def _file_function(filename, |
---|
151 | quantities=None, |
---|
152 | interpolation_points=None, |
---|
153 | domain_starttime=None, |
---|
154 | time_thinning=1, |
---|
155 | verbose=False): |
---|
156 | """Internal function |
---|
157 | |
---|
158 | See file_function for documentatiton |
---|
159 | """ |
---|
160 | |
---|
161 | |
---|
162 | assert type(filename) == type(''),\ |
---|
163 | 'First argument to File_function must be a string' |
---|
164 | |
---|
165 | try: |
---|
166 | fid = open(filename) |
---|
167 | except Exception, e: |
---|
168 | msg = 'File "%s" could not be opened: Error="%s"'\ |
---|
169 | %(filename, e) |
---|
170 | raise msg |
---|
171 | |
---|
172 | line = fid.readline() |
---|
173 | fid.close() |
---|
174 | |
---|
175 | |
---|
176 | if line[:3] == 'CDF': |
---|
177 | return get_netcdf_file_function(filename, |
---|
178 | quantities, |
---|
179 | interpolation_points, |
---|
180 | domain_starttime, |
---|
181 | time_thinning=time_thinning, |
---|
182 | verbose=verbose) |
---|
183 | else: |
---|
184 | raise 'Must be a NetCDF File' |
---|
185 | |
---|
186 | |
---|
187 | |
---|
188 | def get_netcdf_file_function(filename, |
---|
189 | quantity_names=None, |
---|
190 | interpolation_points=None, |
---|
191 | domain_starttime=None, |
---|
192 | time_thinning=1, |
---|
193 | verbose=False): |
---|
194 | """Read time history of spatial data from NetCDF sww file and |
---|
195 | return a callable object f(t,x,y) |
---|
196 | which will return interpolated values based on the input file. |
---|
197 | |
---|
198 | Model time (domain_starttime) |
---|
199 | will be checked, possibly modified and returned |
---|
200 | |
---|
201 | All times are assumed to be in UTC |
---|
202 | |
---|
203 | See Interpolation function for further documetation |
---|
204 | |
---|
205 | """ |
---|
206 | |
---|
207 | |
---|
208 | #FIXME: Check that model origin is the same as file's origin |
---|
209 | #(both in UTM coordinates) |
---|
210 | #If not - modify those from file to match domain |
---|
211 | #(origin should be passed in) |
---|
212 | #Take this code from e.g. dem2pts in data_manager.py |
---|
213 | #FIXME: Use geo_reference to read and write xllcorner... |
---|
214 | |
---|
215 | |
---|
216 | import time, calendar, types |
---|
217 | from anuga.config import time_format |
---|
218 | from Scientific.IO.NetCDF import NetCDFFile |
---|
219 | from Numeric import array, zeros, Float, alltrue, concatenate, reshape |
---|
220 | |
---|
221 | #Open NetCDF file |
---|
222 | if verbose: print 'Reading', filename |
---|
223 | fid = NetCDFFile(filename, 'r') |
---|
224 | |
---|
225 | if type(quantity_names) == types.StringType: |
---|
226 | quantity_names = [quantity_names] |
---|
227 | |
---|
228 | if quantity_names is None or len(quantity_names) < 1: |
---|
229 | #If no quantities are specified get quantities from file |
---|
230 | #x, y, time are assumed as the independent variables so |
---|
231 | #they are excluded as potentiol quantities |
---|
232 | quantity_names = [] |
---|
233 | for name in fid.variables.keys(): |
---|
234 | if name not in ['x', 'y', 'time']: |
---|
235 | quantity_names.append(name) |
---|
236 | |
---|
237 | if len(quantity_names) < 1: |
---|
238 | msg = 'ERROR: At least one independent value must be specified' |
---|
239 | raise msg |
---|
240 | |
---|
241 | |
---|
242 | if interpolation_points is not None: |
---|
243 | interpolation_points = ensure_absolute(interpolation_points) |
---|
244 | msg = 'Points must by N x 2. I got %d' %interpolation_points.shape[1] |
---|
245 | assert interpolation_points.shape[1] == 2, msg |
---|
246 | |
---|
247 | |
---|
248 | #Now assert that requested quantitites (and the independent ones) |
---|
249 | #are present in file |
---|
250 | missing = [] |
---|
251 | for quantity in ['time'] + quantity_names: |
---|
252 | if not fid.variables.has_key(quantity): |
---|
253 | missing.append(quantity) |
---|
254 | |
---|
255 | if len(missing) > 0: |
---|
256 | msg = 'Quantities %s could not be found in file %s'\ |
---|
257 | %(str(missing), filename) |
---|
258 | fid.close() |
---|
259 | raise Exception, msg |
---|
260 | |
---|
261 | #Decide whether this data has a spatial dimension |
---|
262 | spatial = True |
---|
263 | for quantity in ['x', 'y']: |
---|
264 | if not fid.variables.has_key(quantity): |
---|
265 | spatial = False |
---|
266 | |
---|
267 | if filename[-3:] == 'tms' and spatial is True: |
---|
268 | msg = 'Files of type tms must not contain spatial information' |
---|
269 | raise msg |
---|
270 | |
---|
271 | if filename[-3:] == 'sww' and spatial is False: |
---|
272 | msg = 'Files of type sww must contain spatial information' |
---|
273 | raise msg |
---|
274 | |
---|
275 | #Get first timestep |
---|
276 | try: |
---|
277 | starttime = fid.starttime[0] |
---|
278 | except ValueError: |
---|
279 | msg = 'Could not read starttime from file %s' %filename |
---|
280 | raise msg |
---|
281 | |
---|
282 | #Get variables |
---|
283 | #if verbose: print 'Get variables' |
---|
284 | time = fid.variables['time'][:] |
---|
285 | |
---|
286 | # Get time independent stuff |
---|
287 | if spatial: |
---|
288 | #Get origin |
---|
289 | xllcorner = fid.xllcorner[0] |
---|
290 | yllcorner = fid.yllcorner[0] |
---|
291 | zone = fid.zone[0] |
---|
292 | |
---|
293 | x = fid.variables['x'][:] |
---|
294 | y = fid.variables['y'][:] |
---|
295 | triangles = fid.variables['volumes'][:] |
---|
296 | |
---|
297 | x = reshape(x, (len(x),1)) |
---|
298 | y = reshape(y, (len(y),1)) |
---|
299 | vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array |
---|
300 | |
---|
301 | if interpolation_points is not None: |
---|
302 | #Adjust for georef |
---|
303 | interpolation_points[:,0] -= xllcorner |
---|
304 | interpolation_points[:,1] -= yllcorner |
---|
305 | |
---|
306 | |
---|
307 | |
---|
308 | |
---|
309 | if domain_starttime is not None: |
---|
310 | |
---|
311 | #If domain_startime is *later* than starttime, |
---|
312 | #move time back - relative to domain's time |
---|
313 | if domain_starttime > starttime: |
---|
314 | time = time - domain_starttime + starttime |
---|
315 | |
---|
316 | #FIXME Use method in geo to reconcile |
---|
317 | #if spatial: |
---|
318 | #assert domain.geo_reference.xllcorner == xllcorner |
---|
319 | #assert domain.geo_reference.yllcorner == yllcorner |
---|
320 | #assert domain.geo_reference.zone == zone |
---|
321 | |
---|
322 | if verbose: |
---|
323 | print 'File_function data obtained from: %s' %filename |
---|
324 | print ' References:' |
---|
325 | #print ' Datum ....' #FIXME |
---|
326 | if spatial: |
---|
327 | print ' Lower left corner: [%f, %f]'\ |
---|
328 | %(xllcorner, yllcorner) |
---|
329 | print ' Start time: %f' %starttime |
---|
330 | |
---|
331 | |
---|
332 | # Produce values for desired data points at |
---|
333 | # each timestep for each quantity |
---|
334 | quantities = {} |
---|
335 | for i, name in enumerate(quantity_names): |
---|
336 | quantities[name] = fid.variables[name][:] |
---|
337 | fid.close() |
---|
338 | |
---|
339 | |
---|
340 | #from least_squares import Interpolation_function |
---|
341 | from anuga.fit_interpolate.interpolate import Interpolation_function |
---|
342 | |
---|
343 | if not spatial: |
---|
344 | vertex_coordinates = triangles = interpolation_points = None |
---|
345 | |
---|
346 | |
---|
347 | # Return Interpolation_function instance as well as |
---|
348 | # starttime for use to possible modify that of domain |
---|
349 | return Interpolation_function(time, |
---|
350 | quantities, |
---|
351 | quantity_names, |
---|
352 | vertex_coordinates, |
---|
353 | triangles, |
---|
354 | interpolation_points, |
---|
355 | time_thinning=time_thinning, |
---|
356 | verbose=verbose), starttime |
---|
357 | |
---|
358 | # NOTE (Ole): Caching Interpolation function is too slow as |
---|
359 | # the very long parameters need to be hashed. |
---|
360 | |
---|
361 | |
---|
362 | |
---|
363 | |
---|
364 | |
---|
365 | |
---|
366 | def multiple_replace(text, dictionary): |
---|
367 | """Multiple replace of words in text |
---|
368 | |
---|
369 | text: String to be modified |
---|
370 | dictionary: Mapping of words that are to be substituted |
---|
371 | |
---|
372 | Python Cookbook 3.14 page 88 and page 90 |
---|
373 | """ |
---|
374 | |
---|
375 | import re |
---|
376 | |
---|
377 | #Create a regular expression from all of the dictionary keys |
---|
378 | #matching only entire words |
---|
379 | regex = re.compile(r'\b'+ \ |
---|
380 | r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \ |
---|
381 | r'\b' ) |
---|
382 | |
---|
383 | #For each match, lookup the corresponding value in the dictionary |
---|
384 | return regex.sub(lambda match: dictionary[match.group(0)], text) |
---|
385 | |
---|
386 | |
---|
387 | |
---|
388 | |
---|
389 | def apply_expression_to_dictionary(expression, dictionary):#dictionary): |
---|
390 | """Apply arbitrary expression to values of dictionary |
---|
391 | |
---|
392 | Given an expression in terms of the keys, replace key by the |
---|
393 | corresponding values and evaluate. |
---|
394 | |
---|
395 | |
---|
396 | expression: Arbitrary, e.g. arithmetric, expression relating keys |
---|
397 | from dictionary. |
---|
398 | |
---|
399 | dictionary: Mapping between symbols in expression and objects that |
---|
400 | will be evaluated by expression. |
---|
401 | Values in dictionary must support operators given in |
---|
402 | expression e.g. by overloading |
---|
403 | |
---|
404 | due to a limitation with Numeric, this can not evaluate 0/0 |
---|
405 | In general, the user can fix by adding 1e-30 to the numerator. |
---|
406 | SciPy core can handle this situation. |
---|
407 | """ |
---|
408 | |
---|
409 | import types |
---|
410 | import re |
---|
411 | |
---|
412 | assert isinstance(expression, basestring) |
---|
413 | assert type(dictionary) == types.DictType |
---|
414 | |
---|
415 | #Convert dictionary values to textual representations suitable for eval |
---|
416 | D = {} |
---|
417 | for key in dictionary: |
---|
418 | D[key] = 'dictionary["%s"]' %key |
---|
419 | |
---|
420 | #Perform substitution of variables |
---|
421 | expression = multiple_replace(expression, D) |
---|
422 | |
---|
423 | #Evaluate and return |
---|
424 | try: |
---|
425 | return eval(expression) |
---|
426 | except NameError, e: |
---|
427 | msg = 'Expression "%s" could not be evaluated: %s' %(expression, e) |
---|
428 | raise NameError, msg |
---|
429 | except ValueError, e: |
---|
430 | msg = 'Expression "%s" could not be evaluated: %s' %(expression, e) |
---|
431 | raise ValueError, msg |
---|
432 | |
---|
433 | |
---|
434 | |
---|
435 | |
---|
436 | #################################### |
---|
437 | ####OBSOLETE STUFF |
---|
438 | |
---|
439 | |
---|
440 | def angle(v1, v2): |
---|
441 | """Temporary Interface to new location""" |
---|
442 | |
---|
443 | import anuga.utilities.numerical_tools as NT |
---|
444 | |
---|
445 | msg = 'angle has moved from util.py. ' |
---|
446 | msg += 'Please use "from anuga.utilities.numerical_tools import angle"' |
---|
447 | warn(msg, DeprecationWarning) |
---|
448 | |
---|
449 | return NT.angle(v1, v2) |
---|
450 | |
---|
451 | def anglediff(v0, v1): |
---|
452 | """Temporary Interface to new location""" |
---|
453 | |
---|
454 | import anuga.utilities.numerical_tools as NT |
---|
455 | |
---|
456 | msg = 'anglediff has moved from util.py. ' |
---|
457 | msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"' |
---|
458 | warn(msg, DeprecationWarning) |
---|
459 | |
---|
460 | return NT.anglediff(v0, v1) |
---|
461 | |
---|
462 | |
---|
463 | def mean(x): |
---|
464 | """Temporary Interface to new location""" |
---|
465 | |
---|
466 | import anuga.utilities.numerical_tools as NT |
---|
467 | |
---|
468 | msg = 'mean has moved from util.py. ' |
---|
469 | msg += 'Please use "from anuga.utilities.numerical_tools import mean"' |
---|
470 | warn(msg, DeprecationWarning) |
---|
471 | |
---|
472 | return NT.mean(x) |
---|
473 | |
---|
474 | def point_on_line(*args, **kwargs): |
---|
475 | """Temporary Interface to new location""" |
---|
476 | |
---|
477 | msg = 'point_on_line has moved from util.py. ' |
---|
478 | msg += 'Please use "from anuga.utilities.polygon import point_on_line"' |
---|
479 | warn(msg, DeprecationWarning) |
---|
480 | |
---|
481 | return utilities.polygon.point_on_line(*args, **kwargs) |
---|
482 | |
---|
483 | def inside_polygon(*args, **kwargs): |
---|
484 | """Temporary Interface to new location""" |
---|
485 | |
---|
486 | print 'inside_polygon has moved from util.py. ', |
---|
487 | print 'Please use "from anuga.utilities.polygon import inside_polygon"' |
---|
488 | |
---|
489 | return utilities.polygon.inside_polygon(*args, **kwargs) |
---|
490 | |
---|
491 | def outside_polygon(*args, **kwargs): |
---|
492 | """Temporary Interface to new location""" |
---|
493 | |
---|
494 | print 'outside_polygon has moved from util.py. ', |
---|
495 | print 'Please use "from anuga.utilities.polygon import outside_polygon"' |
---|
496 | |
---|
497 | return utilities.polygon.outside_polygon(*args, **kwargs) |
---|
498 | |
---|
499 | |
---|
500 | def separate_points_by_polygon(*args, **kwargs): |
---|
501 | """Temporary Interface to new location""" |
---|
502 | |
---|
503 | print 'separate_points_by_polygon has moved from util.py. ', |
---|
504 | print 'Please use "from anuga.utilities.polygon import separate_points_by_polygon"' |
---|
505 | |
---|
506 | return utilities.polygon.separate_points_by_polygon(*args, **kwargs) |
---|
507 | |
---|
508 | |
---|
509 | |
---|
510 | def read_polygon(*args, **kwargs): |
---|
511 | """Temporary Interface to new location""" |
---|
512 | |
---|
513 | print 'read_polygon has moved from util.py. ', |
---|
514 | print 'Please use "from anuga.utilities.polygon import read_polygon"' |
---|
515 | |
---|
516 | return utilities.polygon.read_polygon(*args, **kwargs) |
---|
517 | |
---|
518 | |
---|
519 | def populate_polygon(*args, **kwargs): |
---|
520 | """Temporary Interface to new location""" |
---|
521 | |
---|
522 | print 'populate_polygon has moved from util.py. ', |
---|
523 | print 'Please use "from anuga.utilities.polygon import populate_polygon"' |
---|
524 | |
---|
525 | return utilities.polygon.populate_polygon(*args, **kwargs) |
---|
526 | |
---|
527 | ##################### end of obsolete stuff ? ############ |
---|
528 | |
---|
529 | def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', |
---|
530 | print_to_screen=False, verbose=False): |
---|
531 | """Temporary Interface to new location""" |
---|
532 | from anuga.shallow_water.data_manager import start_screen_catcher as dm_start_screen_catcher |
---|
533 | |
---|
534 | print 'start_screen_catcher has moved from util.py. ', |
---|
535 | print 'Please use "from anuga.shallow_water.data_manager import start_screen_catcher"' |
---|
536 | |
---|
537 | return dm_start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', |
---|
538 | print_to_screen=False, verbose=False) |
---|
539 | |
---|
540 | def get_revision_number(): |
---|
541 | """Get the version number of the SVN |
---|
542 | NOTE: This requires that the command svn is on the system PATH |
---|
543 | (simply aliasing svn to the binary will not work) |
---|
544 | """ |
---|
545 | |
---|
546 | # Create dummy info |
---|
547 | #info = 'Revision: Version info could not be obtained.' |
---|
548 | #info += 'A command line version of svn must be availbable ' |
---|
549 | #info += 'on the system PATH, access to the subversion ' |
---|
550 | #info += 'repository is necessary and the output must ' |
---|
551 | #info += 'contain a line starting with "Revision:"' |
---|
552 | |
---|
553 | |
---|
554 | #FIXME (Ole): Change this so that svn info is attempted first. |
---|
555 | # If that fails, try to read a stored file with that same info (this would be created by e.g. the release script). Failing that, throw an exception. |
---|
556 | |
---|
557 | #FIXME (Ole): Move this and store_version_info to utilities |
---|
558 | |
---|
559 | |
---|
560 | try: |
---|
561 | from anuga.stored_version_info import version_info |
---|
562 | except: |
---|
563 | msg = 'No version info stored and command "svn" is not ' |
---|
564 | msg += 'recognised on the system PATH.\n\n' |
---|
565 | msg += 'If ANUGA has been installed from a distribution e.g. as ' |
---|
566 | msg += 'obtained from SourceForge,\n' |
---|
567 | msg += 'the version info should be ' |
---|
568 | msg += 'available in the automatically generated file ' |
---|
569 | msg += 'stored_version_info.py\n' |
---|
570 | msg += 'in the anuga root directory.\n' |
---|
571 | msg += 'If run from a Subversion sandpit, ' |
---|
572 | msg += 'ANUGA will try to obtain the version info ' |
---|
573 | msg += 'by using the command: "svn info".\n' |
---|
574 | msg += 'In this case, make sure svn is accessible on the system path. ' |
---|
575 | msg += 'Simply aliasing svn to the binary will not work. ' |
---|
576 | msg += 'Good luck!' |
---|
577 | |
---|
578 | # No file available - try using Subversion |
---|
579 | try: |
---|
580 | # The null stuff is so this section fails quitly. |
---|
581 | # This could cause the svn info command to fail due to |
---|
582 | # the redirection being bad on some platforms. |
---|
583 | # If that occurs then change this code. |
---|
584 | if sys.platform[0:3] == 'win': |
---|
585 | fid = os.popen('svn info 2> null') |
---|
586 | else: |
---|
587 | fid = os.popen('svn info 2>/dev/null') |
---|
588 | |
---|
589 | except: |
---|
590 | raise Exception(msg) |
---|
591 | else: |
---|
592 | #print 'Got version from svn' |
---|
593 | version_info = fid.read() |
---|
594 | |
---|
595 | if version_info == '': |
---|
596 | raise Exception(msg) |
---|
597 | else: |
---|
598 | pass |
---|
599 | #print 'Got version from file' |
---|
600 | |
---|
601 | |
---|
602 | for line in version_info.split('\n'): |
---|
603 | if line.startswith('Revision:'): |
---|
604 | break |
---|
605 | |
---|
606 | fields = line.split(':') |
---|
607 | msg = 'Keyword "Revision" was not found anywhere in text: %s' %version_info |
---|
608 | assert fields[0].startswith('Revision'), msg |
---|
609 | |
---|
610 | try: |
---|
611 | revision_number = int(fields[1]) |
---|
612 | except: |
---|
613 | msg = 'Revision number must be an integer. I got %s' %fields[1] |
---|
614 | msg += 'Check that the command svn is on the system path' |
---|
615 | raise Exception(msg) |
---|
616 | |
---|
617 | return revision_number |
---|
618 | |
---|
619 | |
---|
620 | def store_version_info(destination_path='.', verbose=False): |
---|
621 | """Obtain current version from Subversion and store it. |
---|
622 | |
---|
623 | Title: store_version_info() |
---|
624 | |
---|
625 | Author: Ole Nielsen (Ole.Nielsen@ga.gov.au) |
---|
626 | |
---|
627 | CreationDate: January 2006 |
---|
628 | |
---|
629 | Description: |
---|
630 | This function obtains current version from Subversion and stores it |
---|
631 | is a Python file named 'stored_version_info.py' for use with |
---|
632 | get_version_info() |
---|
633 | |
---|
634 | If svn is not available on the system PATH, an Exception is thrown |
---|
635 | """ |
---|
636 | |
---|
637 | # Note (Ole): This function should not be unit tested as it will only |
---|
638 | # work when running out of the sandpit. End users downloading the |
---|
639 | # ANUGA distribution would see a failure. |
---|
640 | # |
---|
641 | # FIXME: This function should really only be used by developers ( |
---|
642 | # (e.g. for creating new ANUGA releases), so maybe it should move |
---|
643 | # to somewhere else. |
---|
644 | |
---|
645 | import config |
---|
646 | |
---|
647 | try: |
---|
648 | fid = os.popen('svn info') |
---|
649 | except: |
---|
650 | msg = 'Command "svn" is not recognised on the system PATH' |
---|
651 | raise Exception(msg) |
---|
652 | else: |
---|
653 | txt = fid.read() |
---|
654 | fid.close() |
---|
655 | |
---|
656 | |
---|
657 | # Determine absolute filename |
---|
658 | if destination_path[-1] != os.sep: |
---|
659 | destination_path += os.sep |
---|
660 | |
---|
661 | filename = destination_path + config.version_filename |
---|
662 | |
---|
663 | fid = open(filename, 'w') |
---|
664 | |
---|
665 | docstring = 'Stored version info.\n\n' |
---|
666 | docstring += 'This file provides the version for distributions ' |
---|
667 | docstring += 'that are not accessing Subversion directly.\n' |
---|
668 | docstring += 'The file is automatically generated and should not ' |
---|
669 | docstring += 'be modified manually.\n' |
---|
670 | fid.write('"""%s"""\n\n' %docstring) |
---|
671 | |
---|
672 | fid.write('version_info = """\n%s"""' %txt) |
---|
673 | fid.close() |
---|
674 | |
---|
675 | |
---|
676 | if verbose is True: |
---|
677 | print 'Version info stored to %s' %filename |
---|
678 | |
---|
679 | |
---|
680 | def sww2timeseries(swwfiles, |
---|
681 | gauge_filename, |
---|
682 | production_dirs, |
---|
683 | report = None, |
---|
684 | reportname = None, |
---|
685 | plot_quantity = None, |
---|
686 | generate_fig = False, |
---|
687 | surface = None, |
---|
688 | time_min = None, |
---|
689 | time_max = None, |
---|
690 | time_thinning = 1, |
---|
691 | time_unit = None, |
---|
692 | title_on = None, |
---|
693 | use_cache = False, |
---|
694 | verbose = False): |
---|
695 | |
---|
696 | """ Read sww file and plot the time series for the |
---|
697 | prescribed quantities at defined gauge locations and |
---|
698 | prescribed time range. |
---|
699 | |
---|
700 | Input variables: |
---|
701 | |
---|
702 | swwfiles - dictionary of sww files with label_ids (used in |
---|
703 | generating latex output. It will be part of |
---|
704 | the directory name of file_loc (typically the timestamp). |
---|
705 | Helps to differentiate latex files for different simulations |
---|
706 | for a particular scenario. |
---|
707 | - assume that all conserved quantities have been stored |
---|
708 | - assume each sww file has been simulated with same timestep |
---|
709 | |
---|
710 | gauge_filename - name of file containing gauge data |
---|
711 | - easting, northing, name , elevation? |
---|
712 | - OR (this is not yet done) |
---|
713 | - structure which can be converted to a Numeric array, |
---|
714 | such as a geospatial data object |
---|
715 | |
---|
716 | production_dirs - A list of list, example {20061101_121212: '1 in 10000', |
---|
717 | 'boundaries': 'urs boundary'} |
---|
718 | this will use the second part as the label and the first part |
---|
719 | as the ? |
---|
720 | #FIXME: Is it a list or a dictionary |
---|
721 | # This is probably obsolete by now |
---|
722 | |
---|
723 | report - if True, then write figures to report_figures directory in |
---|
724 | relevant production directory |
---|
725 | - if False, figures are already stored with sww file |
---|
726 | - default to False |
---|
727 | |
---|
728 | reportname - name for report if wishing to generate report |
---|
729 | |
---|
730 | plot_quantity - list containing quantities to plot, they must |
---|
731 | be the name of an existing quantity or one of |
---|
732 | the following possibilities |
---|
733 | - possibilities: |
---|
734 | - stage; 'stage' |
---|
735 | - depth; 'depth' |
---|
736 | - speed; calculated as absolute momentum |
---|
737 | (pointwise) divided by depth; 'speed' |
---|
738 | - bearing; calculated as the angle of the momentum |
---|
739 | vector (xmomentum, ymomentum) from the North; 'bearing' |
---|
740 | - absolute momentum; calculated as |
---|
741 | sqrt(xmomentum^2 + ymomentum^2); 'momentum' |
---|
742 | - x momentum; 'xmomentum' |
---|
743 | - y momentum; 'ymomentum' |
---|
744 | - default will be ['stage', 'speed', 'bearing'] |
---|
745 | |
---|
746 | generate_fig - if True, generate figures as well as csv files of quantities |
---|
747 | - if False, csv files created only |
---|
748 | |
---|
749 | surface - if True, then generate solution surface with 3d plot |
---|
750 | and save to current working directory |
---|
751 | - default = False |
---|
752 | |
---|
753 | time_min - beginning of user defined time range for plotting purposes |
---|
754 | - default will be first available time found in swwfile |
---|
755 | |
---|
756 | time_max - end of user defined time range for plotting purposes |
---|
757 | - default will be last available time found in swwfile |
---|
758 | |
---|
759 | title_on - if True, export standard graphics with title |
---|
760 | - if False, export standard graphics without title |
---|
761 | |
---|
762 | |
---|
763 | |
---|
764 | Output: |
---|
765 | |
---|
766 | - time series data stored in .csv for later use if required. |
---|
767 | Name = gauges_timeseries followed by gauge name |
---|
768 | - latex file will be generated in same directory as where script is |
---|
769 | run (usually production scenario directory. |
---|
770 | Name = latexoutputlabel_id.tex |
---|
771 | |
---|
772 | Other important information: |
---|
773 | |
---|
774 | It is assumed that the used has stored all the conserved quantities |
---|
775 | and elevation during the scenario run, i.e. |
---|
776 | ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
777 | If this has not occurred then sww2timeseries will not work. |
---|
778 | |
---|
779 | |
---|
780 | Usage example |
---|
781 | texname = sww2timeseries({project.boundary_name + '.sww': ''}, |
---|
782 | project.polygons_dir + sep + 'boundary_extent.csv', |
---|
783 | project.anuga_dir, |
---|
784 | report = False, |
---|
785 | plot_quantity = ['stage', 'speed', 'bearing'], |
---|
786 | time_min = None, |
---|
787 | time_max = None, |
---|
788 | title_on = True, |
---|
789 | verbose = True) |
---|
790 | |
---|
791 | """ |
---|
792 | |
---|
793 | |
---|
794 | k = _sww2timeseries(swwfiles, |
---|
795 | gauge_filename, |
---|
796 | production_dirs, |
---|
797 | report, |
---|
798 | reportname, |
---|
799 | plot_quantity, |
---|
800 | generate_fig, |
---|
801 | surface, |
---|
802 | time_min, |
---|
803 | time_max, |
---|
804 | time_thinning, |
---|
805 | time_unit, |
---|
806 | title_on, |
---|
807 | use_cache, |
---|
808 | verbose) |
---|
809 | |
---|
810 | return k |
---|
811 | |
---|
812 | def _sww2timeseries(swwfiles, |
---|
813 | gauge_filename, |
---|
814 | production_dirs, |
---|
815 | report = None, |
---|
816 | reportname = None, |
---|
817 | plot_quantity = None, |
---|
818 | generate_fig = False, |
---|
819 | surface = None, |
---|
820 | time_min = None, |
---|
821 | time_max = None, |
---|
822 | time_thinning = 1, |
---|
823 | time_unit = None, |
---|
824 | title_on = None, |
---|
825 | use_cache = False, |
---|
826 | verbose = False): |
---|
827 | |
---|
828 | assert type(gauge_filename) == type(''),\ |
---|
829 | 'Gauge filename must be a string' |
---|
830 | |
---|
831 | try: |
---|
832 | fid = open(gauge_filename) |
---|
833 | except Exception, e: |
---|
834 | msg = 'File "%s" could not be opened: Error="%s"'\ |
---|
835 | %(gauge_filename, e) |
---|
836 | raise msg |
---|
837 | |
---|
838 | if report is None: |
---|
839 | report = False |
---|
840 | |
---|
841 | if plot_quantity is None: |
---|
842 | plot_quantity = ['depth', 'speed'] |
---|
843 | else: |
---|
844 | assert type(plot_quantity) == list,\ |
---|
845 | 'plot_quantity must be a list' |
---|
846 | check_list(plot_quantity) |
---|
847 | |
---|
848 | if surface is None: |
---|
849 | surface = False |
---|
850 | |
---|
851 | if time_unit is None: |
---|
852 | time_unit = 'hours' |
---|
853 | |
---|
854 | if title_on is None: |
---|
855 | title_on = True |
---|
856 | |
---|
857 | if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename |
---|
858 | gauges, locations, elev = get_gauges_from_file(gauge_filename) |
---|
859 | |
---|
860 | sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
861 | |
---|
862 | file_loc = [] |
---|
863 | f_list = [] |
---|
864 | label_id = [] |
---|
865 | leg_label = [] |
---|
866 | themaxT = 0.0 |
---|
867 | theminT = 0.0 |
---|
868 | |
---|
869 | for swwfile in swwfiles.keys(): |
---|
870 | |
---|
871 | try: |
---|
872 | fid = open(swwfile) |
---|
873 | except Exception, e: |
---|
874 | msg = 'File "%s" could not be opened: Error="%s"'\ |
---|
875 | %(swwfile, e) |
---|
876 | raise msg |
---|
877 | |
---|
878 | print 'swwfile', swwfile |
---|
879 | |
---|
880 | # Extract parent dir name and use as label |
---|
881 | path, _ = os.path.split(swwfile) |
---|
882 | _, label = os.path.split(path) |
---|
883 | |
---|
884 | #print 'label', label |
---|
885 | leg_label.append(label) |
---|
886 | |
---|
887 | |
---|
888 | |
---|
889 | f = file_function(swwfile, |
---|
890 | quantities = sww_quantity, |
---|
891 | interpolation_points = gauges, |
---|
892 | time_thinning = time_thinning, |
---|
893 | verbose = verbose, |
---|
894 | use_cache = use_cache) |
---|
895 | |
---|
896 | # determine which gauges are contained in sww file |
---|
897 | count = 0 |
---|
898 | gauge_index = [] |
---|
899 | for k, g in enumerate(gauges): |
---|
900 | if f(0.0, point_id = k)[2] > 1.0e6: |
---|
901 | count += 1 |
---|
902 | if count == 1: print 'Gauges not contained here:' |
---|
903 | print locations[k] |
---|
904 | else: |
---|
905 | gauge_index.append(k) |
---|
906 | |
---|
907 | if len(gauge_index) > 0: |
---|
908 | print 'Gauges contained here: \n', |
---|
909 | else: |
---|
910 | print 'No gauges contained here. \n' |
---|
911 | for i in range(len(gauge_index)): |
---|
912 | print locations[gauge_index[i]] |
---|
913 | |
---|
914 | index = swwfile.rfind(sep) |
---|
915 | file_loc.append(swwfile[:index+1]) |
---|
916 | label_id.append(swwfiles[swwfile]) |
---|
917 | |
---|
918 | |
---|
919 | f_list.append(f) |
---|
920 | maxT = max(f.get_time()) |
---|
921 | minT = min(f.get_time()) |
---|
922 | if maxT > themaxT: themaxT = maxT |
---|
923 | if minT > theminT: theminT = minT |
---|
924 | |
---|
925 | if time_min is None: |
---|
926 | time_min = theminT # min(T) |
---|
927 | else: |
---|
928 | if time_min < theminT: # min(T): |
---|
929 | msg = 'Minimum time entered not correct - please try again' |
---|
930 | raise Exception, msg |
---|
931 | |
---|
932 | if time_max is None: |
---|
933 | time_max = themaxT # max(T) |
---|
934 | else: |
---|
935 | if time_max > themaxT: # max(T): |
---|
936 | msg = 'Maximum time entered not correct - please try again' |
---|
937 | raise Exception, msg |
---|
938 | |
---|
939 | if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures' |
---|
940 | |
---|
941 | if len(gauge_index) <> 0: |
---|
942 | texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface, |
---|
943 | leg_label, f_list, gauges, locations, elev, gauge_index, |
---|
944 | production_dirs, time_min, time_max, time_unit, |
---|
945 | title_on, label_id, generate_fig, verbose) |
---|
946 | else: |
---|
947 | texfile = '' |
---|
948 | elev_output = [] |
---|
949 | |
---|
950 | return texfile, elev_output |
---|
951 | |
---|
952 | |
---|
953 | |
---|
954 | #Fixme - Use geospatial to read this file - it's an xya file |
---|
955 | #Need to include other information into this filename, so xya + Name - required for report |
---|
956 | def get_gauges_from_file(filename): |
---|
957 | """ Read in gauge information from file |
---|
958 | """ |
---|
959 | from os import sep, getcwd, access, F_OK, mkdir |
---|
960 | fid = open(filename) |
---|
961 | lines = fid.readlines() |
---|
962 | fid.close() |
---|
963 | |
---|
964 | gauges = [] |
---|
965 | gaugelocation = [] |
---|
966 | elev = [] |
---|
967 | |
---|
968 | # Check header information |
---|
969 | line1 = lines[0] |
---|
970 | line11 = line1.split(',') |
---|
971 | |
---|
972 | if isinstance(line11[0],str) is True: |
---|
973 | # We have found text in the first line |
---|
974 | east_index = None |
---|
975 | north_index = None |
---|
976 | name_index = None |
---|
977 | elev_index = None |
---|
978 | for i in range(len(line11)): |
---|
979 | if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'easting': east_index = i |
---|
980 | if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'northing': north_index = i |
---|
981 | if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'name': name_index = i |
---|
982 | if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'elevation': elev_index = i |
---|
983 | |
---|
984 | if east_index < len(line11) and north_index < len(line11): |
---|
985 | pass |
---|
986 | else: |
---|
987 | msg = 'WARNING: %s does not contain correct header information' %(filename) |
---|
988 | msg += 'The header must be: easting, northing, name, elevation' |
---|
989 | raise Exception, msg |
---|
990 | |
---|
991 | if elev_index is None: |
---|
992 | raise Exception |
---|
993 | |
---|
994 | if name_index is None: |
---|
995 | raise Exception |
---|
996 | |
---|
997 | lines = lines[1:] # Remove header from data |
---|
998 | else: |
---|
999 | # No header, assume that this is a simple easting, northing file |
---|
1000 | |
---|
1001 | msg = 'There was no header in file %s and the number of columns is %d' %(filename, len(line11)) |
---|
1002 | msg += '- I was assuming two columns corresponding to Easting and Northing' |
---|
1003 | assert len(line11) == 2, msg |
---|
1004 | |
---|
1005 | east_index = 0 |
---|
1006 | north_index = 1 |
---|
1007 | |
---|
1008 | N = len(lines) |
---|
1009 | elev = [-9999]*N |
---|
1010 | gaugelocation = range(N) |
---|
1011 | |
---|
1012 | # Read in gauge data |
---|
1013 | for line in lines: |
---|
1014 | fields = line.split(',') |
---|
1015 | |
---|
1016 | gauges.append([float(fields[east_index]), float(fields[north_index])]) |
---|
1017 | |
---|
1018 | if len(fields) > 2: |
---|
1019 | elev.append(float(fields[elev_index])) |
---|
1020 | loc = fields[name_index] |
---|
1021 | gaugelocation.append(loc.strip('\n')) |
---|
1022 | |
---|
1023 | return gauges, gaugelocation, elev |
---|
1024 | |
---|
1025 | |
---|
1026 | |
---|
1027 | def check_list(quantity): |
---|
1028 | """ Check that input quantities in quantity list are possible |
---|
1029 | """ |
---|
1030 | all_quantity = ['stage', 'depth', 'momentum', 'xmomentum', |
---|
1031 | 'ymomentum', 'speed', 'bearing', 'elevation'] |
---|
1032 | |
---|
1033 | |
---|
1034 | |
---|
1035 | import sys |
---|
1036 | #if not sys.version.startswith('2.4'): |
---|
1037 | # Backwards compatibility |
---|
1038 | # from sets import Set as set |
---|
1039 | |
---|
1040 | from sets import Set as set |
---|
1041 | |
---|
1042 | for i,j in enumerate(quantity): |
---|
1043 | quantity[i] = quantity[i].lower() |
---|
1044 | p = list(set(quantity).difference(set(all_quantity))) |
---|
1045 | if len(p) <> 0: |
---|
1046 | msg = 'Quantities %s do not exist - please try again' %p |
---|
1047 | raise Exception, msg |
---|
1048 | |
---|
1049 | return |
---|
1050 | |
---|
1051 | def calc_bearing(uh, vh): |
---|
1052 | """ Calculate velocity bearing from North |
---|
1053 | """ |
---|
1054 | from math import atan, degrees |
---|
1055 | |
---|
1056 | angle = degrees(atan(vh/(uh+1.e-15))) |
---|
1057 | if (0 < angle < 90.0): |
---|
1058 | if vh > 0: |
---|
1059 | bearing = 90.0 - abs(angle) |
---|
1060 | if vh < 0: |
---|
1061 | bearing = 270.0 - abs(angle) |
---|
1062 | if (-90 < angle < 0): |
---|
1063 | if vh < 0: |
---|
1064 | bearing = 90.0 - abs(angle) |
---|
1065 | if vh > 0: |
---|
1066 | bearing = 270.0 - abs(angle) |
---|
1067 | if angle == 0: bearing = 0.0 |
---|
1068 | |
---|
1069 | return bearing |
---|
1070 | |
---|
1071 | def generate_figures(plot_quantity, file_loc, report, reportname, surface, |
---|
1072 | leg_label, f_list, gauges, locations, elev, gauge_index, |
---|
1073 | production_dirs, time_min, time_max, time_unit, |
---|
1074 | title_on, label_id, generate_fig, verbose): |
---|
1075 | """ Generate figures based on required quantities and gauges for each sww file |
---|
1076 | """ |
---|
1077 | from math import sqrt, atan, degrees |
---|
1078 | from Numeric import ones, allclose, zeros, Float, ravel |
---|
1079 | from os import sep, altsep, getcwd, mkdir, access, F_OK, environ |
---|
1080 | |
---|
1081 | if generate_fig is True: |
---|
1082 | from pylab import ion, hold, plot, axis, figure, legend, savefig, \ |
---|
1083 | xlabel, ylabel, title, close, subplot |
---|
1084 | |
---|
1085 | if surface is True: |
---|
1086 | import pylab as p1 |
---|
1087 | import mpl3d.mplot3d as p3 |
---|
1088 | |
---|
1089 | if report == True: |
---|
1090 | texdir = getcwd()+sep+'report'+sep |
---|
1091 | if access(texdir,F_OK) == 0: |
---|
1092 | mkdir (texdir) |
---|
1093 | if len(label_id) == 1: |
---|
1094 | label_id1 = label_id[0].replace(sep,'') |
---|
1095 | label_id2 = label_id1.replace('_','') |
---|
1096 | texfile = texdir+reportname+'%s' %(label_id2) |
---|
1097 | texfile2 = reportname+'%s' %(label_id2) |
---|
1098 | texfilename = texfile + '.tex' |
---|
1099 | if verbose: print '\n Latex output printed to %s \n' %texfilename |
---|
1100 | fid = open(texfilename, 'w') |
---|
1101 | else: |
---|
1102 | texfile = texdir+reportname |
---|
1103 | texfile2 = reportname |
---|
1104 | texfilename = texfile + '.tex' |
---|
1105 | if verbose: print '\n Latex output printed to %s \n' %texfilename |
---|
1106 | fid = open(texfilename, 'w') |
---|
1107 | else: |
---|
1108 | texfile = '' |
---|
1109 | texfile2 = '' |
---|
1110 | |
---|
1111 | p = len(f_list) |
---|
1112 | n = [] |
---|
1113 | n0 = 0 |
---|
1114 | for i in range(len(f_list)): |
---|
1115 | n.append(len(f_list[i].get_time())) |
---|
1116 | if n[i] > n0: n0 = n[i] |
---|
1117 | n0 = int(n0) |
---|
1118 | m = len(locations) |
---|
1119 | model_time = zeros((n0,m,p), Float) |
---|
1120 | stages = zeros((n0,m,p), Float) |
---|
1121 | elevations = zeros((n0,m,p), Float) |
---|
1122 | momenta = zeros((n0,m,p), Float) |
---|
1123 | xmom = zeros((n0,m,p), Float) |
---|
1124 | ymom = zeros((n0,m,p), Float) |
---|
1125 | speed = zeros((n0,m,p), Float) |
---|
1126 | bearings = zeros((n0,m,p), Float) |
---|
1127 | depths = zeros((n0,m,p), Float) |
---|
1128 | eastings = zeros((n0,m,p), Float) |
---|
1129 | min_stages = [] |
---|
1130 | max_stages = [] |
---|
1131 | min_momentums = [] |
---|
1132 | max_momentums = [] |
---|
1133 | max_xmomentums = [] |
---|
1134 | max_ymomentums = [] |
---|
1135 | min_xmomentums = [] |
---|
1136 | min_ymomentums = [] |
---|
1137 | max_speeds = [] |
---|
1138 | min_speeds = [] |
---|
1139 | max_depths = [] |
---|
1140 | model_time_plot3d = zeros((n0,m), Float) |
---|
1141 | stages_plot3d = zeros((n0,m), Float) |
---|
1142 | eastings_plot3d = zeros((n0,m),Float) |
---|
1143 | if time_unit is 'mins': scale = 60.0 |
---|
1144 | if time_unit is 'hours': scale = 3600.0 |
---|
1145 | ##### loop over each swwfile ##### |
---|
1146 | for j, f in enumerate(f_list): |
---|
1147 | starttime = f.starttime |
---|
1148 | if verbose: print 'swwfile %d of %d' %(j, len(f_list)) |
---|
1149 | comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv' |
---|
1150 | fid_compare = open(comparefile, 'w') |
---|
1151 | file0 = file_loc[j]+'gauges_t0.csv' |
---|
1152 | fid_0 = open(file0, 'w') |
---|
1153 | ##### loop over each gauge ##### |
---|
1154 | for k in gauge_index: |
---|
1155 | g = gauges[k] |
---|
1156 | if verbose: print 'Gauge %d of %d' %(k, len(gauges)) |
---|
1157 | min_stage = 10 |
---|
1158 | max_stage = 0 |
---|
1159 | max_momentum = max_xmomentum = max_ymomentum = 0 |
---|
1160 | min_momentum = min_xmomentum = min_ymomentum = 100 |
---|
1161 | max_speed = 0 |
---|
1162 | min_speed = 0 |
---|
1163 | max_depth = 0 |
---|
1164 | gaugeloc = str(locations[k]) |
---|
1165 | thisfile = file_loc[j]+sep+'gauges_time_series'+'_'\ |
---|
1166 | +gaugeloc+'.csv' |
---|
1167 | fid_out = open(thisfile, 'w') |
---|
1168 | s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom \n' |
---|
1169 | fid_out.write(s) |
---|
1170 | #### generate quantities ####### |
---|
1171 | for i, t in enumerate(f.get_time()): |
---|
1172 | if time_min <= t <= time_max: |
---|
1173 | w = f(t, point_id = k)[0] |
---|
1174 | z = f(t, point_id = k)[1] |
---|
1175 | uh = f(t, point_id = k)[2] |
---|
1176 | vh = f(t, point_id = k)[3] |
---|
1177 | depth = w-z |
---|
1178 | m = sqrt(uh*uh + vh*vh) |
---|
1179 | if depth < 0.001: |
---|
1180 | vel = 0.0 |
---|
1181 | else: |
---|
1182 | vel = m / (depth + 1.e-6/depth) |
---|
1183 | #bearing = calc_bearing(uh, vh) |
---|
1184 | model_time[i,k,j] = (t + starttime)/scale #t/60.0 |
---|
1185 | stages[i,k,j] = w |
---|
1186 | elevations[i,k,j] = z |
---|
1187 | xmom[i,k,j] = uh |
---|
1188 | ymom[i,k,j] = vh |
---|
1189 | momenta[i,k,j] = m |
---|
1190 | speed[i,k,j] = vel |
---|
1191 | #bearings[i,k,j] = bearing |
---|
1192 | depths[i,k,j] = depth |
---|
1193 | thisgauge = gauges[k] |
---|
1194 | eastings[i,k,j] = thisgauge[0] |
---|
1195 | s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel, z, uh, vh) |
---|
1196 | fid_out.write(s) |
---|
1197 | if t == 0: |
---|
1198 | s = '%.2f, %.2f, %.2f\n' %(g[0], g[1], w) |
---|
1199 | fid_0.write(s) |
---|
1200 | if t/60.0 <= 13920: tindex = i |
---|
1201 | if w > max_stage: max_stage = w |
---|
1202 | if w < min_stage: min_stage = w |
---|
1203 | if m > max_momentum: max_momentum = m |
---|
1204 | if m < min_momentum: min_momentum = m |
---|
1205 | if uh > max_xmomentum: max_xmomentum = uh |
---|
1206 | if vh > max_ymomentum: max_ymomentum = vh |
---|
1207 | if uh < min_xmomentum: min_xmomentum = uh |
---|
1208 | if vh < min_ymomentum: min_ymomentum = vh |
---|
1209 | if vel > max_speed: max_speed = vel |
---|
1210 | if vel < min_speed: min_speed = vel |
---|
1211 | if z > 0 and depth > max_depth: max_depth = depth |
---|
1212 | |
---|
1213 | |
---|
1214 | s = '%.2f, %.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, thisgauge[0], leg_label[j]) |
---|
1215 | fid_compare.write(s) |
---|
1216 | max_stages.append(max_stage) |
---|
1217 | min_stages.append(min_stage) |
---|
1218 | max_momentums.append(max_momentum) |
---|
1219 | max_xmomentums.append(max_xmomentum) |
---|
1220 | max_ymomentums.append(max_ymomentum) |
---|
1221 | min_xmomentums.append(min_xmomentum) |
---|
1222 | min_ymomentums.append(min_ymomentum) |
---|
1223 | min_momentums.append(min_momentum) |
---|
1224 | max_depths.append(max_depth) |
---|
1225 | max_speeds.append(max_speed) |
---|
1226 | min_speeds.append(min_speed) |
---|
1227 | #### finished generating quantities for each swwfile ##### |
---|
1228 | |
---|
1229 | model_time_plot3d[:,:] = model_time[:,:,j] |
---|
1230 | stages_plot3d[:,:] = stages[:,:,j] |
---|
1231 | eastings_plot3d[:,] = eastings[:,:,j] |
---|
1232 | |
---|
1233 | if surface is True: |
---|
1234 | print 'Printing surface figure' |
---|
1235 | for i in range(2): |
---|
1236 | fig = p1.figure(10) |
---|
1237 | ax = p3.Axes3D(fig) |
---|
1238 | if len(gauges) > 80: |
---|
1239 | ax.plot_surface(model_time[:,:,j],eastings[:,:,j],stages[:,:,j]) |
---|
1240 | else: |
---|
1241 | #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j]) |
---|
1242 | ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j])) |
---|
1243 | ax.set_xlabel('time') |
---|
1244 | ax.set_ylabel('x') |
---|
1245 | ax.set_zlabel('stage') |
---|
1246 | fig.add_axes(ax) |
---|
1247 | p1.show() |
---|
1248 | surfacefig = 'solution_surface%s' %leg_label[j] |
---|
1249 | p1.savefig(surfacefig) |
---|
1250 | p1.close() |
---|
1251 | |
---|
1252 | #### finished generating quantities for all swwfiles ##### |
---|
1253 | |
---|
1254 | # x profile for given time |
---|
1255 | if surface is True: |
---|
1256 | figure(11) |
---|
1257 | plot(eastings[tindex,:,j],stages[tindex,:,j]) |
---|
1258 | xlabel('x') |
---|
1259 | ylabel('stage') |
---|
1260 | profilefig = 'solution_xprofile' |
---|
1261 | savefig('profilefig') |
---|
1262 | |
---|
1263 | elev_output = [] |
---|
1264 | if generate_fig is True: |
---|
1265 | depth_axis = axis([starttime/scale, time_max/scale, -0.1, max(max_depths)*1.1]) |
---|
1266 | stage_axis = axis([starttime/scale, time_max/scale, min(min_stages), max(max_stages)*1.1]) |
---|
1267 | vel_axis = axis([starttime/scale, time_max/scale, min(min_speeds), max(max_speeds)*1.1]) |
---|
1268 | mom_axis = axis([starttime/scale, time_max/scale, min(min_momentums), max(max_momentums)*1.1]) |
---|
1269 | xmom_axis = axis([starttime/scale, time_max/scale, min(min_xmomentums), max(max_xmomentums)*1.1]) |
---|
1270 | ymom_axis = axis([starttime/scale, time_max/scale, min(min_ymomentums), max(max_ymomentums)*1.1]) |
---|
1271 | cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] |
---|
1272 | nn = len(plot_quantity) |
---|
1273 | no_cols = 2 |
---|
1274 | |
---|
1275 | if len(label_id) > 1: graphname_report = [] |
---|
1276 | pp = 1 |
---|
1277 | div = 11. |
---|
1278 | cc = 0 |
---|
1279 | for k in gauge_index: |
---|
1280 | g = gauges[k] |
---|
1281 | count1 = 0 |
---|
1282 | if report == True and len(label_id) > 1: |
---|
1283 | s = '\\begin{figure}[ht] \n \\centering \n \\begin{tabular}{cc} \n' |
---|
1284 | fid.write(s) |
---|
1285 | if len(label_id) > 1: graphname_report = [] |
---|
1286 | #### generate figures for each gauge #### |
---|
1287 | for j, f in enumerate(f_list): |
---|
1288 | ion() |
---|
1289 | hold(True) |
---|
1290 | count = 0 |
---|
1291 | where1 = 0 |
---|
1292 | where2 = 0 |
---|
1293 | word_quantity = '' |
---|
1294 | if report == True and len(label_id) == 1: |
---|
1295 | s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n' |
---|
1296 | fid.write(s) |
---|
1297 | |
---|
1298 | for which_quantity in plot_quantity: |
---|
1299 | count += 1 |
---|
1300 | where1 += 1 |
---|
1301 | figure(count, frameon = False) |
---|
1302 | if which_quantity == 'depth': |
---|
1303 | plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1304 | units = 'm' |
---|
1305 | axis(depth_axis) |
---|
1306 | if which_quantity == 'stage': |
---|
1307 | if elevations[0,k,j] <= 0: |
---|
1308 | plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1309 | axis(stage_axis) |
---|
1310 | else: |
---|
1311 | plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1312 | #axis(depth_axis) |
---|
1313 | units = 'm' |
---|
1314 | if which_quantity == 'momentum': |
---|
1315 | plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1316 | axis(mom_axis) |
---|
1317 | units = 'm^2 / sec' |
---|
1318 | if which_quantity == 'xmomentum': |
---|
1319 | plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1320 | axis(xmom_axis) |
---|
1321 | units = 'm^2 / sec' |
---|
1322 | if which_quantity == 'ymomentum': |
---|
1323 | plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1324 | axis(ymom_axis) |
---|
1325 | units = 'm^2 / sec' |
---|
1326 | if which_quantity == 'speed': |
---|
1327 | plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
1328 | axis(vel_axis) |
---|
1329 | units = 'm / sec' |
---|
1330 | if which_quantity == 'bearing': |
---|
1331 | due_east = 90.0*ones(shape(model_time[0:n[j]-1,k,j],Float)) |
---|
1332 | due_west = 270.0*ones(shape(model_time[0:n[j]-1,k,j],Float)) |
---|
1333 | plot(model_time[0:n[j]-1,k,j], bearings, '-', |
---|
1334 | model_time[0:n[j]-1,k,j], due_west, '-.', |
---|
1335 | model_time[0:n[j]-1,k,j], due_east, '-.') |
---|
1336 | units = 'degrees from North' |
---|
1337 | ax = axis([time_min, time_max, 0.0, 360.0]) |
---|
1338 | legend(('Bearing','West','East')) |
---|
1339 | |
---|
1340 | if time_unit is 'mins': xlabel('time (mins)') |
---|
1341 | if time_unit is 'hours': xlabel('time (hours)') |
---|
1342 | #if which_quantity == 'stage' and elevations[0:n[j]-1,k,j] > 0: |
---|
1343 | # ylabel('%s (%s)' %('depth', units)) |
---|
1344 | #else: |
---|
1345 | # ylabel('%s (%s)' %(which_quantity, units)) |
---|
1346 | #ylabel('%s (%s)' %('wave height', units)) |
---|
1347 | ylabel('%s (%s)' %(which_quantity, units)) |
---|
1348 | if len(label_id) > 1: legend((leg_label),loc='upper right') |
---|
1349 | |
---|
1350 | #gaugeloc1 = gaugeloc.replace(' ','') |
---|
1351 | #gaugeloc2 = gaugeloc1.replace('_','') |
---|
1352 | gaugeloc2 = str(locations[k]).replace(' ','') |
---|
1353 | graphname = '%sgauge%s_%s' %(file_loc[j], |
---|
1354 | gaugeloc2, |
---|
1355 | which_quantity) |
---|
1356 | |
---|
1357 | if report == True and len(label_id) > 1: |
---|
1358 | figdir = getcwd()+sep+'report_figures'+sep |
---|
1359 | if access(figdir,F_OK) == 0 : |
---|
1360 | mkdir (figdir) |
---|
1361 | latex_file_loc = figdir.replace(sep,altsep) |
---|
1362 | graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory |
---|
1363 | graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file |
---|
1364 | graphname_report.append(graphname_report_input) |
---|
1365 | |
---|
1366 | savefig(graphname_latex) # save figures in production directory for report generation |
---|
1367 | |
---|
1368 | if report == True: |
---|
1369 | |
---|
1370 | figdir = getcwd()+sep+'report_figures'+sep |
---|
1371 | if access(figdir,F_OK) == 0 : |
---|
1372 | mkdir (figdir) |
---|
1373 | latex_file_loc = figdir.replace(sep,altsep) |
---|
1374 | |
---|
1375 | if len(label_id) == 1: |
---|
1376 | graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory |
---|
1377 | graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file |
---|
1378 | s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png') |
---|
1379 | fid.write(s) |
---|
1380 | if where1 % 2 == 0: |
---|
1381 | s = '\\\\ \n' |
---|
1382 | where1 = 0 |
---|
1383 | else: |
---|
1384 | s = '& \n' |
---|
1385 | fid.write(s) |
---|
1386 | savefig(graphname_latex) |
---|
1387 | |
---|
1388 | if title_on == True: |
---|
1389 | title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2)) |
---|
1390 | #title('Gauge %s (MOST elevation %.2f, ANUGA elevation %.2f)' %(gaugeloc2, elevations[10,k,0], elevations[10,k,1] )) |
---|
1391 | |
---|
1392 | savefig(graphname) # save figures with sww file |
---|
1393 | |
---|
1394 | if report == True and len(label_id) == 1: |
---|
1395 | for i in range(nn-1): |
---|
1396 | if nn > 2: |
---|
1397 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1398 | word_quantity += 'depth' + ', ' |
---|
1399 | else: |
---|
1400 | word_quantity += plot_quantity[i] + ', ' |
---|
1401 | else: |
---|
1402 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1403 | word_quantity += 'depth' + ', ' |
---|
1404 | else: |
---|
1405 | word_quantity += plot_quantity[i] |
---|
1406 | |
---|
1407 | if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0: |
---|
1408 | word_quantity += ' and ' + 'depth' |
---|
1409 | else: |
---|
1410 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
1411 | caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' ')) |
---|
1412 | if elev[k] == 0.0: |
---|
1413 | caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j]) |
---|
1414 | east = gauges[0] |
---|
1415 | north = gauges[1] |
---|
1416 | elev_output.append([locations[k],east,north,elevations[0,k,j]]) |
---|
1417 | label = '%sgauge%s' %(label_id2, gaugeloc2) |
---|
1418 | s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label) |
---|
1419 | fid.write(s) |
---|
1420 | cc += 1 |
---|
1421 | if cc % 6 == 0: fid.write('\\clearpage \n') |
---|
1422 | savefig(graphname_latex) |
---|
1423 | |
---|
1424 | if report == True and len(label_id) > 1: |
---|
1425 | for i in range(nn-1): |
---|
1426 | if nn > 2: |
---|
1427 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1428 | word_quantity += 'depth' + ',' |
---|
1429 | else: |
---|
1430 | word_quantity += plot_quantity[i] + ', ' |
---|
1431 | else: |
---|
1432 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
1433 | word_quantity += 'depth' |
---|
1434 | else: |
---|
1435 | word_quantity += plot_quantity[i] |
---|
1436 | where1 = 0 |
---|
1437 | count1 += 1 |
---|
1438 | index = j*len(plot_quantity) |
---|
1439 | for which_quantity in plot_quantity: |
---|
1440 | where1 += 1 |
---|
1441 | s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report[index], '.png') |
---|
1442 | index += 1 |
---|
1443 | fid.write(s) |
---|
1444 | if where1 % 2 == 0: |
---|
1445 | s = '\\\\ \n' |
---|
1446 | where1 = 0 |
---|
1447 | else: |
---|
1448 | s = '& \n' |
---|
1449 | fid.write(s) |
---|
1450 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
1451 | label = 'gauge%s' %(gaugeloc2) |
---|
1452 | caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) |
---|
1453 | if elev[k] == 0.0: |
---|
1454 | caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j]) |
---|
1455 | thisgauge = gauges[k] |
---|
1456 | east = thisgauge[0] |
---|
1457 | north = thisgauge[1] |
---|
1458 | elev_output.append([locations[k],east,north,elevations[0,k,j]]) |
---|
1459 | |
---|
1460 | s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label) |
---|
1461 | fid.write(s) |
---|
1462 | if float((k+1)/div - pp) == 0.: |
---|
1463 | fid.write('\\clearpage \n') |
---|
1464 | pp += 1 |
---|
1465 | |
---|
1466 | #### finished generating figures ### |
---|
1467 | |
---|
1468 | close('all') |
---|
1469 | |
---|
1470 | return texfile2, elev_output |
---|
1471 | |
---|
1472 | # FIXME (DSG): Add unit test, make general, not just 2 files, |
---|
1473 | # but any number of files. |
---|
1474 | def copy_code_files(dir_name, filename1, filename2): |
---|
1475 | """Temporary Interface to new location""" |
---|
1476 | |
---|
1477 | from anuga.shallow_water.data_manager import \ |
---|
1478 | copy_code_files as dm_copy_code_files |
---|
1479 | print 'copy_code_files has moved from util.py. ', |
---|
1480 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1481 | import copy_code_files"' |
---|
1482 | |
---|
1483 | return dm_copy_code_files(dir_name, filename1, filename2) |
---|
1484 | |
---|
1485 | |
---|
1486 | def add_directories(root_directory, directories): |
---|
1487 | """ |
---|
1488 | Add the first directory in directories to root_directory. |
---|
1489 | Then add the second |
---|
1490 | directory to the first directory and so on. |
---|
1491 | |
---|
1492 | Return the path of the final directory. |
---|
1493 | |
---|
1494 | This is handy for specifying and creating a directory |
---|
1495 | where data will go. |
---|
1496 | """ |
---|
1497 | dir = root_directory |
---|
1498 | for new_dir in directories: |
---|
1499 | dir = os.path.join(dir, new_dir) |
---|
1500 | if not access(dir,F_OK): |
---|
1501 | mkdir (dir) |
---|
1502 | return dir |
---|
1503 | |
---|
1504 | def get_data_from_file(filename,separator_value = ','): |
---|
1505 | """Temporary Interface to new location""" |
---|
1506 | from anuga.shallow_water.data_manager import \ |
---|
1507 | get_data_from_file as dm_get_data_from_file |
---|
1508 | print 'get_data_from_file has moved from util.py' |
---|
1509 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1510 | import get_data_from_file"' |
---|
1511 | |
---|
1512 | return dm_get_data_from_file(filename,separator_value = ',') |
---|
1513 | |
---|
1514 | def store_parameters(verbose=False,**kwargs): |
---|
1515 | """Temporary Interface to new location""" |
---|
1516 | |
---|
1517 | from anuga.shallow_water.data_manager \ |
---|
1518 | import store_parameters as dm_store_parameters |
---|
1519 | print 'store_parameters has moved from util.py.' |
---|
1520 | print 'Please use "from anuga.shallow_water.data_manager \ |
---|
1521 | import store_parameters"' |
---|
1522 | |
---|
1523 | return dm_store_parameters(verbose=False,**kwargs) |
---|
1524 | |
---|
1525 | def remove_lone_verts(verts, triangles, number_of_full_nodes=None): |
---|
1526 | """ |
---|
1527 | Removes vertices that are not associated with any triangles. |
---|
1528 | |
---|
1529 | verts is a list/array of points |
---|
1530 | triangles is a list of 3 element tuples. |
---|
1531 | Each tuple represents a triangle. |
---|
1532 | |
---|
1533 | number_of_full_nodes relate to parallelism when a mesh has an |
---|
1534 | extra layer of ghost points. |
---|
1535 | |
---|
1536 | """ |
---|
1537 | verts = ensure_numeric(verts) |
---|
1538 | triangles = ensure_numeric(triangles) |
---|
1539 | |
---|
1540 | N = len(verts) |
---|
1541 | |
---|
1542 | # initialise the array to easily find the index of the first loner |
---|
1543 | loners=arange(2*N, N, -1) # if N=3 [6,5,4] |
---|
1544 | for t in triangles: |
---|
1545 | for vert in t: |
---|
1546 | loners[vert]= vert # all non-loners will have loners[i]=i |
---|
1547 | #print loners |
---|
1548 | |
---|
1549 | lone_start = 2*N - max(loners) # The index of the first loner |
---|
1550 | |
---|
1551 | if lone_start-1 == N: |
---|
1552 | # no loners |
---|
1553 | pass |
---|
1554 | elif min(loners[lone_start:N]) > N: |
---|
1555 | # All the loners are at the end of the vert array |
---|
1556 | verts = verts[0:lone_start] |
---|
1557 | else: |
---|
1558 | # change the loners list so it can be used to modify triangles |
---|
1559 | # Remove the loners from verts |
---|
1560 | # Could've used X=compress(less(loners,N),loners) |
---|
1561 | # verts=take(verts,X) to Remove the loners from verts |
---|
1562 | # but I think it would use more memory |
---|
1563 | new_i = 0 |
---|
1564 | for i in range(lone_start, N): |
---|
1565 | if loners[i] >= N: |
---|
1566 | # Loner! |
---|
1567 | pass |
---|
1568 | else: |
---|
1569 | loners[i] = new_i |
---|
1570 | verts[new_i] = verts[i] |
---|
1571 | new_i += 1 |
---|
1572 | verts = verts[0:new_i] |
---|
1573 | |
---|
1574 | # Modify the triangles |
---|
1575 | #print "loners", loners |
---|
1576 | #print "triangles before", triangles |
---|
1577 | triangles = choose(triangles,loners) |
---|
1578 | #print "triangles after", triangles |
---|
1579 | return verts, triangles |
---|
1580 | |
---|
1581 | |
---|
1582 | def get_centroid_values(x, triangles): |
---|
1583 | """Compute centroid values from vertex values |
---|
1584 | |
---|
1585 | x: Values at vertices of triangular mesh |
---|
1586 | triangles: Nx3 integer array pointing to vertex information |
---|
1587 | for each of the N triangels. Elements of triangles are |
---|
1588 | indices into x |
---|
1589 | """ |
---|
1590 | |
---|
1591 | |
---|
1592 | xc = zeros(triangles.shape[0], Float) # Space for centroid info |
---|
1593 | |
---|
1594 | for k in range(triangles.shape[0]): |
---|
1595 | # Indices of vertices |
---|
1596 | i0 = triangles[k][0] |
---|
1597 | i1 = triangles[k][1] |
---|
1598 | i2 = triangles[k][2] |
---|
1599 | |
---|
1600 | xc[k] = (x[i0] + x[i1] + x[i2])/3 |
---|
1601 | |
---|
1602 | |
---|
1603 | return xc |
---|
1604 | |
---|
1605 | def make_plots_from_csv_file(directories_dic={}, |
---|
1606 | output_dir='', |
---|
1607 | base_name=None, |
---|
1608 | plot_numbers='0', |
---|
1609 | quantities=['Stage'], |
---|
1610 | extra_plot_name='', |
---|
1611 | assess_all_csv_files=True, |
---|
1612 | create_latex=False, |
---|
1613 | verbose=False): |
---|
1614 | |
---|
1615 | """Read in csv files that have the right header information and |
---|
1616 | plot time series such as Stage, Speed, etc. Will also plot several |
---|
1617 | time series on one plot. Filenames must follow this convention, |
---|
1618 | <base_name><plot_number>.csv eg gauge_timeseries3.csv |
---|
1619 | |
---|
1620 | For example: |
---|
1621 | if "directories_dic" defines 4 directories and in each directories |
---|
1622 | there is a csv files corresponding to the right "plot_numbers", |
---|
1623 | this will create a plot with 4 lines one for each directory AND |
---|
1624 | one plot for each "quantities". |
---|
1625 | |
---|
1626 | Inputs: |
---|
1627 | directories_dic: dictionary of directory with values (plot |
---|
1628 | legend name for directory), (start time of |
---|
1629 | the time series) and the (value to add to |
---|
1630 | stage if needed). For example |
---|
1631 | {dir1:['Anuga_ons',5000, 0], |
---|
1632 | dir2:['b_emoth',5000,1.5], |
---|
1633 | dir3:['b_ons',5000,1.5]} |
---|
1634 | |
---|
1635 | output_dir: directory for the plot outputs |
---|
1636 | |
---|
1637 | base_name: common name to all the csv files to be read |
---|
1638 | |
---|
1639 | plot_numbers: a String list of numbers to plot. For example |
---|
1640 | [0-4,10,15-17] will read and attempt to plot |
---|
1641 | the follow 0,1,2,3,4,10,15,16,17 |
---|
1642 | quantities: Currently limited to "Stage", "Speed", and |
---|
1643 | "Momentum", should be changed to incorporate |
---|
1644 | any quantity read from CSV file.... |
---|
1645 | |
---|
1646 | extra_plot_name: A string that is appended to the end of the |
---|
1647 | output filename. |
---|
1648 | |
---|
1649 | assess_all_csv_files: if true it will read ALL csv file with |
---|
1650 | "base_name", regardless of 'plot_numbers' |
---|
1651 | and determine a uniform set of axes for |
---|
1652 | Stage, Speed and Momentum. IF FALSE it |
---|
1653 | will only read the csv file within the |
---|
1654 | 'plot_numbers' |
---|
1655 | |
---|
1656 | create_latex: NOT IMPLEMENTED YET!! sorry Jane.... |
---|
1657 | |
---|
1658 | OUTPUTS: None, it saves the plots to |
---|
1659 | <output_dir><base_name><plot_number><extra_plot_name>.png |
---|
1660 | |
---|
1661 | """ |
---|
1662 | import pylab# import plot, xlabel, ylabel, savefig, \ |
---|
1663 | #ion, hold, axis, close, figure, legend |
---|
1664 | from os import sep |
---|
1665 | from anuga.shallow_water.data_manager import \ |
---|
1666 | get_all_files_with_extension, csv2dict |
---|
1667 | #find all the files that meet the specs |
---|
1668 | #FIXME remove all references to word that such as Stage |
---|
1669 | #(with a Capital letter) could use the header info from |
---|
1670 | #the dictionary that is returned from the csv2dict this could |
---|
1671 | #be very good and very flexable.... but this works for now! |
---|
1672 | |
---|
1673 | #FIXME plot all available data from the csv file, currently will |
---|
1674 | #only plot Stage,Speed and Momentum and only if the header is |
---|
1675 | #named like wise |
---|
1676 | |
---|
1677 | # quantities_dic={'time':'time (hour)'} |
---|
1678 | # if 'time' in quantities: |
---|
1679 | # quantities_dic['time']='time (hour)' |
---|
1680 | # if 'Time' in quantities: |
---|
1681 | # quantities_dic['Time']='time (hour)' |
---|
1682 | # if 'stage' in quantities: |
---|
1683 | # quantities_dic['stage']='wave height (m)' |
---|
1684 | # if 'Stage' in quantities: |
---|
1685 | # quantities_dic['Stage']='wave height (m)' |
---|
1686 | # if 'speed' in quantities: |
---|
1687 | # quantities_dic['speed']='speed (m/s)' |
---|
1688 | # if 'Speed' in quantities: |
---|
1689 | # quantities_dic['Speed']='speed (m/s)' |
---|
1690 | # if 'stage' or 'Stage' in quantities: |
---|
1691 | # quantities_dic['stage']='speed (m/s)' |
---|
1692 | # if 'stage' or 'Stage' in quantities: |
---|
1693 | # quantities_dic['stage']='speed (m/s)' |
---|
1694 | |
---|
1695 | seconds_in_hour = 3600 |
---|
1696 | time_label = 'time (hour)' |
---|
1697 | stage_label = 'wave height (m)' |
---|
1698 | speed_label = 'speed (m/s)' |
---|
1699 | momentum_label = 'momentum (m^2/sec)' |
---|
1700 | |
---|
1701 | if extra_plot_name != '': |
---|
1702 | extra_plot_name='_'+extra_plot_name |
---|
1703 | |
---|
1704 | #finds all the files that fit the specs and return a list of them |
---|
1705 | #so to help find a uniform max and min for the plots... |
---|
1706 | list_filenames=[] |
---|
1707 | if verbose: print 'Determining files to access for axes ranges \n' |
---|
1708 | for i,directory in enumerate(directories_dic.keys()): |
---|
1709 | list_filenames.append(get_all_files_with_extension(directory, |
---|
1710 | base_name,'.csv')) |
---|
1711 | # print 'list_filenames',list_filenames |
---|
1712 | |
---|
1713 | #use all the files to get the values for the plot axis |
---|
1714 | max_st=max_sp=max_mom=min_st=min_sp=min_mom=max_t=min_t=0. |
---|
1715 | max_start_time= 0. |
---|
1716 | min_start_time = 100000 |
---|
1717 | |
---|
1718 | |
---|
1719 | new_plot_numbers=[] |
---|
1720 | #change plot_numbers to list, eg ['0-4','10'] |
---|
1721 | #to ['0','1','2','3','4','10'] |
---|
1722 | for i, num_string in enumerate(plot_numbers): |
---|
1723 | if '-' in num_string: |
---|
1724 | start = int(num_string[:num_string.rfind('-')]) |
---|
1725 | end = int(num_string[num_string.rfind('-')+1:])+1 |
---|
1726 | for x in range(start, end): |
---|
1727 | new_plot_numbers.append(str(x)) |
---|
1728 | else: |
---|
1729 | new_plot_numbers.append(num_string) |
---|
1730 | # print 'new_plot_numbers',new_plot_numbers |
---|
1731 | |
---|
1732 | if verbose: print 'Determining uniform axes \n' |
---|
1733 | #this entire loop is to determine the min and max range for the |
---|
1734 | #axes of the plot |
---|
1735 | for i, directory in enumerate(directories_dic.keys()): |
---|
1736 | |
---|
1737 | if assess_all_csv_files==False: |
---|
1738 | which_csv_to_assess = new_plot_numbers |
---|
1739 | else: |
---|
1740 | which_csv_to_assess = list_filenames[i] |
---|
1741 | |
---|
1742 | for j, filename in enumerate(which_csv_to_assess): |
---|
1743 | if assess_all_csv_files==False: |
---|
1744 | dir_filename=directory+base_name+sep+filename |
---|
1745 | else: |
---|
1746 | dir_filename=join(directory,filename) |
---|
1747 | |
---|
1748 | attribute_dic, title_index_dic = csv2dict(dir_filename+ |
---|
1749 | '.csv') |
---|
1750 | |
---|
1751 | directory_start_time = directories_dic[directory][1] |
---|
1752 | directory_add_tide = directories_dic[directory][2] |
---|
1753 | |
---|
1754 | time = [float(x) for x in attribute_dic["Time"]] |
---|
1755 | min_t, max_t = get_min_max_values(time,min_t,max_t) |
---|
1756 | |
---|
1757 | stage = [float(x) for x in attribute_dic["Stage"]] |
---|
1758 | stage =array(stage)+directory_add_tide |
---|
1759 | min_st, max_st = get_min_max_values(stage,min_st,max_st) |
---|
1760 | |
---|
1761 | speed = [float(x) for x in attribute_dic["Speed"]] |
---|
1762 | min_sp, max_sp = get_min_max_values(speed,min_sp,max_sp) |
---|
1763 | |
---|
1764 | momentum = [float(x) for x in attribute_dic["Momentum"]] |
---|
1765 | min_mom, max_mom = get_min_max_values(momentum, |
---|
1766 | min_mom, |
---|
1767 | max_mom) |
---|
1768 | |
---|
1769 | # print 'min_sp, max_sp',min_sp, max_sp, |
---|
1770 | # print directory_start_time |
---|
1771 | if min_start_time > directory_start_time: |
---|
1772 | min_start_time = directory_start_time |
---|
1773 | if max_start_time < directory_start_time: |
---|
1774 | max_start_time = directory_start_time |
---|
1775 | # print 'start_time' , max_start_time, min_start_time |
---|
1776 | |
---|
1777 | stage_axis = (min_start_time/seconds_in_hour, |
---|
1778 | (max_t+max_start_time)/seconds_in_hour, |
---|
1779 | min_st, max_st) |
---|
1780 | speed_axis = (min_start_time/seconds_in_hour, |
---|
1781 | (max_t+max_start_time)/seconds_in_hour, |
---|
1782 | min_sp, max_sp) |
---|
1783 | momentum_axis = (min_start_time/seconds_in_hour, |
---|
1784 | (max_t+max_start_time)/seconds_in_hour, |
---|
1785 | min_mom, max_mom) |
---|
1786 | # ion() |
---|
1787 | pylab.hold() |
---|
1788 | |
---|
1789 | |
---|
1790 | cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k'] |
---|
1791 | |
---|
1792 | # if verbose: print 'Now start to plot \n' |
---|
1793 | |
---|
1794 | i_max = len(directories_dic.keys()) |
---|
1795 | legend_list_dic =[] |
---|
1796 | legend_list =[] |
---|
1797 | for i, directory in enumerate(directories_dic.keys()): |
---|
1798 | if verbose: print'Plotting in %s', directory |
---|
1799 | for j, number in enumerate(new_plot_numbers): |
---|
1800 | if verbose: print'Starting %s',base_name+number |
---|
1801 | directory_name = directories_dic[directory][0] |
---|
1802 | directory_start_time = directories_dic[directory][1] |
---|
1803 | directory_add_tide = directories_dic[directory][2] |
---|
1804 | |
---|
1805 | #create an if about the start time and tide hieght |
---|
1806 | #if they don't exist |
---|
1807 | file=directory+sep+base_name+number |
---|
1808 | #print 'i %s,j %s, number %s, file %s' %(i,j,number,file) |
---|
1809 | attribute_dic, title_index_dic = csv2dict(file+'.csv') |
---|
1810 | #get data from dict in to list |
---|
1811 | t = [float(x) for x in attribute_dic["Time"]] |
---|
1812 | |
---|
1813 | #do maths to list by changing to array |
---|
1814 | t=(array(t)+directory_start_time)/seconds_in_hour |
---|
1815 | stage = [float(x) for x in attribute_dic["Stage"]] |
---|
1816 | speed = [float(x) for x in attribute_dic["Speed"]] |
---|
1817 | momentum = [float(x) for x in attribute_dic["Momentum"]] |
---|
1818 | |
---|
1819 | #Can add tide so plots are all on the same tide height |
---|
1820 | stage =array(stage)+directory_add_tide |
---|
1821 | |
---|
1822 | #finds the maximum elevation |
---|
1823 | max_ele=-100000 |
---|
1824 | min_ele=100000 |
---|
1825 | elevation = [float(x) for x in attribute_dic["Elevation"]] |
---|
1826 | |
---|
1827 | min_ele, max_ele = get_min_max_values(elevation, |
---|
1828 | min_ele, |
---|
1829 | max_ele) |
---|
1830 | if min_ele != max_ele: |
---|
1831 | print "Note! Elevation changes in %s" %dir_filename |
---|
1832 | # print 'min_ele, max_ele',min_ele, max_ele |
---|
1833 | |
---|
1834 | |
---|
1835 | #populates the legend_list_dic with dir_name and the elevation |
---|
1836 | if i==0: |
---|
1837 | legend_list_dic.append({directory_name:max_ele}) |
---|
1838 | else: |
---|
1839 | legend_list_dic[j][directory_name]=max_ele |
---|
1840 | |
---|
1841 | # creates a list for the legend after at "legend_dic" has been fully populated |
---|
1842 | # only runs on the last iteration for all the gauges(csv) files |
---|
1843 | # empties the list before creating it |
---|
1844 | if i==i_max-1: |
---|
1845 | legend_list=[] |
---|
1846 | for k, l in legend_list_dic[j].iteritems(): |
---|
1847 | legend_list.append('%s (elevation = %sm)'%(k,l)) |
---|
1848 | #print k,l, legend_list_dic[j] |
---|
1849 | |
---|
1850 | if "Stage" in quantities: |
---|
1851 | pylab.figure(100+j) |
---|
1852 | pylab.plot(t, stage, c = cstr[i], linewidth=1) |
---|
1853 | pylab.xlabel(time_label) |
---|
1854 | pylab.ylabel(stage_label) |
---|
1855 | pylab.axis(stage_axis) |
---|
1856 | pylab.legend(legend_list,loc='upper right') |
---|
1857 | figname = '%sstage_%s%s.png' %(output_dir+sep, |
---|
1858 | base_name+number, |
---|
1859 | extra_plot_name) |
---|
1860 | pylab.savefig(figname) |
---|
1861 | if "Speed" in quantities: |
---|
1862 | pylab.figure(200+j) |
---|
1863 | pylab.plot(t, speed, c = cstr[i], linewidth=1) |
---|
1864 | pylab.xlabel(time_label) |
---|
1865 | pylab.ylabel(speed_label) |
---|
1866 | pylab.axis(speed_axis) |
---|
1867 | pylab.legend(legend_list,loc='upper right') |
---|
1868 | figname = '%sspeed_%s%s.png' %(output_dir+sep, |
---|
1869 | base_name+number, |
---|
1870 | extra_plot_name) |
---|
1871 | pylab.savefig(figname) |
---|
1872 | if "Momentum" in quantities: |
---|
1873 | pylab.figure(300+j) |
---|
1874 | pylab.plot(t, momentum, c = cstr[i], linewidth=1) |
---|
1875 | pylab.xlabel(time_label) |
---|
1876 | pylab.ylabel(momentum_label) |
---|
1877 | pylab.axis(momentum_axis) |
---|
1878 | pylab.legend(legend_list,loc='upper right') |
---|
1879 | figname = '%smomentum_%s%s.png' %(output_dir+sep, |
---|
1880 | base_name+number, |
---|
1881 | extra_plot_name) |
---|
1882 | pylab.savefig(figname) |
---|
1883 | if verbose: print 'Closing all plots' |
---|
1884 | pylab.close('all') |
---|
1885 | del pylab |
---|
1886 | if verbose: print 'Finished closing plots' |
---|
1887 | |
---|
1888 | def get_min_max_values(list=None,min1=100,max1=-100): |
---|
1889 | """ Returns the min and max of the list it was provided. |
---|
1890 | NOTE: default min and max may need to change depeending on |
---|
1891 | your list |
---|
1892 | """ |
---|
1893 | if list == None: print 'List must be provided' |
---|
1894 | # min = max_list = 0 |
---|
1895 | if max(list) > max1: |
---|
1896 | max1 = max(list) |
---|
1897 | if min(list) < min1: |
---|
1898 | min1 = min(list) |
---|
1899 | |
---|
1900 | return min1, max1 |
---|
1901 | |
---|
1902 | |
---|
1903 | |
---|
1904 | |
---|
1905 | |
---|
1906 | |
---|
1907 | |
---|