1 | """Gauge functions |
---|
2 | |
---|
3 | Functions for converting gauge and sww files into timeseries plots |
---|
4 | |
---|
5 | |
---|
6 | Copyright 2010 |
---|
7 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou, James Hudson |
---|
8 | Geoscience Australia |
---|
9 | """ |
---|
10 | |
---|
11 | import numpy as num |
---|
12 | |
---|
13 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
14 | |
---|
15 | import os |
---|
16 | |
---|
17 | from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep, getcwd |
---|
18 | from os.path import exists, split, join |
---|
19 | import anuga.utilities.log as log |
---|
20 | |
---|
21 | from math import sqrt |
---|
22 | |
---|
23 | ## |
---|
24 | # @brief Take gauge readings, given a gauge file and a sww mesh |
---|
25 | # |
---|
26 | # Use this function to take a timeseries sample, given a list of gauge points |
---|
27 | # @param sww_file sww file to use as input |
---|
28 | # @param gauge_file gauge file as input, containing list of gauge points to sample |
---|
29 | # @param out_name output file prefix |
---|
30 | # @param quantities which quantities in the sww file we want to export |
---|
31 | # @param verbose show extra logging information for debug purposes |
---|
32 | # @param use_cache cache requests if possible, for speed |
---|
33 | # @param output_centroids Set to true to output the values at the centroid of the mesh triangle |
---|
34 | def gauge_sww2csv(sww_file, |
---|
35 | gauge_file, |
---|
36 | out_name='gauge_', |
---|
37 | quantities=['stage', 'depth', 'elevation', |
---|
38 | 'xmomentum', 'ymomentum'], |
---|
39 | verbose=False, |
---|
40 | use_cache=True, |
---|
41 | output_centroids=False): |
---|
42 | """ |
---|
43 | |
---|
44 | Inputs: |
---|
45 | NOTE: if using csv2timeseries_graphs after creating csv file, |
---|
46 | it is essential to export quantities 'depth' and 'elevation'. |
---|
47 | 'depth' is good to analyse gauges on land and elevation is used |
---|
48 | automatically by csv2timeseries_graphs in the legend. |
---|
49 | |
---|
50 | sww_file: path to any sww file |
---|
51 | |
---|
52 | gauge_file: Assumes that it follows this format |
---|
53 | name, easting, northing, elevation |
---|
54 | point1, 100.3, 50.2, 10.0 |
---|
55 | point2, 10.3, 70.3, 78.0 |
---|
56 | |
---|
57 | NOTE: order of column can change but names eg 'easting', 'elevation' |
---|
58 | must be the same! ALL lowercaps! |
---|
59 | |
---|
60 | out_name: prefix for output file name (default is 'gauge_') |
---|
61 | |
---|
62 | Outputs: |
---|
63 | one file for each gauge/point location in the points file. They |
---|
64 | will be named with this format in the same directory as the 'sww_file' |
---|
65 | <out_name><name>.csv |
---|
66 | eg gauge_point1.csv if <out_name> not supplied |
---|
67 | myfile_2_point1.csv if <out_name> ='myfile_2_' |
---|
68 | |
---|
69 | They will all have a header |
---|
70 | |
---|
71 | Usage: sww2csv_gauges(sww_file='test1.sww', |
---|
72 | quantities = ['stage', 'elevation','depth','bearing'], |
---|
73 | gauge_file='gauge.txt') |
---|
74 | |
---|
75 | Interpolate the quantities at a given set of locations, given |
---|
76 | an sww file. |
---|
77 | The results are written to a csv file. |
---|
78 | |
---|
79 | In the future let points be a points file. |
---|
80 | And the user choose the quantities. |
---|
81 | |
---|
82 | This is currently quite specific. |
---|
83 | If it needs to be more general, change things. |
---|
84 | |
---|
85 | This is really returning speed, not velocity. |
---|
86 | """ |
---|
87 | |
---|
88 | from csv import reader,writer |
---|
89 | from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN |
---|
90 | import string |
---|
91 | from anuga.shallow_water.data_manager import get_all_swwfiles |
---|
92 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
93 | |
---|
94 | assert type(gauge_file) == type(''), 'Gauge filename must be a string' |
---|
95 | assert type(out_name) == type(''), 'Output filename prefix must be a string' |
---|
96 | |
---|
97 | try: |
---|
98 | point_reader = reader(file(gauge_file)) |
---|
99 | except Exception, e: |
---|
100 | msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e) |
---|
101 | raise msg |
---|
102 | |
---|
103 | if verbose: log.critical('Gauges obtained from: %s' % gauge_file) |
---|
104 | |
---|
105 | point_reader = reader(file(gauge_file)) |
---|
106 | points = [] |
---|
107 | point_name = [] |
---|
108 | |
---|
109 | # read point info from file |
---|
110 | for i,row in enumerate(point_reader): |
---|
111 | # read header and determine the column numbers to read correctly. |
---|
112 | if i==0: |
---|
113 | for j,value in enumerate(row): |
---|
114 | if value.strip()=='easting':easting=j |
---|
115 | if value.strip()=='northing':northing=j |
---|
116 | if value.strip()=='name':name=j |
---|
117 | if value.strip()=='elevation':elevation=j |
---|
118 | else: |
---|
119 | points.append([float(row[easting]),float(row[northing])]) |
---|
120 | point_name.append(row[name]) |
---|
121 | |
---|
122 | #convert to array for file_function |
---|
123 | points_array = num.array(points,num.float) |
---|
124 | |
---|
125 | points_array = ensure_absolute(points_array) |
---|
126 | |
---|
127 | dir_name, base = os.path.split(sww_file) |
---|
128 | |
---|
129 | #need to get current directory so when path and file |
---|
130 | #are "joined" below the directory is correct |
---|
131 | if dir_name == '': |
---|
132 | dir_name =getcwd() |
---|
133 | |
---|
134 | if access(sww_file,R_OK): |
---|
135 | if verbose: log.critical('File %s exists' % sww_file) |
---|
136 | else: |
---|
137 | msg = 'File "%s" could not be opened: no read permission' % sww_file |
---|
138 | raise msg |
---|
139 | |
---|
140 | sww_files = get_all_swwfiles(look_in_dir=dir_name, |
---|
141 | base_name=base, |
---|
142 | verbose=verbose) |
---|
143 | |
---|
144 | # fudge to get SWW files in 'correct' order, oldest on the left |
---|
145 | sww_files.sort() |
---|
146 | |
---|
147 | if verbose: |
---|
148 | log.critical('sww files=%s' % sww_files) |
---|
149 | |
---|
150 | #to make all the quantities lower case for file_function |
---|
151 | quantities = [quantity.lower() for quantity in quantities] |
---|
152 | |
---|
153 | # what is quantities are needed from sww file to calculate output quantities |
---|
154 | # also |
---|
155 | |
---|
156 | core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
157 | |
---|
158 | gauge_file = out_name |
---|
159 | |
---|
160 | heading = [quantity for quantity in quantities] |
---|
161 | heading.insert(0,'time') |
---|
162 | heading.insert(1,'hours') |
---|
163 | |
---|
164 | #create a list of csv writers for all the points and write header |
---|
165 | points_writer = [] |
---|
166 | for point_i,point in enumerate(points): |
---|
167 | points_writer.append(writer(file(dir_name + sep + gauge_file |
---|
168 | + point_name[point_i] + '.csv', "wb"))) |
---|
169 | points_writer[point_i].writerow(heading) |
---|
170 | |
---|
171 | if verbose: log.critical('Writing csv files') |
---|
172 | |
---|
173 | quake_offset_time = None |
---|
174 | |
---|
175 | for sww_file in sww_files: |
---|
176 | sww_file = join(dir_name, sww_file+'.sww') |
---|
177 | callable_sww = file_function(sww_file, |
---|
178 | quantities=core_quantities, |
---|
179 | interpolation_points=points_array, |
---|
180 | verbose=verbose, |
---|
181 | use_cache=use_cache, |
---|
182 | output_centroids = output_centroids) |
---|
183 | |
---|
184 | if quake_offset_time is None: |
---|
185 | quake_offset_time = callable_sww.starttime |
---|
186 | |
---|
187 | for time in callable_sww.get_time(): |
---|
188 | for point_i, point in enumerate(points_array): |
---|
189 | #add domain starttime to relative time. |
---|
190 | quake_time = time + quake_offset_time |
---|
191 | points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug |
---|
192 | point_quantities = callable_sww(time, point_i) # __call__ is overridden |
---|
193 | |
---|
194 | for quantity in quantities: |
---|
195 | if quantity == NAN: |
---|
196 | log.critical('quantity does not exist in %s' |
---|
197 | % callable_sww.get_name) |
---|
198 | else: |
---|
199 | #core quantities that are exported from the interpolator |
---|
200 | if quantity == 'stage': |
---|
201 | points_list.append(point_quantities[0]) |
---|
202 | |
---|
203 | if quantity == 'elevation': |
---|
204 | points_list.append(point_quantities[1]) |
---|
205 | |
---|
206 | if quantity == 'xmomentum': |
---|
207 | points_list.append(point_quantities[2]) |
---|
208 | |
---|
209 | if quantity == 'ymomentum': |
---|
210 | points_list.append(point_quantities[3]) |
---|
211 | |
---|
212 | #derived quantities that are calculated from the core ones |
---|
213 | if quantity == 'depth': |
---|
214 | points_list.append(point_quantities[0] |
---|
215 | - point_quantities[1]) |
---|
216 | |
---|
217 | if quantity == 'momentum': |
---|
218 | momentum = sqrt(point_quantities[2]**2 |
---|
219 | + point_quantities[3]**2) |
---|
220 | points_list.append(momentum) |
---|
221 | |
---|
222 | if quantity == 'speed': |
---|
223 | #if depth is less than 0.001 then speed = 0.0 |
---|
224 | if point_quantities[0] - point_quantities[1] < 0.001: |
---|
225 | vel = 0.0 |
---|
226 | else: |
---|
227 | if point_quantities[2] < 1.0e6: |
---|
228 | momentum = sqrt(point_quantities[2]**2 |
---|
229 | + point_quantities[3]**2) |
---|
230 | vel = momentum / (point_quantities[0] |
---|
231 | - point_quantities[1]) |
---|
232 | else: |
---|
233 | momentum = 0 |
---|
234 | vel = 0 |
---|
235 | |
---|
236 | points_list.append(vel) |
---|
237 | |
---|
238 | if quantity == 'bearing': |
---|
239 | points_list.append(calc_bearing(point_quantities[2], |
---|
240 | point_quantities[3])) |
---|
241 | if quantity == 'xcentroid': |
---|
242 | points_list.append(callable_sww.centroids[point_i][0]) |
---|
243 | |
---|
244 | if quantity == 'ycentroid': |
---|
245 | points_list.append(callable_sww.centroids[point_i][1]) |
---|
246 | |
---|
247 | points_writer[point_i].writerow(points_list) |
---|
248 | |
---|
249 | ## |
---|
250 | # @brief Read a .sww file and plot the time series. |
---|
251 | # @param swwfiles Dictionary of .sww files. |
---|
252 | # @param gauge_filename Name of gauge file. |
---|
253 | # @param production_dirs ?? |
---|
254 | # @param report If True, write figures to report directory. |
---|
255 | # @param reportname Name of generated report (if any). |
---|
256 | # @param plot_quantity List containing quantities to plot. |
---|
257 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
258 | # @param surface If True, then generate solution surface with 3d plot. |
---|
259 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
260 | # @param time_max End of user defined time range for plotting purposes. |
---|
261 | # @param time_thinning ?? |
---|
262 | # @param time_unit ?? |
---|
263 | # @param title_on If True, export standard graphics with title. |
---|
264 | # @param use_cache If True, use caching. |
---|
265 | # @param verbose If True, this function is verbose. |
---|
266 | def gauge_sww2timeseries(swwfiles, |
---|
267 | gauge_filename, |
---|
268 | production_dirs, |
---|
269 | report=None, |
---|
270 | reportname=None, |
---|
271 | plot_quantity=None, |
---|
272 | generate_fig=False, |
---|
273 | surface=None, |
---|
274 | time_min=None, |
---|
275 | time_max=None, |
---|
276 | time_thinning=1, |
---|
277 | time_unit=None, |
---|
278 | title_on=None, |
---|
279 | use_cache=False, |
---|
280 | verbose=False, |
---|
281 | output_centroids=False): |
---|
282 | """ Read sww file and plot the time series for the |
---|
283 | prescribed quantities at defined gauge locations and |
---|
284 | prescribed time range. |
---|
285 | |
---|
286 | Input variables: |
---|
287 | |
---|
288 | swwfiles - dictionary of sww files with label_ids (used in |
---|
289 | generating latex output. It will be part of |
---|
290 | the directory name of file_loc (typically the timestamp). |
---|
291 | Helps to differentiate latex files for different |
---|
292 | simulations for a particular scenario. |
---|
293 | - assume that all conserved quantities have been stored |
---|
294 | - assume each sww file has been simulated with same timestep |
---|
295 | |
---|
296 | gauge_filename - name of file containing gauge data |
---|
297 | - easting, northing, name , elevation? |
---|
298 | - OR (this is not yet done) |
---|
299 | - structure which can be converted to a numeric array, |
---|
300 | such as a geospatial data object |
---|
301 | |
---|
302 | production_dirs - A list of list, example {20061101_121212: '1 in 10000', |
---|
303 | 'boundaries': 'urs boundary'} |
---|
304 | this will use the second part as the label and the |
---|
305 | first part as the ? |
---|
306 | #FIXME: Is it a list or a dictionary |
---|
307 | # This is probably obsolete by now |
---|
308 | |
---|
309 | report - if True, then write figures to report_figures directory in |
---|
310 | relevant production directory |
---|
311 | - if False, figures are already stored with sww file |
---|
312 | - default to False |
---|
313 | |
---|
314 | reportname - name for report if wishing to generate report |
---|
315 | |
---|
316 | plot_quantity - list containing quantities to plot, they must |
---|
317 | be the name of an existing quantity or one of |
---|
318 | the following possibilities |
---|
319 | - possibilities: |
---|
320 | - stage; 'stage' |
---|
321 | - depth; 'depth' |
---|
322 | - speed; calculated as absolute momentum |
---|
323 | (pointwise) divided by depth; 'speed' |
---|
324 | - bearing; calculated as the angle of the momentum |
---|
325 | vector (xmomentum, ymomentum) from the North; 'bearing' |
---|
326 | - absolute momentum; calculated as |
---|
327 | sqrt(xmomentum^2 + ymomentum^2); 'momentum' |
---|
328 | - x momentum; 'xmomentum' |
---|
329 | - y momentum; 'ymomentum' |
---|
330 | - default will be ['stage', 'speed', 'bearing'] |
---|
331 | |
---|
332 | generate_fig - if True, generate figures as well as csv file |
---|
333 | - if False, csv files created only |
---|
334 | |
---|
335 | surface - if True, then generate solution surface with 3d plot |
---|
336 | and save to current working directory |
---|
337 | - default = False |
---|
338 | |
---|
339 | time_min - beginning of user defined time range for plotting purposes |
---|
340 | - default will be first available time found in swwfile |
---|
341 | |
---|
342 | time_max - end of user defined time range for plotting purposes |
---|
343 | - default will be last available time found in swwfile |
---|
344 | |
---|
345 | title_on - if True, export standard graphics with title |
---|
346 | - if False, export standard graphics without title |
---|
347 | |
---|
348 | |
---|
349 | Output: |
---|
350 | |
---|
351 | - time series data stored in .csv for later use if required. |
---|
352 | Name = gauges_timeseries followed by gauge name |
---|
353 | - latex file will be generated in same directory as where script is |
---|
354 | run (usually production scenario directory. |
---|
355 | Name = latexoutputlabel_id.tex |
---|
356 | |
---|
357 | Other important information: |
---|
358 | |
---|
359 | It is assumed that the used has stored all the conserved quantities |
---|
360 | and elevation during the scenario run, i.e. |
---|
361 | ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
362 | If this has not occurred then sww2timeseries will not work. |
---|
363 | |
---|
364 | |
---|
365 | Usage example |
---|
366 | texname = sww2timeseries({project.boundary_name + '.sww': ''}, |
---|
367 | project.polygons_dir + sep + 'boundary_extent.csv', |
---|
368 | project.anuga_dir, |
---|
369 | report = False, |
---|
370 | plot_quantity = ['stage', 'speed', 'bearing'], |
---|
371 | time_min = None, |
---|
372 | time_max = None, |
---|
373 | title_on = True, |
---|
374 | verbose = True) |
---|
375 | |
---|
376 | """ |
---|
377 | |
---|
378 | msg = 'NOTE: A new function is available to create csv files from sww ' |
---|
379 | msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util' |
---|
380 | msg += ' PLUS another new function to create graphs from csv files called ' |
---|
381 | msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util' |
---|
382 | log.critical(msg) |
---|
383 | |
---|
384 | k = _sww2timeseries(swwfiles, |
---|
385 | gauge_filename, |
---|
386 | production_dirs, |
---|
387 | report, |
---|
388 | reportname, |
---|
389 | plot_quantity, |
---|
390 | generate_fig, |
---|
391 | surface, |
---|
392 | time_min, |
---|
393 | time_max, |
---|
394 | time_thinning, |
---|
395 | time_unit, |
---|
396 | title_on, |
---|
397 | use_cache, |
---|
398 | verbose, |
---|
399 | output_centroids = output_centroids) |
---|
400 | return k |
---|
401 | |
---|
402 | |
---|
403 | ## |
---|
404 | # @brief Read a .sww file and plot the time series. |
---|
405 | # @param swwfiles Dictionary of .sww files. |
---|
406 | # @param gauge_filename Name of gauge file. |
---|
407 | # @param production_dirs ?? |
---|
408 | # @param report If True, write figures to report directory. |
---|
409 | # @param reportname Name of generated report (if any). |
---|
410 | # @param plot_quantity List containing quantities to plot. |
---|
411 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
412 | # @param surface If True, then generate solution surface with 3d plot. |
---|
413 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
414 | # @param time_max End of user defined time range for plotting purposes. |
---|
415 | # @param time_thinning ?? |
---|
416 | # @param time_unit ?? |
---|
417 | # @param title_on If True, export standard graphics with title. |
---|
418 | # @param use_cache If True, use caching. |
---|
419 | # @param verbose If True, this function is verbose. |
---|
420 | def _sww2timeseries(swwfiles, |
---|
421 | gauge_filename, |
---|
422 | production_dirs, |
---|
423 | report = None, |
---|
424 | reportname = None, |
---|
425 | plot_quantity = None, |
---|
426 | generate_fig = False, |
---|
427 | surface = None, |
---|
428 | time_min = None, |
---|
429 | time_max = None, |
---|
430 | time_thinning = 1, |
---|
431 | time_unit = None, |
---|
432 | title_on = None, |
---|
433 | use_cache = False, |
---|
434 | verbose = False, |
---|
435 | output_centroids = False): |
---|
436 | |
---|
437 | # FIXME(Ole): Shouldn't print statements here be governed by verbose? |
---|
438 | assert type(gauge_filename) == type(''), 'Gauge filename must be a string' |
---|
439 | |
---|
440 | try: |
---|
441 | fid = open(gauge_filename) |
---|
442 | except Exception, e: |
---|
443 | msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e) |
---|
444 | raise msg |
---|
445 | |
---|
446 | if report is None: |
---|
447 | report = False |
---|
448 | |
---|
449 | if plot_quantity is None: |
---|
450 | plot_quantity = ['depth', 'speed'] |
---|
451 | else: |
---|
452 | assert type(plot_quantity) == list, 'plot_quantity must be a list' |
---|
453 | check_list(plot_quantity) |
---|
454 | |
---|
455 | if surface is None: |
---|
456 | surface = False |
---|
457 | |
---|
458 | if time_unit is None: |
---|
459 | time_unit = 'hours' |
---|
460 | |
---|
461 | if title_on is None: |
---|
462 | title_on = True |
---|
463 | |
---|
464 | if verbose: log.critical('Gauges obtained from: %s' % gauge_filename) |
---|
465 | |
---|
466 | gauges, locations, elev = get_gauges_from_file(gauge_filename) |
---|
467 | |
---|
468 | sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
469 | |
---|
470 | file_loc = [] |
---|
471 | f_list = [] |
---|
472 | label_id = [] |
---|
473 | leg_label = [] |
---|
474 | themaxT = 0.0 |
---|
475 | theminT = 0.0 |
---|
476 | |
---|
477 | for swwfile in swwfiles.keys(): |
---|
478 | try: |
---|
479 | fid = open(swwfile) |
---|
480 | except Exception, e: |
---|
481 | msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e) |
---|
482 | raise msg |
---|
483 | |
---|
484 | if verbose: |
---|
485 | log.critical('swwfile = %s' % swwfile) |
---|
486 | |
---|
487 | # Extract parent dir name and use as label |
---|
488 | path, _ = os.path.split(swwfile) |
---|
489 | _, label = os.path.split(path) |
---|
490 | |
---|
491 | leg_label.append(label) |
---|
492 | |
---|
493 | f = file_function(swwfile, |
---|
494 | quantities = sww_quantity, |
---|
495 | interpolation_points = gauges, |
---|
496 | time_thinning = time_thinning, |
---|
497 | verbose = verbose, |
---|
498 | use_cache = use_cache, |
---|
499 | output_centroids = output_centroids) |
---|
500 | |
---|
501 | # determine which gauges are contained in sww file |
---|
502 | count = 0 |
---|
503 | gauge_index = [] |
---|
504 | for k, g in enumerate(gauges): |
---|
505 | if f(0.0, point_id = k)[2] > 1.0e6: |
---|
506 | count += 1 |
---|
507 | if count == 1: log.critical('Gauges not contained here:') |
---|
508 | log.critical(locations[k]) |
---|
509 | else: |
---|
510 | gauge_index.append(k) |
---|
511 | |
---|
512 | if len(gauge_index) > 0: |
---|
513 | log.critical('Gauges contained here:') |
---|
514 | else: |
---|
515 | log.critical('No gauges contained here.') |
---|
516 | for i in range(len(gauge_index)): |
---|
517 | log.critical(locations[gauge_index[i]]) |
---|
518 | |
---|
519 | index = swwfile.rfind(sep) |
---|
520 | file_loc.append(swwfile[:index+1]) |
---|
521 | label_id.append(swwfiles[swwfile]) |
---|
522 | |
---|
523 | f_list.append(f) |
---|
524 | maxT = max(f.get_time()) |
---|
525 | minT = min(f.get_time()) |
---|
526 | if maxT > themaxT: themaxT = maxT |
---|
527 | if minT > theminT: theminT = minT |
---|
528 | |
---|
529 | if time_min is None: |
---|
530 | time_min = theminT # min(T) |
---|
531 | else: |
---|
532 | if time_min < theminT: # min(T): |
---|
533 | msg = 'Minimum time entered not correct - please try again' |
---|
534 | raise Exception, msg |
---|
535 | |
---|
536 | if time_max is None: |
---|
537 | time_max = themaxT # max(T) |
---|
538 | else: |
---|
539 | if time_max > themaxT: # max(T): |
---|
540 | msg = 'Maximum time entered not correct - please try again' |
---|
541 | raise Exception, msg |
---|
542 | |
---|
543 | if verbose and len(gauge_index) > 0: |
---|
544 | log.critical('Inputs OK - going to generate figures') |
---|
545 | |
---|
546 | if len(gauge_index) <> 0: |
---|
547 | texfile, elev_output = \ |
---|
548 | generate_figures(plot_quantity, file_loc, report, reportname, |
---|
549 | surface, leg_label, f_list, gauges, locations, |
---|
550 | elev, gauge_index, production_dirs, time_min, |
---|
551 | time_max, time_unit, title_on, label_id, |
---|
552 | generate_fig, verbose) |
---|
553 | else: |
---|
554 | texfile = '' |
---|
555 | elev_output = [] |
---|
556 | |
---|
557 | return texfile, elev_output |
---|
558 | |
---|
559 | |
---|
560 | ## |
---|
561 | # @brief Read gauge info from a file. |
---|
562 | # @param filename The name of the file to read. |
---|
563 | # @return A (gauges, gaugelocation, elev) tuple. |
---|
564 | def get_from_file(filename): |
---|
565 | """ Read in gauge information from file |
---|
566 | """ |
---|
567 | |
---|
568 | from os import sep, getcwd, access, F_OK, mkdir |
---|
569 | |
---|
570 | # Get data from the gauge file |
---|
571 | fid = open(filename) |
---|
572 | lines = fid.readlines() |
---|
573 | fid.close() |
---|
574 | |
---|
575 | gauges = [] |
---|
576 | gaugelocation = [] |
---|
577 | elev = [] |
---|
578 | |
---|
579 | # Check header information |
---|
580 | line1 = lines[0] |
---|
581 | line11 = line1.split(',') |
---|
582 | |
---|
583 | if isinstance(line11[0], str) is True: |
---|
584 | # We have found text in the first line |
---|
585 | east_index = None |
---|
586 | north_index = None |
---|
587 | name_index = None |
---|
588 | elev_index = None |
---|
589 | |
---|
590 | for i in range(len(line11)): |
---|
591 | if line11[i].strip().lower() == 'easting': east_index = i |
---|
592 | if line11[i].strip().lower() == 'northing': north_index = i |
---|
593 | if line11[i].strip().lower() == 'name': name_index = i |
---|
594 | if line11[i].strip().lower() == 'elevation': elev_index = i |
---|
595 | |
---|
596 | if east_index < len(line11) and north_index < len(line11): |
---|
597 | pass |
---|
598 | else: |
---|
599 | msg = 'WARNING: %s does not contain correct header information' \ |
---|
600 | % filename |
---|
601 | msg += 'The header must be: easting, northing, name, elevation' |
---|
602 | raise Exception, msg |
---|
603 | |
---|
604 | if elev_index is None: |
---|
605 | raise Exception |
---|
606 | |
---|
607 | if name_index is None: |
---|
608 | raise Exception |
---|
609 | |
---|
610 | lines = lines[1:] # Remove header from data |
---|
611 | else: |
---|
612 | # No header, assume that this is a simple easting, northing file |
---|
613 | |
---|
614 | msg = 'There was no header in file %s and the number of columns is %d' \ |
---|
615 | % (filename, len(line11)) |
---|
616 | msg += '- was assuming two columns corresponding to Easting and Northing' |
---|
617 | assert len(line11) == 2, msg |
---|
618 | |
---|
619 | east_index = 0 |
---|
620 | north_index = 1 |
---|
621 | |
---|
622 | N = len(lines) |
---|
623 | elev = [-9999]*N |
---|
624 | gaugelocation = range(N) |
---|
625 | |
---|
626 | # Read in gauge data |
---|
627 | for line in lines: |
---|
628 | fields = line.split(',') |
---|
629 | |
---|
630 | gauges.append([float(fields[east_index]), float(fields[north_index])]) |
---|
631 | |
---|
632 | if len(fields) > 2: |
---|
633 | elev.append(float(fields[elev_index])) |
---|
634 | loc = fields[name_index] |
---|
635 | gaugelocation.append(loc.strip('\n')) |
---|
636 | |
---|
637 | return gauges, gaugelocation, elev |
---|