source: trunk/anuga_core/source/anuga/file_conversion/file_conversion.py @ 7814

Last change on this file since 7814 was 7814, checked in by hudson, 13 years ago

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.

File size: 11.5 KB
Line 
1""" Conversion routines.
2    ANUGA needs to deal with many different file formats, and this
3    module provides routines for easily converting between them.
4   
5    These routines are necessarily high level, sitting above the various
6    ANUGA modules. They take a file as input, and output a file.
7"""
8
9#non ANUGA imports
10from Scientific.IO.NetCDF import NetCDFFile
11import numpy as num
12import os.path
13import os
14
15#ANUGA imports
16from anuga.coordinate_transforms.geo_reference import Geo_reference, \
17     write_NetCDF_georeference, ensure_geo_reference
18from anuga.abstract_2d_finite_volumes.pmesh2domain import \
19     pmesh_to_domain_instance
20from anuga.utilities.numerical_tools import ensure_numeric
21from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
22                            netcdf_float
23
24
25#shallow water imports
26from anuga.file.sww import Read_sww, Write_sww
27from anuga.shallow_water.shallow_water_domain import Domain
28from anuga.shallow_water.shallow_water_domain import Domain
29
30
31def sww2obj(filename, size):
32    """ Convert netcdf based data output to obj
33
34        Convert SWW data to OBJ data.
35        basefilename Stem of filename, needs size and extension added.
36        size The number of lines to write.
37    """
38
39    from Scientific.IO.NetCDF import NetCDFFile
40
41    if filename[-4:] != '.sww':
42        raise IOError('Output file %s should be of type .sww.' % sww_file)
43
44    basefilename = filename[:-4]
45
46    # Get NetCDF
47    nc_fname = create_filename('.', basefilename, 'sww', size)
48    log.critical('Reading from %s' % nc_fname)
49    fid = NetCDFFile(nc_fname, netcdf_mode_r)  #Open existing file for read
50
51    # Get the variables
52    x = fid.variables['x']
53    y = fid.variables['y']
54    z = fid.variables['elevation']
55    time = fid.variables['time']
56    stage = fid.variables['stage']
57
58    M = size  #Number of lines
59    xx = num.zeros((M,3), num.float)
60    yy = num.zeros((M,3), num.float)
61    zz = num.zeros((M,3), num.float)
62
63    for i in range(M):
64        for j in range(3):
65            xx[i,j] = x[i+j*M]
66            yy[i,j] = y[i+j*M]
67            zz[i,j] = z[i+j*M]
68
69    # Write obj for bathymetry
70    FN = create_filename('.', basefilename, 'obj', size)
71    write_obj(FN,xx,yy,zz)
72
73    # Now read all the data with variable information, combine with
74    # x,y info and store as obj
75    for k in range(len(time)):
76        t = time[k]
77        log.critical('Processing timestep %f' % t)
78
79        for i in range(M):
80            for j in range(3):
81                zz[i,j] = stage[k,i+j*M]
82
83        #Write obj for variable data
84        #FN = create_filename(basefilename, 'obj', size, time=t)
85        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
86        write_obj(FN, xx, yy, zz)
87
88##
89# @brief Convert time-series text file to TMS file.
90# @param filename
91# @param quantity_names
92# @param time_as_seconds
93def timefile2netcdf(file_text, quantity_names=None, time_as_seconds=False):
94    """Template for converting typical text files with time series to
95    NetCDF tms file.
96
97    The file format is assumed to be either two fields separated by a comma:
98
99        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
100
101    E.g
102
103      31/08/04 00:00:00, 1.328223 0 0
104      31/08/04 00:15:00, 1.292912 0 0
105
106    or time (seconds), value0 value1 value2 ...
107
108      0.0, 1.328223 0 0
109      0.1, 1.292912 0 0
110
111    will provide a time dependent function f(t) with three attributes
112
113    filename is assumed to be the rootname with extenisons .txt and .sww
114    """
115
116    import time, calendar
117    from anuga.config import time_format
118    from anuga.utilities.numerical_tools import ensure_numeric
119
120    if file_text[-4:] != '.txt':
121        raise IOError('Input file %s should be of type .txt.' % file_text)
122
123    filename = file_text[:-4]
124    fid = open(file_text)
125    line = fid.readline()
126    fid.close()
127
128    fields = line.split(',')
129    msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \
130          % file_text
131    assert len(fields) == 2, msg
132
133    if not time_as_seconds:
134        try:
135            starttime = calendar.timegm(time.strptime(fields[0], time_format))
136        except ValueError:
137            msg = 'First field in file %s must be' % file_text
138            msg += ' date-time with format %s.\n' % time_format
139            msg += 'I got %s instead.' % fields[0]
140            raise DataTimeError, msg
141    else:
142        try:
143            starttime = float(fields[0])
144        except Error:
145            msg = "Bad time format"
146            raise DataTimeError, msg
147
148    # Split values
149    values = []
150    for value in fields[1].split():
151        values.append(float(value))
152
153    q = ensure_numeric(values)
154
155    msg = 'ERROR: File must contain at least one independent value'
156    assert len(q.shape) == 1, msg
157
158    # Read times proper
159    from anuga.config import time_format
160    import time, calendar
161
162    fid = open(file_text)
163    lines = fid.readlines()
164    fid.close()
165
166    N = len(lines)
167    d = len(q)
168
169    T = num.zeros(N, num.float)       # Time
170    Q = num.zeros((N, d), num.float)  # Values
171
172    for i, line in enumerate(lines):
173        fields = line.split(',')
174        if not time_as_seconds:
175            realtime = calendar.timegm(time.strptime(fields[0], time_format))
176        else:
177            realtime = float(fields[0])
178        T[i] = realtime - starttime
179
180        for j, value in enumerate(fields[1].split()):
181            Q[i, j] = float(value)
182
183    msg = 'File %s must list time as a monotonuosly ' % filename
184    msg += 'increasing sequence'
185    assert num.alltrue(T[1:] - T[:-1] > 0), msg
186
187    #Create NetCDF file
188    from Scientific.IO.NetCDF import NetCDFFile
189
190    fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
191
192    fid.institution = 'Geoscience Australia'
193    fid.description = 'Time series'
194
195    #Reference point
196    #Start time in seconds since the epoch (midnight 1/1/1970)
197    #FIXME: Use Georef
198    fid.starttime = starttime
199
200    # dimension definitions
201    #fid.createDimension('number_of_volumes', self.number_of_volumes)
202    #fid.createDimension('number_of_vertices', 3)
203
204    fid.createDimension('number_of_timesteps', len(T))
205
206    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
207
208    fid.variables['time'][:] = T
209
210    for i in range(Q.shape[1]):
211        try:
212            name = quantity_names[i]
213        except:
214            name = 'Attribute%d' % i
215
216        fid.createVariable(name, netcdf_float, ('number_of_timesteps',))
217        fid.variables[name][:] = Q[:,i]
218
219    fid.close()
220   
221   
222
223##
224# @brief
225# @param filename
226# @param verbose
227def tsh2sww(filename, verbose=False):
228    """
229    to check if a tsh/msh file 'looks' good.
230    """
231
232    if filename[-4:] != '.tsh' and filename[-4:] != '.msh':
233        raise IOError('Input file %s should be .tsh or .msh.' % name_out)
234
235    if verbose == True: log.critical('Creating domain from %s' % filename)
236
237    domain = pmesh_to_domain_instance(filename, Domain)
238
239    if verbose == True: log.critical("Number of triangles = %s" % len(domain))
240
241    domain.smooth = True
242    domain.format = 'sww'   #Native netcdf visualisation format
243    file_path, filename = path.split(filename)
244    filename, ext = path.splitext(filename)
245    domain.set_name(filename)
246    domain.reduction = mean
247
248    if verbose == True: log.critical("file_path = %s" % file_path)
249
250    if file_path == "":
251        file_path = "."
252    domain.set_datadir(file_path)
253
254    if verbose == True:
255        log.critical("Output written to %s%s%s.%s"
256                     % (domain.get_datadir(), sep, domain.get_name(),
257                        domain.format))
258
259    sww = SWW_file(domain)
260    sww.store_connectivity()
261    sww.store_timestep()
262
263
264
265##
266# @brief Return max&min indexes (for slicing) of area specified.
267# @param latitudes_ref ??
268# @param longitudes_ref ??
269# @param minlat Minimum latitude of specified area.
270# @param maxlat Maximum latitude of specified area.
271# @param minlon Minimum longitude of specified area.
272# @param maxlon Maximum longitude of specified area.
273# @return Tuple (lat_min_index, lat_max_index, lon_min_index, lon_max_index)
274def get_min_max_indices(latitudes_ref, longitudes_ref,
275                         minlat=None, maxlat=None,
276                         minlon=None, maxlon=None):
277    """
278    Return max, min indexes (for slicing) of the lat and long arrays to cover
279    the area specified with min/max lat/long.
280
281    Think of the latitudes and longitudes describing a 2d surface.
282    The area returned is, if possible, just big enough to cover the
283    inputed max/min area. (This will not be possible if the max/min area
284    has a section outside of the latitudes/longitudes area.)
285
286    asset  longitudes are sorted,
287    long - from low to high (west to east, eg 148 - 151)
288    assert latitudes are sorted, ascending or decending
289    """
290
291    latitudes = latitudes_ref[:]
292    longitudes = longitudes_ref[:]
293
294    latitudes = ensure_numeric(latitudes)
295    longitudes = ensure_numeric(longitudes)
296
297    assert num.allclose(num.sort(longitudes), longitudes)
298
299    #print latitudes[0],longitudes[0]
300    #print len(latitudes),len(longitudes)
301    #print latitudes[len(latitudes)-1],longitudes[len(longitudes)-1]
302
303    lat_ascending = True
304    if not num.allclose(num.sort(latitudes), latitudes):
305        lat_ascending = False
306        # reverse order of lat, so it's in ascending order
307        latitudes = latitudes[::-1]
308        assert num.allclose(num.sort(latitudes), latitudes)
309
310    largest_lat_index = len(latitudes)-1
311
312    #Cut out a smaller extent.
313    if minlat == None:
314        lat_min_index = 0
315    else:
316        lat_min_index = num.searchsorted(latitudes, minlat)-1
317        if lat_min_index < 0:
318            lat_min_index = 0
319
320    if maxlat == None:
321        lat_max_index = largest_lat_index #len(latitudes)
322    else:
323        lat_max_index = num.searchsorted(latitudes, maxlat)
324        if lat_max_index > largest_lat_index:
325            lat_max_index = largest_lat_index
326
327    if minlon == None:
328        lon_min_index = 0
329    else:
330        lon_min_index = num.searchsorted(longitudes, minlon)-1
331        if lon_min_index < 0:
332            lon_min_index = 0
333
334    if maxlon == None:
335        lon_max_index = len(longitudes)
336    else:
337        lon_max_index = num.searchsorted(longitudes, maxlon)
338
339    # Reversing the indexes, if the lat array is decending
340    if lat_ascending is False:
341        lat_min_index, lat_max_index = largest_lat_index - lat_max_index, \
342                                       largest_lat_index - lat_min_index
343    lat_max_index = lat_max_index + 1 # taking into account how slicing works
344    lon_max_index = lon_max_index + 1 # taking into account how slicing works
345
346    return lat_min_index, lat_max_index, lon_min_index, lon_max_index
347
348
349##
350# @brief
351# @param filename
352# @param x
353# @param y
354# @param z
355def write_obj(filename, x, y, z):
356    """Store x,y,z vectors into filename (obj format).
357
358       Vectors are assumed to have dimension (M,3) where
359       M corresponds to the number elements.
360       triangles are assumed to be disconnected
361
362       The three numbers in each vector correspond to three vertices,
363
364       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
365    """
366
367    import os.path
368
369    root, ext = os.path.splitext(filename)
370    if ext == '.obj':
371        FN = filename
372    else:
373        FN = filename + '.obj'
374
375    outfile = open(FN, 'wb')
376    outfile.write("# Triangulation as an obj file\n")
377
378    M, N = x.shape
379    assert N == 3  #Assuming three vertices per element
380
381    for i in range(M):
382        for j in range(N):
383            outfile.write("v %f %f %f\n" % (x[i,j], y[i,j], z[i,j]))
384
385    for i in range(M):
386        base = i * N
387        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
388
389    outfile.close()
390
Note: See TracBrowser for help on using the repository browser.