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

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

Cleaned up unit tests to use API.

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