source: inundation/ga/storm_surge/pyvolution/data_manager.py @ 280

Last change on this file since 280 was 280, checked in by ole, 21 years ago

fiddle fiddle
File size: 20.6 KB
Line 
1"""Functions to store and retrieve data for the Shallow Water Wave equation.
2There are two kinds of data
3
4  1: Constant data: Vertex coordinates and field values
5  2: Variable data: Conserved quantities
6
7All data is assumed to reside at vertex locations. 
8
9"""
10
11from Numeric import concatenate       
12
13
14def make_filename(s):
15    """Transform argument string into a suitable filename
16    """
17
18    s = s.strip()
19    s = s.replace(' ', '_')
20    #s = s.replace('(', '{')
21    #s = s.replace(')', '}')
22    s = s.replace('(', '')
23    s = s.replace(')', '')
24    s = s.replace('__', '_')               
25       
26    return s
27
28
29def check_dir(path, verbose=None):
30    """Check that specified path exists.
31    If path does not exist it will be created if possible   
32
33    USAGE:
34       checkdir(path, verbose):
35
36    ARGUMENTS:
37        path -- Directory
38        verbose -- Flag verbose output (default: None)
39
40    RETURN VALUE:
41        Verified path including trailing separator
42       
43    """
44
45    import os, sys
46    import os.path
47
48    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
49        unix = 0
50    else:
51        unix = 1
52
53
54    if path[-1] != os.sep: 
55        path = path + os.sep  # Add separator for directories
56
57    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
58    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
59        try:
60            exitcode=os.mkdir(path)
61
62            # Change access rights if possible
63            #
64            if unix:
65                exitcode=os.system('chmod 775 '+path)
66            else:
67                pass  # FIXME: What about acces rights under Windows?
68           
69            if verbose: print 'MESSAGE: Directory', path, 'created.'
70         
71        except:
72            print 'WARNING: Directory', path, 'could not be created.'
73            if unix:
74                path = '/tmp/'
75            else:
76                path = 'C:'
77               
78            print 'Using directory %s instead' %path
79
80    return(path)
81
82   
83
84def del_dir(path):
85    """Recursively delete directory path and all its contents
86    """
87
88    import os
89
90    if os.path.isdir(path):
91        for file in os.listdir(path):
92            X = os.path.join(path, file)
93
94
95            if os.path.isdir(X) and not os.path.islink(X):
96                del_dir(X) 
97            else:
98                try:
99                    os.remove(X)
100                except:
101                    print "Could not remove file %s" %X
102
103        os.rmdir(path)
104
105
106       
107def create_filename(filename, format, size, time=None):
108
109    import os
110    from config import data_dir
111
112           
113    FN = check_dir(data_dir) + filename + '_size%d' %size
114       
115    if time is not None:
116        FN += '_time%.2f' %time
117       
118    FN += '.' + format
119       
120    return FN
121   
122
123def get_files(filename, format, size):
124    """Get all file (names) with gven name, size and format
125    """
126
127    import glob
128
129    import os
130    from config import data_dir
131
132           
133    dir = check_dir(data_dir)
134
135    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
136    return glob.glob(pattern)
137
138#Native checkpoint format.
139#Information needed to recreate a state is preserved
140def cpt_variable_writer(filename, t, v0, v1, v2):
141    """Store all conserved quantities to file
142    """
143
144    M, N = v0.shape
145       
146    FN = create_filename(filename, 'cpt', M, t)
147    #print 'Writing to %s' %FN
148   
149    fid = open(FN, 'w')
150    for i in range(M):
151        for j in range(N):       
152            fid.write('%.16e ' %v0[i,j])           
153        for j in range(N):
154            fid.write('%.16e ' %v1[i,j])                       
155        for j in range(N):
156            fid.write('%.16e ' %v2[i,j])                       
157           
158        fid.write('\n')
159    fid.close()
160
161
162def cpt_variable_reader(filename, t, v0, v1, v2):
163    """Store all conserved quantities to file
164    """
165
166    M, N = v0.shape
167       
168    FN = create_filename(filename, 'cpt', M, t)
169    #print 'Reading from %s' %FN
170   
171    fid = open(FN)
172
173
174    for i in range(M):
175        values = fid.readline().split() #Get one line
176
177        for j in range(N):
178            v0[i,j] = float(values[j])
179            v1[i,j] = float(values[3+j])
180            v2[i,j] = float(values[6+j])
181   
182    fid.close()
183
184def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
185    """Writes x,y,z,z,z coordinates of triangles constituting the bed
186    elevation.
187    Not in use pt
188    """
189
190    M, N = v0.shape
191
192    print X0
193    import sys; sys.exit() 
194    FN = create_filename(filename, 'cpt', M)
195    print 'Writing to %s' %FN
196   
197    fid = open(FN, 'w')
198    for i in range(M):
199        for j in range(2):       
200            fid.write('%.16e ' %X0[i,j])   #x, y
201        for j in range(N):                   
202            fid.write('%.16e ' %v0[i,j])       #z,z,z,
203           
204        for j in range(2):                   
205            fid.write('%.16e ' %X1[i,j])   #x, y
206        for j in range(N):                               
207            fid.write('%.16e ' %v1[i,j]) 
208       
209        for j in range(2):                               
210            fid.write('%.16e ' %X2[i,j])   #x, y
211        for j in range(N):                                           
212            fid.write('%.16e ' %v2[i,j]) 
213           
214        fid.write('\n')
215    fid.close()
216
217
218
219#Function for storing out to e.g. visualisation
220def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
221    """Writes x,y,z coordinates of triangles constituting the bed elevation.
222    """
223
224    M, N = v0.shape
225       
226    FN = create_filename(filename, 'dat', M)
227    #print 'Writing to %s' %FN
228   
229    fid = open(FN, 'w')
230    for i in range(M):
231        for j in range(2):       
232            fid.write('%f ' %X0[i,j])   #x, y
233        fid.write('%f ' %v0[i,0])       #z
234           
235        for j in range(2):                   
236            fid.write('%f ' %X1[i,j])   #x, y
237        fid.write('%f ' %v1[i,0])       #z
238       
239        for j in range(2):                               
240            fid.write('%f ' %X2[i,j])   #x, y
241        fid.write('%f ' %v2[i,0])       #z
242           
243        fid.write('\n')
244    fid.close()
245
246
247
248def dat_variable_writer(filename, t, v0, v1, v2):
249    """Store water height to file
250    """
251
252    M, N = v0.shape
253       
254    FN = create_filename(filename, 'dat', M, t)
255    #print 'Writing to %s' %FN
256   
257    fid = open(FN, 'w')
258    for i in range(M):
259        fid.write('%.4f ' %v0[i,0])
260        fid.write('%.4f ' %v1[i,0])
261        fid.write('%.4f ' %v2[i,0])       
262
263        fid.write('\n')
264    fid.close()
265
266
267#Function for storing out to e.g. visualisation
268class Data_format:
269    """Generic interface to data formats
270    """
271
272       
273    def __init__(self, domain, extension, mode = 'w'):
274        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
275                                        '''   'w' (write)'''+\
276                                        '''   'r' (read)''' +\
277                                        '''   'a' (append)'''   
278
279        #Create filename
280        #FIXME: This may change if separate timesteps are stored in separate files.
281        self.filename = create_filename(domain.get_name(), extension, len(domain))
282        self.timestep = 0
283        self.number_of_volumes = len(domain)
284        self.domain = domain
285   
286   
287
288    def store_all(self):
289        """Writes x,y,z coordinates of triangles constituting
290        the bed elevation.
291        """
292
293        msg = 'Method should be overridden in concrete class'
294        raise msg
295
296   
297    def store_timestep(self, t, V0, V1, V2):
298        """Store time, water heights (and momentums) to file
299        """
300
301        msg = 'Method should be overridden in concrete class'
302        raise msg       
303   
304    #FIXME: Should we have a general set_precision function?
305   
306
307#Function for storing out to e.g. visualisation
308class Data_format_sww(Data_format):
309    """Interface to native NetCDF format (.sww)
310    """
311
312       
313    def __init__(self, domain, mode = 'w'):
314        from Scientific.IO.NetCDF import NetCDFFile
315        from Numeric import Int, Float, Float32
316
317        self.precision = Float32 #Use single precision
318       
319        Data_format.__init__(self, domain, 'sww', mode)
320
321
322        # NetCDF file definition
323        fid = NetCDFFile(self.filename, mode)
324       
325        if mode == 'w':
326           
327            #Create new file
328            fid.institution = 'Geoscience Australia'
329            if domain.smooth:
330                fid.smoothing = 'Yes'
331            else:   
332                fid.smoothing = 'No'
333
334            fid.order = domain.default_order   
335
336            # dimension definitions
337            fid.createDimension('number_of_volumes', self.number_of_volumes) 
338            fid.createDimension('number_of_vertices', 3)
339
340            if domain.smooth is True:
341                fid.createDimension('number_of_points', len(domain.vertexlist))
342            else:   
343                fid.createDimension('number_of_points', 3*self.number_of_volumes)
344               
345            fid.createDimension('number_of_timesteps', None) #extensible
346
347            # variable definitions
348            fid.createVariable('x', self.precision, ('number_of_points',))
349            fid.createVariable('y', self.precision, ('number_of_points',))
350            fid.createVariable('z', self.precision, ('number_of_points',))
351       
352            fid.createVariable('volumes', Int, ('number_of_volumes',
353                                                'number_of_vertices'))
354
355            fid.createVariable('time', self.precision,
356                               ('number_of_timesteps',))
357       
358            fid.createVariable('stage', self.precision,
359                               ('number_of_timesteps',
360                                'number_of_points'))
361
362
363        #Close
364        fid.close()               
365
366    def store_all(self):
367        """Specialisation of store all for net CDF format
368
369        Writes x,y,z coordinates of triangles constituting
370        the bed elevation.
371        """
372
373        from Scientific.IO.NetCDF import NetCDFFile
374
375        from Numeric import concatenate
376
377        domain = self.domain
378
379        #Get NetCDF
380        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
381       
382        # Get the variables
383        x = fid.variables['x']
384        y = fid.variables['y']
385        z = fid.variables['z']
386       
387        volumes = fid.variables['volumes']
388
389        # Get X, Y and bed elevation (index=0)
390        X,Y,A,V = domain.get_vertex_values(xy=True,
391                                           value_array='field_values',
392                                           indices = (0,),
393                                           precision = self.precision)
394
395
396        x[:] = X
397        y[:] = Y
398        z[:] = A[:,0]
399
400        volumes[:] = V
401       
402        #Close
403        fid.close()
404
405
406    def store_timestep(self):
407        """Store time, water heights (and momentums) to file
408        """
409        from Scientific.IO.NetCDF import NetCDFFile
410        from time import sleep
411               
412        #Get NetCDF
413        retries = 0
414        file_open = False
415        while not file_open and retries < 10:
416            try:       
417                fid = NetCDFFile(self.filename, 'a') #Open existing file
418            except IOError:
419                #This could happen if someone was reading the file.
420                #In that case, wait a while and try again
421                msg = 'Warning (store_timestep): File %s could not be opened'\
422                %self.filename
423                msg += ' - trying again'
424                print msg
425                retries += 1
426                sleep(1)
427            else:
428               file_open = True
429               
430        if not file_open:
431            msg = 'File %s could not be opened for append' %self.filename
432            raise msg
433               
434           
435        domain = self.domain
436       
437        # Get the variables
438        time = fid.variables['time']
439        stage = fid.variables['stage']       
440        i = len(time)
441
442        #Store stage
443        time[i] = self.domain.time
444
445        # Get stage (index = 0)
446        A,V = domain.get_vertex_values(xy=False,
447                                       value_array='conserved_quantities',
448                                       indices = (0,),
449                                       precision = self.precision)
450
451       
452        stage[i,:] = A[:,0]  #Only store stage
453
454        #Flush and close
455        fid.sync()
456        fid.close()
457
458
459#Function for storing out to e.g. visualisation
460class Data_format_xya(Data_format):
461    """Generic interface to data formats
462    """
463
464
465    def __init__(self, domain, mode = 'w'):
466        from Scientific.IO.NetCDF import NetCDFFile
467        from Numeric import Int, Float, Float32
468
469        self.precision = Float32 #Use single precision
470       
471        Data_format.__init__(self, domain, 'xya', mode)
472       
473   
474   
475
476    def store_all(self):
477        """Specialisation of store all for xya format
478
479        Writes x,y,z coordinates of triangles constituting
480        the bed elevation.
481        """
482
483        from Numeric import concatenate
484
485        domain = self.domain
486       
487        fd = open(self.filename, 'w')
488
489
490        if domain.smooth is True:
491            number_of_points =  len(domain.vertexlist)
492        else:   
493            number_of_points = 3*self.number_of_volumes
494
495        numVertAttrib = 3 #Three attributes is what is assumed by the xya format
496
497        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
498                 " # <vertex #> <x> <y> [attributes]" + "\n")
499
500
501        # Get X, Y, bed elevation and friction (index=0,1)
502        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
503                                           indices = (0,1), precision = self.precision)
504
505        bed_eles = A[:,0]
506        fricts = A[:,1]
507
508        # Get stage (index=0)
509        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
510                                       indices = (0,), precision = self.precision)
511
512        stages = B[:,0]
513
514        #<vertex #> <x> <y> [attributes]
515        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
516                                               stages, fricts):
517
518            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
519            fd.write(s)
520           
521        #close
522        fd.close()
523
524
525       
526
527
528    def store_timestep(self, t, V0, V1, V2):
529        """Store time, water heights (and momentums) to file
530        """
531        pass
532   
533
534#Auxiliary
535def write_obj(filename,x,y,z):
536    """Store x,y,z vectors into filename (obj format)
537       Vectors are assumed to have dimension (M,3) where
538       M corresponds to the number elements.
539       triangles are assumed to be disconnected
540
541       The three numbers in each vector correspond to three vertices,
542
543       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
544       
545    """
546    #print 'Writing obj to %s' % filename
547
548    import os.path
549
550    root, ext = os.path.splitext(filename)
551    if ext == '.obj':
552        FN = filename
553    else:
554        FN = filename + '.obj'
555       
556
557    outfile = open(FN, 'wb')
558    outfile.write("# Triangulation as an obj file\n")
559
560    M, N = x.shape
561    assert N==3  #Assuming three vertices per element
562
563    for i in range(M):
564        for j in range(N):
565            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
566
567    for i in range(M):
568        base = i*N
569        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
570
571    outfile.close()
572
573
574
575#Conversion routines
576def sww2obj(basefilename, size):
577    """Convert netcdf based data output to obj
578    """
579    from Scientific.IO.NetCDF import NetCDFFile       
580
581    from Numeric import Float, zeros
582
583    #Get NetCDF
584    FN = create_filename(basefilename, 'sww', size)
585    print 'Reading from ', FN
586    fid = NetCDFFile(FN, 'r')  #Open existing file for read
587
588       
589    # Get the variables
590    x = fid.variables['x']
591    y = fid.variables['y']
592    z = fid.variables['z']       
593    time = fid.variables['time']
594    stage = fid.variables['stage']
595       
596    M = size  #Number of lines
597    xx = zeros((M,3), Float)
598    yy = zeros((M,3), Float)
599    zz = zeros((M,3), Float)
600
601    for i in range(M):
602        for j in range(3):
603            xx[i,j] = x[i+j*M]
604            yy[i,j] = y[i+j*M]
605            zz[i,j] = z[i+j*M]
606       
607    #Write obj for bathymetry
608    FN = create_filename(basefilename, 'obj', size)       
609    write_obj(FN,xx,yy,zz)
610
611
612    #Now read all the data with variable information, combine with
613    #x,y info and store as obj
614
615    for k in range(len(time)):
616        t = time[k]
617        print 'Processing timestep %f' %t
618
619        for i in range(M):
620            for j in range(3):
621                zz[i,j] = stage[k,i+j*M]
622
623
624        #Write obj for variable data
625        #FN = create_filename(basefilename, 'obj', size, time=t)
626        FN = create_filename(basefilename[:5], 'obj', size, time=t)   
627        write_obj(FN,xx,yy,zz)
628
629
630def dat2obj(basefilename):
631    """Convert line based data output to obj
632    """
633
634    import glob, os
635    from config import data_dir
636   
637   
638    #Get bathymetry and x,y's
639    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
640   
641    from Numeric import zeros, Float
642
643    M = len(lines)  #Number of lines
644    x = zeros((M,3), Float)
645    y = zeros((M,3), Float)
646    z = zeros((M,3), Float)
647
648    ##i = 0
649    for i, line in enumerate(lines):
650        tokens = line.split()
651        values = map(float,tokens)       
652
653        for j in range(3):
654            x[i,j] = values[j*3]
655            y[i,j] = values[j*3+1]
656            z[i,j] = values[j*3+2]       
657       
658        ##i += 1
659
660
661    #Write obj for bathymetry
662    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
663
664
665    #Now read all the data files with variable information, combine with
666    #x,y info
667    #and store as obj
668
669    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
670
671    for filename in files:
672        print 'Processing %s' % filename
673
674        lines = open(data_dir+os.sep+filename,'r').readlines()
675        assert len(lines) == M
676        root, ext = os.path.splitext(filename)
677
678        #Get time from filename
679        i0 = filename.find('_time=')
680        if i0 == -1:
681            #Skip bathymetry file
682            continue
683       
684        i0 += 6  #Position where time starts           
685        i1 = filename.find('.dat')         
686
687        if i1 > i0:
688            t = float(filename[i0:i1])
689        else:
690            raise 'Hmmmm'
691               
692       
693       
694        ##i = 0
695        for i, line in enumerate(lines):
696            tokens = line.split()
697            values = map(float,tokens)       
698
699            for j in range(3):
700                z[i,j] = values[j]       
701       
702            ##i += 1
703
704        #Write obj for variable data
705        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
706
707
708def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
709    """Read netcdf filename1, pick timesteps first:step:last and save to
710    nettcdf file filename2
711    """
712    from Scientific.IO.NetCDF import NetCDFFile       
713   
714    #Get NetCDF
715    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
716    outfile = NetCDFFile(filename2, 'w')  #Open new file
717   
718
719    #Copy dimensions
720    for d in infile.dimensions:
721        outfile.createDimension(d, infile.dimensions[d])
722
723    for name in infile.variables:
724        var = infile.variables[name]
725        outfile.createVariable(name, var.typecode(), var.dimensions)
726       
727   
728    #Copy the static variables
729    for name in infile.variables:
730        if name == 'time' or name == 'stage':
731            pass
732        else:
733            #Copy
734            outfile.variables[name][:] = infile.variables[name][:]           
735
736    #Copy selected timesteps
737    time = infile.variables['time']
738    stage = infile.variables['stage']
739
740    newtime = outfile.variables['time']
741    newstage = outfile.variables['stage']   
742
743    if last is None:
744        last = len(time)
745
746    selection = range(first, last, step)
747    for i, j in enumerate(selection):
748        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
749        newtime[i] = time[j]
750        newstage[i,:] = stage[j,:]       
751       
752    #Close
753    infile.close()
754    outfile.close()
755
756
757#Get data objects
758def get_dataobject(domain, mode='w'):
759    """Return instance of class of given format using filename
760    """
761
762    cls = eval('Data_format_%s' %domain.format)
763    return cls(domain, mode)
764
765
766
767
Note: See TracBrowser for help on using the repository browser.