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

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

urs2sww has an extra urs_ungridded2sww function.

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