source: trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py @ 7858

Last change on this file since 7858 was 7858, checked in by hudson, 14 years ago

Refactorings to increase code quality, fixed missing log import from sww_interrogate.

File size: 19.6 KB
Line 
1'''
2    Operations to extract information from an SWW file.
3'''
4
5import os
6import anuga.utilities.log as log
7import numpy as num
8
9from anuga.utilities.file_utils import get_all_swwfiles
10from anuga.coordinate_transforms.geo_reference import Geo_reference
11from anuga.abstract_2d_finite_volumes.util import file_function
12from anuga.geometry.polygon import is_inside_polygon
13from anuga.file.sww import get_mesh_and_quantities_from_file
14from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
15
16##
17# @brief Get values for quantities interpolated to polyline midpoints from SWW.
18# @param filename Path to file to read.
19# @param quantity_names Quantity names to get.
20# @param polyline Representation of desired cross-section.
21# @param verbose True if this function is to be verbose.
22# @return (segments, i_func) where segments is a list of Triangle_intersection
23#         instances and i_func is an instance of Interpolation_function.
24# @note For 'polyline' assume absolute UTM coordinates.
25def get_interpolated_quantities_at_polyline_midpoints(filename,
26                                                      quantity_names=None,
27                                                      polyline=None,
28                                                      verbose=False):
29    """Get values for quantities interpolated to polyline midpoints from SWW
30
31    Input:
32        filename - Name of sww file
33        quantity_names - Names of quantities to load
34        polyline: Representation of desired cross section - it may contain
35                  multiple sections allowing for complex shapes. Assume
36                  absolute UTM coordinates.
37                  Format [[x0, y0], [x1, y1], ...]
38
39    Output:
40        segments: list of instances of class Triangle_intersection
41        interpolation_function: Instance of class Interpolation_function
42
43
44      This function is used by get_flow_through_cross_section and
45      get_energy_through_cross_section
46    """
47
48    from anuga.fit_interpolate.interpolate import Interpolation_function
49
50    # Get mesh and quantities from sww file
51    X = get_mesh_and_quantities_from_file(filename,
52                                          quantities=quantity_names,
53                                          verbose=verbose)
54    mesh, quantities, time = X
55
56    # Find all intersections and associated triangles.
57    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
58
59    # Get midpoints
60    interpolation_points = segment_midpoints(segments)
61
62    # Interpolate
63    if verbose:
64        log.critical('Interpolating - total number of interpolation points = %d'
65                     % len(interpolation_points))
66
67    I = Interpolation_function(time,
68                               quantities,
69                               quantity_names=quantity_names,
70                               vertex_coordinates=mesh.nodes,
71                               triangles=mesh.triangles,
72                               interpolation_points=interpolation_points,
73                               verbose=verbose)
74
75    return segments, I
76
77
78##
79# @brief Obtain flow (m^3/s) perpendicular to specified cross section.
80# @param filename Path to file to read.
81# @param polyline Representation of desired cross-section.
82# @param verbose Trie if this function is to be verbose.
83# @return (time, Q) where time and Q are lists of time and flow respectively.
84def get_flow_through_cross_section(filename, polyline, verbose=False):
85    """Obtain flow (m^3/s) perpendicular to specified cross section.
86
87    Inputs:
88        filename: Name of sww file
89        polyline: Representation of desired cross section - it may contain
90                  multiple sections allowing for complex shapes. Assume
91                  absolute UTM coordinates.
92                  Format [[x0, y0], [x1, y1], ...]
93
94    Output:
95        time: All stored times in sww file
96        Q: Hydrograph of total flow across given segments for all stored times.
97
98    The normal flow is computed for each triangle intersected by the polyline
99    and added up.  Multiple segments at different angles are specified the
100    normal flows may partially cancel each other.
101
102    The typical usage of this function would be to get flow through a channel,
103    and the polyline would then be a cross section perpendicular to the flow.
104    """
105
106    quantity_names =['elevation',
107                     'stage',
108                     'xmomentum',
109                     'ymomentum']
110
111    # Get values for quantities at each midpoint of poly line from sww file
112    X = get_interpolated_quantities_at_polyline_midpoints(filename,
113                                                          quantity_names=\
114                                                              quantity_names,
115                                                          polyline=polyline,
116                                                          verbose=verbose)
117    segments, interpolation_function = X
118
119    # Get vectors for time and interpolation_points
120    time = interpolation_function.time
121    interpolation_points = interpolation_function.interpolation_points
122
123    if verbose: log.critical('Computing hydrograph')
124
125    # Compute hydrograph
126    Q = []
127    for t in time:
128        total_flow = 0
129        for i in range(len(interpolation_points)):
130            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
131            normal = segments[i].normal
132
133            # Inner product of momentum vector with segment normal [m^2/s]
134            normal_momentum = uh*normal[0] + vh*normal[1]
135
136            # Flow across this segment [m^3/s]
137            segment_flow = normal_momentum * segments[i].length
138
139            # Accumulate
140            total_flow += segment_flow
141
142        # Store flow at this timestep
143        Q.append(total_flow)
144
145
146    return time, Q
147
148
149##
150# @brief Get average energy across a cross-section.
151# @param filename Path to file of interest.
152# @param polyline Representation of desired cross-section.
153# @param kind Select energy to compute: 'specific' or 'total'.
154# @param verbose True if this function is to be verbose.
155# @return (time, E) where time and E are lists of timestep and energy.
156def get_energy_through_cross_section(filename,
157                                     polyline,
158                                     kind='total',
159                                     verbose=False):
160    """Obtain average energy head [m] across specified cross section.
161
162    Inputs:
163        polyline: Representation of desired cross section - it may contain
164                  multiple sections allowing for complex shapes. Assume
165                  absolute UTM coordinates.
166                  Format [[x0, y0], [x1, y1], ...]
167        kind:     Select which energy to compute.
168                  Options are 'specific' and 'total' (default)
169
170    Output:
171        E: Average energy [m] across given segments for all stored times.
172
173    The average velocity is computed for each triangle intersected by
174    the polyline and averaged weighted by segment lengths.
175
176    The typical usage of this function would be to get average energy of
177    flow in a channel, and the polyline would then be a cross section
178    perpendicular to the flow.
179
180    #FIXME (Ole) - need name for this energy reflecting that its dimension
181    is [m].
182    """
183
184    from anuga.config import g, epsilon, velocity_protection as h0
185
186    quantity_names =['elevation',
187                     'stage',
188                     'xmomentum',
189                     'ymomentum']
190
191    # Get values for quantities at each midpoint of poly line from sww file
192    X = get_interpolated_quantities_at_polyline_midpoints(filename,
193                                                          quantity_names=\
194                                                              quantity_names,
195                                                          polyline=polyline,
196                                                          verbose=verbose)
197    segments, interpolation_function = X
198
199    # Get vectors for time and interpolation_points
200    time = interpolation_function.time
201    interpolation_points = interpolation_function.interpolation_points
202
203    if verbose: log.critical('Computing %s energy' % kind)
204
205    # Compute total length of polyline for use with weighted averages
206    total_line_length = 0.0
207    for segment in segments:
208        total_line_length += segment.length
209
210    # Compute energy
211    E = []
212    for t in time:
213        average_energy = 0.0
214        for i, p in enumerate(interpolation_points):
215            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
216
217            # Depth
218            h = depth = stage-elevation
219
220            # Average velocity across this segment
221            if h > epsilon:
222                # Use protection against degenerate velocities
223                u = uh / (h + h0/h)
224                v = vh / (h + h0/h)
225            else:
226                u = v = 0.0
227
228            speed_squared = u*u + v*v
229            kinetic_energy = 0.5 * speed_squared / g
230
231            if kind == 'specific':
232                segment_energy = depth + kinetic_energy
233            elif kind == 'total':
234                segment_energy = stage + kinetic_energy
235            else:
236                msg = 'Energy kind must be either "specific" or "total". '
237                msg += 'I got %s' % kind
238
239            # Add to weighted average
240            weigth = segments[i].length / total_line_length
241            average_energy += segment_energy * weigth
242
243        # Store energy at this timestep
244        E.append(average_energy)
245
246    return time, E
247
248
249##
250# @brief Return highest elevation where depth > 0.
251# @param filename Path to SWW file of interest.
252# @param polygon If specified resrict to points inside this polygon.
253# @param time_interval If specified resrict to within the time specified.
254# @param verbose True if this function is  to be verbose.
255def get_maximum_inundation_elevation(filename,
256                                     polygon=None,
257                                     time_interval=None,
258                                     verbose=False):
259    """Return highest elevation where depth > 0
260
261    Usage:
262    max_runup = get_maximum_inundation_elevation(filename,
263                                                 polygon=None,
264                                                 time_interval=None,
265                                                 verbose=False)
266
267    filename is a NetCDF sww file containing ANUGA model output.
268    Optional arguments polygon and time_interval restricts the maximum
269    runup calculation
270    to a points that lie within the specified polygon and time interval.
271
272    If no inundation is found within polygon and time_interval the return value
273    is None signifying "No Runup" or "Everything is dry".
274
275    See general function get_maximum_inundation_data for details.
276    """
277
278    runup, _ = get_maximum_inundation_data(filename,
279                                           polygon=polygon,
280                                           time_interval=time_interval,
281                                           verbose=verbose)
282    return runup
283
284
285##
286# @brief Return location of highest elevation where h > 0
287# @param filename Path to SWW file to read.
288# @param polygon If specified resrict to points inside this polygon.
289# @param time_interval If specified resrict to within the time specified.
290# @param verbose True if this function is  to be verbose.
291def get_maximum_inundation_location(filename,
292                                    polygon=None,
293                                    time_interval=None,
294                                    verbose=False):
295    """Return location of highest elevation where h > 0
296
297    Usage:
298    max_runup_location = get_maximum_inundation_location(filename,
299                                                         polygon=None,
300                                                         time_interval=None,
301                                                         verbose=False)
302
303    filename is a NetCDF sww file containing ANUGA model output.
304    Optional arguments polygon and time_interval restricts the maximum
305    runup calculation
306    to a points that lie within the specified polygon and time interval.
307
308    If no inundation is found within polygon and time_interval the return value
309    is None signifying "No Runup" or "Everything is dry".
310
311    See general function get_maximum_inundation_data for details.
312    """
313
314    _, max_loc = get_maximum_inundation_data(filename,
315                                             polygon=polygon,
316                                             time_interval=time_interval,
317                                             verbose=verbose)
318    return max_loc
319
320
321##
322# @brief Compute maximum run up height from SWW file.
323# @param filename Path to SWW file to read.
324# @param polygon If specified resrict to points inside this polygon.
325# @param time_interval If specified resrict to within the time specified.
326# @param use_centroid_values
327# @param verbose True if this function is to be verbose.
328# @return (maximal_runup, maximal_runup_location)
329def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
330                                use_centroid_values=False,
331                                verbose=False):
332    """Compute maximum run up height from sww file.
333
334    Usage:
335    runup, location = get_maximum_inundation_data(filename,
336                                                  polygon=None,
337                                                  time_interval=None,
338                                                  verbose=False)
339
340    Algorithm is as in get_maximum_inundation_elevation from
341    shallow_water_domain except that this function works with the sww file and
342    computes the maximal runup height over multiple timesteps.
343
344    Optional arguments polygon and time_interval restricts the maximum runup
345    calculation to a points that lie within the specified polygon and time
346    interval.
347
348    Polygon is assumed to be in (absolute) UTM coordinates in the same zone
349    as domain.
350
351    If no inundation is found within polygon and time_interval the return value
352    is None signifying "No Runup" or "Everything is dry".
353    """
354
355    # We are using nodal values here as that is what is stored in sww files.
356
357    # Water depth below which it is considered to be 0 in the model
358    # FIXME (Ole): Allow this to be specified as a keyword argument as well
359
360    from anuga.geometry.polygon import inside_polygon
361    from anuga.config import minimum_allowed_height
362    from Scientific.IO.NetCDF import NetCDFFile
363
364    dir, base = os.path.split(filename)
365
366    iterate_over = get_all_swwfiles(dir, base)
367
368    # Read sww file
369    if verbose: log.critical('Reading from %s' % filename)
370    # FIXME: Use general swwstats (when done)
371
372    maximal_runup = None
373    maximal_runup_location = None
374
375    for _, swwfile in enumerate (iterate_over):
376        # Read sww file
377        filename = os.path.join(dir, swwfile+'.sww')
378
379        if verbose: log.critical('Reading from %s' % filename)
380        # FIXME: Use general swwstats (when done)
381
382        fid = NetCDFFile(filename)
383
384        # Get geo_reference
385        # sww files don't have to have a geo_ref
386        try:
387            geo_reference = Geo_reference(NetCDFObject=fid)
388        except AttributeError:
389            geo_reference = Geo_reference() # Default georef object
390
391        xllcorner = geo_reference.get_xllcorner()
392        yllcorner = geo_reference.get_yllcorner()
393
394        # Get extent
395        volumes = fid.variables['volumes'][:]
396        x = fid.variables['x'][:] + xllcorner
397        y = fid.variables['y'][:] + yllcorner
398
399        # Get the relevant quantities (Convert from single precison)
400        elevation = num.array(fid.variables['elevation'][:], num.float)
401        stage = num.array(fid.variables['stage'][:], num.float)
402
403        # Here's where one could convert nodal information to centroid
404        # information but is probably something we need to write in C.
405        # Here's a Python thought which is NOT finished!!!
406        if use_centroid_values is True:
407            x = get_centroid_values(x, volumes)
408            y = get_centroid_values(y, volumes)
409            elevation = get_centroid_values(elevation, volumes)
410
411        # Spatial restriction
412        if polygon is not None:
413            msg = 'polygon must be a sequence of points.'
414            assert len(polygon[0]) == 2, msg
415            # FIXME (Ole): Make a generic polygon input check in polygon.py
416            # and call it here
417            points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis],
418                                                            y[:, num.newaxis]),
419                                                            axis=1))
420            point_indices = inside_polygon(points, polygon)
421
422            # Restrict quantities to polygon
423            elevation = num.take(elevation, point_indices, axis=0)
424            stage = num.take(stage, point_indices, axis=1)
425
426            # Get info for location of maximal runup
427            points_in_polygon = num.take(points, point_indices, axis=0)
428
429            x = points_in_polygon[:,0]
430            y = points_in_polygon[:,1]
431        else:
432            # Take all points
433            point_indices = num.arange(len(x))
434
435        # Temporal restriction
436        time = fid.variables['time'][:]
437        all_timeindices = num.arange(len(time))
438        if time_interval is not None:
439            msg = 'time_interval must be a sequence of length 2.'
440            assert len(time_interval) == 2, msg
441            msg = 'time_interval %s must not be decreasing.' % time_interval
442            assert time_interval[1] >= time_interval[0], msg
443            msg = 'Specified time interval [%.8f:%.8f] ' % tuple(time_interval)
444            msg += 'must does not match model time interval: [%.8f, %.8f]\n' \
445                   % (time[0], time[-1])
446            if time_interval[1] < time[0]: raise ValueError(msg)
447            if time_interval[0] > time[-1]: raise ValueError(msg)
448
449            # Take time indices corresponding to interval (& is bitwise AND)
450            timesteps = num.compress((time_interval[0] <= time) \
451                                     & (time <= time_interval[1]),
452                                     all_timeindices)
453
454            msg = 'time_interval %s did not include any model timesteps.' \
455                  % time_interval
456            assert not num.alltrue(timesteps == 0), msg
457        else:
458            # Take them all
459            timesteps = all_timeindices
460
461        fid.close()
462
463        # Compute maximal runup for each timestep
464        #maximal_runup = None
465        #maximal_runup_location = None
466        #maximal_runups = [None]
467        #maximal_runup_locations = [None]
468
469        for i in timesteps:
470            if use_centroid_values is True:
471                stage_i = get_centroid_values(stage[i,:], volumes)
472            else:
473                stage_i = stage[i,:]
474
475            depth = stage_i - elevation
476
477            # Get wet nodes i.e. nodes with depth>0 within given region
478            # and timesteps
479            wet_nodes = num.compress(depth > minimum_allowed_height,
480                                     num.arange(len(depth)))
481
482            if num.alltrue(wet_nodes == 0):
483                runup = None
484            else:
485                # Find maximum elevation among wet nodes
486                wet_elevation = num.take(elevation, wet_nodes, axis=0)
487                runup_index = num.argmax(wet_elevation)
488                runup = max(wet_elevation)
489                assert wet_elevation[runup_index] == runup       # Must be True
490
491            if runup > maximal_runup:
492                maximal_runup = runup      # works even if maximal_runup is None
493
494                # Record location
495                wet_x = num.take(x, wet_nodes, axis=0)
496                wet_y = num.take(y, wet_nodes, axis=0)
497                maximal_runup_location =    [wet_x[runup_index], \
498                                            wet_y[runup_index]]
499
500    return maximal_runup, maximal_runup_location
501
Note: See TracBrowser for help on using the repository browser.