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

Last change on this file since 8149 was 8149, checked in by wilsonr, 13 years ago

Removed '@brief' comments.

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