source: trunk/anuga_core/source/anuga/file_conversion/urs2sww.py @ 7770

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

Broke up data_manager into more modules - some failing tests.

File size: 8.8 KB
Line 
1import numpy as num
2from Scientific.IO.NetCDF import NetCDFFile
3
4from anuga.file.urs import Read_urs
5
6from anuga.coordinate_transforms.redfearn import redfearn, \
7     convert_from_latlon_to_utm
8
9from anuga.geospatial_data.geospatial_data import ensure_absolute, \
10                                                    Geospatial_data
11
12from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
13                            NORTH_VELOCITY_LABEL
14                           
15from anuga.utilities.numerical_tools import ensure_numeric                           
16
17from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
18                            netcdf_float
19
20from sww_file import Read_sww, Write_sww
21
22
23################################################################################
24# CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
25################################################################################
26
27##
28# @brief Convert URS file(s) (wave prop) to an SWW file.
29# @param basename_in Stem of the input filenames.
30# @param basename_out Path to the output SWW file.
31# @param verbose True if this function is to be verbose.
32# @param mint
33# @param maxt
34# @param mean_stage
35# @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing).
36# @param hole_points_UTM
37# @param zscale
38# @note Also convert latitude and longitude to UTM. All coordinates are
39#       assumed to be given in the GDA94 datum.
40# @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux'
41#       added for relative height, x-velocity and y-velocity respectively.
42def urs2sww(basename_in='o', basename_out=None, verbose=False,
43                      mint=None, maxt=None,
44                      mean_stage=0,
45                      origin=None,
46                      hole_points_UTM=None,
47                      zscale=1):
48    """
49    Convert URS C binary format for wave propagation to
50    sww format native to abstract_2d_finite_volumes.
51
52    Specify only basename_in and read files of the form
53    basefilename-z-mux, basefilename-e-mux and
54    basefilename-n-mux containing relative height,
55    x-velocity and y-velocity, respectively.
56
57    Also convert latitude and longitude to UTM. All coordinates are
58    assumed to be given in the GDA94 datum. The latitude and longitude
59    information is assumed ungridded grid.
60
61    min's and max's: If omitted - full extend is used.
62    To include a value min ans max may equal it.
63    Lat and lon are assumed to be in decimal degrees.
64
65    origin is a 3-tuple with geo referenced
66    UTM coordinates (zone, easting, northing)
67    It will be the origin of the sww file. This shouldn't be used,
68    since all of anuga should be able to handle an arbitary origin.
69    The mux point info is NOT relative to this origin.
70
71    URS C binary format has data organised as TIME, LONGITUDE, LATITUDE
72    which means that latitude is the fastest
73    varying dimension (row major order, so to speak)
74
75    In URS C binary the latitudes and longitudes are in assending order.
76
77    Note, interpolations of the resulting sww file will be different
78    from results of urs2sww.  This is due to the interpolation
79    function used, and the different grid structure between urs2sww
80    and this function.
81
82    Interpolating data that has an underlying gridded source can
83    easily end up with different values, depending on the underlying
84    mesh.
85
86    consider these 4 points
87    50  -50
88
89    0     0
90
91    The grid can be
92     -
93    |\|   A
94     -
95     or;
96      -
97     |/|  B
98      -
99
100    If a point is just below the center of the midpoint, it will have a
101    +ve value in grid A and a -ve value in grid B.
102    """
103
104    from anuga.mesh_engine.mesh_engine import NoTrianglesError
105    from anuga.pmesh.mesh import Mesh
106
107    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
108                basename_in + EAST_VELOCITY_LABEL,
109                basename_in + NORTH_VELOCITY_LABEL]
110    quantities = ['HA','UA','VA']
111
112    # instantiate urs_points of the three mux files.
113    mux = {}
114    for quantity, file in map(None, quantities, files_in):
115        mux[quantity] = Read_urs(file)
116
117    # Could check that the depth is the same. (hashing)
118
119    # handle to a mux file to do depth stuff
120    a_mux = mux[quantities[0]]
121
122    # Convert to utm
123    lat = a_mux.lonlatdep[:,1]
124    long = a_mux.lonlatdep[:,0]
125    points_utm, zone = convert_from_latlon_to_utm(latitudes=lat,
126                                                  longitudes=long)
127
128    elevation = a_mux.lonlatdep[:,2] * -1
129
130    # grid (create a mesh from the selected points)
131    # This mesh has a problem.  Triangles are streched over ungridded areas.
132    # If these areas could be described as holes in pmesh, that would be great.
133
134    # I can't just get the user to selection a point in the middle.
135    # A boundary is needed around these points.
136    # But if the zone of points is obvious enough auto-segment should do
137    # a good boundary.
138    mesh = Mesh()
139    mesh.add_vertices(points_utm)
140    mesh.auto_segment(smooth_indents=True, expand_pinch=True)
141
142    # To try and avoid alpha shape 'hugging' too much
143    mesh.auto_segment(mesh.shape.get_alpha() * 1.1)
144    if hole_points_UTM is not None:
145        point = ensure_absolute(hole_points_UTM)
146        mesh.add_hole(point[0], point[1])
147
148    try:
149        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
150    except NoTrianglesError:
151        # This is a bit of a hack, going in and changing the data structure.
152        mesh.holes = []
153        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
154
155    mesh_dic = mesh.Mesh2MeshList()
156
157    #mesh.export_mesh_file(basename_in + '_168.tsh')
158    #import sys; sys.exit()
159    # These are the times of the mux file
160    mux_times = []
161    for i in range(a_mux.time_step_count):
162        mux_times.append(a_mux.time_step * i)
163    (mux_times_start_i, mux_times_fin_i) = read_time_from_mux(mux_times, mint, maxt)
164    times = mux_times[mux_times_start_i:mux_times_fin_i]
165
166    if mux_times_start_i == mux_times_fin_i:
167        # Close the mux files
168        for quantity, file in map(None, quantities, files_in):
169            mux[quantity].close()
170        msg = "Due to mint and maxt there's no time info in the boundary SWW."
171        raise Exception, msg
172
173    # If this raise is removed there is currently no downstream errors
174
175    points_utm=ensure_numeric(points_utm)
176    assert num.alltrue(ensure_numeric(mesh_dic['generatedpointlist'])
177                       == ensure_numeric(points_utm))
178
179    volumes = mesh_dic['generatedtrianglelist']
180
181    # Write sww intro and grid stuff.
182    if basename_out is None:
183        swwname = basename_in + '.sww'
184    else:
185        swwname = basename_out + '.sww'
186
187    if verbose: log.critical('Output to %s' % swwname)
188
189    outfile = NetCDFFile(swwname, netcdf_mode_w)
190
191    # For a different way of doing this, check out tsh2sww
192    # work out sww_times and the index range this covers
193    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
194    sww.store_header(outfile, times, len(volumes), len(points_utm),
195                     verbose=verbose, sww_precision=netcdf_float)
196    outfile.mean_stage = mean_stage
197    outfile.zscale = zscale
198
199    sww.store_triangulation(outfile, points_utm, volumes,
200                            zone, 
201                            new_origin=origin,
202                            verbose=verbose)
203    sww.store_static_quantities(outfile, elevation=elevation)
204
205    if verbose: log.critical('Converting quantities')
206
207    # Read in a time slice from each mux file and write it to the SWW file
208    j = 0
209    for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
210        if j >= mux_times_start_i and j < mux_times_fin_i:
211            stage = zscale*ha + mean_stage
212            h = stage - elevation
213            xmomentum = ua*h
214            ymomentum = -1 * va * h # -1 since in mux files south is positive.
215            sww.store_quantities(outfile,
216                                 slice_index=j-mux_times_start_i,
217                                 verbose=verbose,
218                                 stage=stage,
219                                 xmomentum=xmomentum,
220                                 ymomentum=ymomentum,
221                                 sww_precision=num.float)
222        j += 1
223
224    if verbose: sww.verbose_quantities(outfile)
225
226    outfile.close()
227
228    # Do some conversions while writing the sww file
229
230
231def read_time_from_mux(mux_times, mint, maxt):
232    """
233        Read a list of mux times.
234        Return start and finish times which lie within the passed time period.
235    """
236
237    if mint == None:
238        mux_times_start_i = 0
239    else:
240        mux_times_start_i = num.searchsorted(mux_times, mint)
241
242    if maxt == None:
243        mux_times_fin_i = len(mux_times)
244    else:
245        maxt += 0.5 # so if you specify a time where there is
246                    # data that time will be included
247        mux_times_fin_i = num.searchsorted(mux_times, maxt)
248
249    return mux_times_start_i, mux_times_fin_i
250
Note: See TracBrowser for help on using the repository browser.