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

Last change on this file since 1097 was 1090, checked in by prow, 20 years ago

Removing tracing.
More testing.
File size: 54.1 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. Stored once
5  2: Variable data: Conserved quantities. Stored once per timestep.
6
7All data is assumed to reside at vertex locations.
8
9
10Formats used within ANUGA:
11
12.sww: Netcdf format for storing model output
13
14.xya: ASCII format for storing arbitrary points and associated attributes
15.pts: NetCDF format for storing arbitrary points and associated attributes
16
17.asc: ASCII foramt of regular DEMs as output from ArcView
18.dem: NetCDF representation of regular DEM data
19
20.tsh: ASCII format for storing meshes and associated boundary and region info
21
22.nc: Native ferret NetCDF format
23
24A typical dataflow can be described as follows
25
26Manually created files:
27ASC, PRJ:     Digital elevation models (gridded)
28TSH:          Triangular meshes (e.g. dreated from pmesh)
29NC            Model outputs for use as boundary conditions (e.g from MOST)
30
31
32AUTOMATICALLY CREATED FILES:
33
34ASC, PRJ  ->  DEM  ->  PTS: Conversion of DEM's to native pts file
35
36NC -> SWW: Conversion of MOST bundary files to boundary sww
37
38PTS + TSH -> TSH with elevation: Least squares fit
39
40TSH -> SWW:  Conversion of TSH to sww viewable using Swollen
41
42TSH + Boundary SWW -> SWW: SImluation using pyvolution
43
44
45"""
46
47from Numeric import concatenate
48
49
50def make_filename(s):
51    """Transform argument string into a suitable filename
52    """
53
54    s = s.strip()
55    s = s.replace(' ', '_')
56    s = s.replace('(', '')
57    s = s.replace(')', '')
58    s = s.replace('__', '_')
59
60    return s
61
62
63def check_dir(path, verbose=None):
64    """Check that specified path exists.
65    If path does not exist it will be created if possible
66
67    USAGE:
68       checkdir(path, verbose):
69
70    ARGUMENTS:
71        path -- Directory
72        verbose -- Flag verbose output (default: None)
73
74    RETURN VALUE:
75        Verified path including trailing separator
76
77    """
78
79    import os, sys
80    import os.path
81
82    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
83        unix = 0
84    else:
85        unix = 1
86
87
88    if path[-1] != os.sep:
89        path = path + os.sep  # Add separator for directories
90
91    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
92    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
93        try:
94            exitcode=os.mkdir(path)
95
96            # Change access rights if possible
97            #
98            if unix:
99                exitcode=os.system('chmod 775 '+path)
100            else:
101                pass  # FIXME: What about acces rights under Windows?
102
103            if verbose: print 'MESSAGE: Directory', path, 'created.'
104
105        except:
106            print 'WARNING: Directory', path, 'could not be created.'
107            if unix:
108                path = '/tmp/'
109            else:
110                path = 'C:'
111
112            print 'Using directory %s instead' %path
113
114    return(path)
115
116
117
118def del_dir(path):
119    """Recursively delete directory path and all its contents
120    """
121
122    import os
123
124    if os.path.isdir(path):
125        for file in os.listdir(path):
126            X = os.path.join(path, file)
127
128
129            if os.path.isdir(X) and not os.path.islink(X):
130                del_dir(X)
131            else:
132                try:
133                    os.remove(X)
134                except:
135                    print "Could not remove file %s" %X
136
137        os.rmdir(path)
138
139
140
141def create_filename(datadir, filename, format, size=None, time=None):
142
143    import os
144    #from config import data_dir
145
146    FN = check_dir(datadir) + filename
147
148    if size is not None:
149        FN += '_size%d' %size
150
151    if time is not None:
152        FN += '_time%.2f' %time
153
154    FN += '.' + format
155    return FN
156
157
158def get_files(datadir, filename, format, size):
159    """Get all file (names) with gven name, size and format
160    """
161
162    import glob
163
164    import os
165    #from config import data_dir
166
167    dir = check_dir(datadir)
168
169    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
170    return glob.glob(pattern)
171
172
173
174#Generic class for storing output to e.g. visualisation or checkpointing
175class Data_format:
176    """Generic interface to data formats
177    """
178
179
180    def __init__(self, domain, extension, mode = 'w'):
181        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
182                                        '''   'w' (write)'''+\
183                                        '''   'r' (read)''' +\
184                                        '''   'a' (append)'''
185
186        #Create filename
187        #self.filename = create_filename(domain.get_datadir(),
188        #domain.get_name(), extension, len(domain))
189
190
191        self.filename = create_filename(domain.get_datadir(),
192                        domain.get_name(), extension)
193
194        #print 'F', self.filename
195        self.timestep = 0
196        self.number_of_volumes = len(domain)
197        self.domain = domain
198
199
200        #FIXME: Should we have a general set_precision function?
201
202
203
204#Class for storing output to e.g. visualisation
205class Data_format_sww(Data_format):
206    """Interface to native NetCDF format (.sww)
207    """
208
209
210    def __init__(self, domain, mode = 'w',\
211                 max_size = 2000000000,recursion=False):
212        from Scientific.IO.NetCDF import NetCDFFile
213        from Numeric import Int, Float, Float32
214
215        self.precision = Float32 #Use single precision
216        if hasattr(domain,'max_size'):
217            self.max_size = domain.max_size#file size max is 2Gig
218        else:
219            self.max_size = max_size
220        self.recursion = recursion
221        self.mode = mode
222
223        Data_format.__init__(self, domain, 'sww', mode)
224
225
226        # NetCDF file definition
227        fid = NetCDFFile(self.filename, mode)
228
229        if mode == 'w':
230
231            #Create new file
232            fid.institution = 'Geoscience Australia'
233            fid.description = 'Output from pyvolution suitable for plotting'
234
235            if domain.smooth:
236                fid.smoothing = 'Yes'
237            else:
238                fid.smoothing = 'No'
239
240            fid.order = domain.default_order
241
242            #Reference point
243            #Start time in seconds since the epoch (midnight 1/1/1970)
244            fid.starttime = domain.starttime
245            fid.xllcorner = domain.xllcorner
246            fid.yllcorner = domain.yllcorner
247            fid.zone = str(domain.zone) #FIXME: ?
248
249            # dimension definitions
250            fid.createDimension('number_of_volumes', self.number_of_volumes)
251            fid.createDimension('number_of_vertices', 3)
252
253            if domain.smooth is True:
254                fid.createDimension('number_of_points', len(domain.vertexlist))
255            else:
256                fid.createDimension('number_of_points', 3*self.number_of_volumes)
257
258            fid.createDimension('number_of_timesteps', None) #extensible
259
260            # variable definitions
261            fid.createVariable('x', self.precision, ('number_of_points',))
262            fid.createVariable('y', self.precision, ('number_of_points',))
263            fid.createVariable('elevation', self.precision, ('number_of_points',))
264
265            #FIXME: Backwards compatibility
266            fid.createVariable('z', self.precision, ('number_of_points',))
267            #################################
268
269            fid.createVariable('volumes', Int, ('number_of_volumes',
270                                                'number_of_vertices'))
271
272            fid.createVariable('time', self.precision,
273                               ('number_of_timesteps',))
274
275            fid.createVariable('stage', self.precision,
276                               ('number_of_timesteps',
277                                'number_of_points'))
278
279            fid.createVariable('xmomentum', self.precision,
280                               ('number_of_timesteps',
281                                'number_of_points'))
282
283            fid.createVariable('ymomentum', self.precision,
284                               ('number_of_timesteps',
285                                'number_of_points'))
286
287        #Close
288        fid.close()
289
290
291    def store_connectivity(self):
292        """Specialisation of store_connectivity for net CDF format
293
294        Writes x,y,z coordinates of triangles constituting
295        the bed elevation.
296        """
297
298        from Scientific.IO.NetCDF import NetCDFFile
299
300        from Numeric import concatenate, Int
301
302        domain = self.domain
303
304        #Get NetCDF
305        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
306
307        # Get the variables
308        x = fid.variables['x']
309        y = fid.variables['y']
310        z = fid.variables['elevation']
311
312        volumes = fid.variables['volumes']
313
314        # Get X, Y and bed elevation Z
315        Q = domain.quantities['elevation']
316        X,Y,Z,V = Q.get_vertex_values(xy=True,
317                                      precision = self.precision)
318
319
320
321        x[:] = X.astype(self.precision)
322        y[:] = Y.astype(self.precision)
323        z[:] = Z.astype(self.precision)
324
325        #FIXME: Backwards compatibility
326        z = fid.variables['z']
327        z[:] = Z.astype(self.precision)
328        ################################
329
330        volumes[:] = V.astype(volumes.typecode())
331
332        #Close
333        fid.close()
334
335    def store_timestep(self, names):
336        """Store time and named quantities to file
337        """
338        from Scientific.IO.NetCDF import NetCDFFile
339        import types
340        from time import sleep
341        from os import stat
342
343
344        #Get NetCDF
345        retries = 0
346        file_open = False
347        while not file_open and retries < 10:
348            try:
349                fid = NetCDFFile(self.filename, 'a') #Open existing file
350            except IOError:
351                #This could happen if someone was reading the file.
352                #In that case, wait a while and try again
353                msg = 'Warning (store_timestep): File %s could not be opened'\
354                %self.filename
355                msg += ' - trying again'
356                print msg
357                retries += 1
358                sleep(1)
359            else:
360               file_open = True
361
362        if not file_open:
363            msg = 'File %s could not be opened for append' %self.filename
364            raise msg
365
366
367
368        #Check to see if the file is already too big:
369        time = fid.variables['time']
370        i = len(time)+1
371        file_size = stat(self.filename)[6]
372        file_size_increase =  file_size/i
373        if file_size + file_size_increase > self.max_size*(2**self.recursion):
374            #in order to get the file name and start time correct,
375            #I change the domian.filename and domain.starttime.
376            #This is the only way to do this without changing
377            #other modules (I think).
378
379            #write a filename addon that won't break swollens reader
380            #(10.sww is bad)
381            filename_ext = '_time_%s'%self.domain.time           
382            filename_ext = filename_ext.replace('.', '_')
383            #remember the old filename, then give domain a
384            #name with the extension
385            old_domain_filename = self.domain.filename
386            if not self.recursion:
387                self.domain.filename = self.domain.filename+filename_ext
388
389            #change the domain starttime to the current time
390            old_domain_starttime = self.domain.starttime
391            self.domain.starttime = self.domain.time
392
393            #build a new data_structure.
394            next_data_structure=\
395                Data_format_sww(self.domain, mode=self.mode,\
396                                max_size = self.max_size,\
397                                recursion = self.recursion+1)
398            if not self.recursion:
399                print '    file_size = %s'%file_size
400                print '    saving file to %s'%next_data_structure.filename
401            #set up the new data_structure
402            self.domain.writer = next_data_structure
403
404            #FIXME - could be cleaner to use domain.store_timestep etc.
405            next_data_structure.store_connectivity()
406            next_data_structure.store_timestep(names)
407            fid.sync()
408            fid.close()
409
410            #restore the old starttime and filename
411            self.domain.starttime = old_domain_starttime
412            self.domain.filename = old_domain_filename
413        else:
414            self.recursion = False
415            domain = self.domain
416
417            # Get the variables
418            time = fid.variables['time']
419            stage = fid.variables['stage']
420            xmomentum = fid.variables['xmomentum']
421            ymomentum = fid.variables['ymomentum']
422            i = len(time)
423
424            #Store time
425            time[i] = self.domain.time
426
427
428            if type(names) not in [types.ListType, types.TupleType]:
429                names = [names]
430
431            for name in names:
432                # Get quantity
433                Q = domain.quantities[name]
434                A,V = Q.get_vertex_values(xy=False,
435                                      precision = self.precision)
436
437                #FIXME: Make this general (see below)
438                if name == 'stage':
439                    stage[i,:] = A.astype(self.precision)
440                elif name == 'xmomentum':
441                    xmomentum[i,:] = A.astype(self.precision)
442                elif name == 'ymomentum':
443                    ymomentum[i,:] = A.astype(self.precision)
444
445                #As in....
446                #eval( name + '[i,:] = A.astype(self.precision)' )
447                #FIXME: But we need a UNIT test for that before refactoring
448
449
450
451            #Flush and close
452            fid.sync()
453            fid.close()
454
455
456
457#Class for handling checkpoints data
458class Data_format_cpt(Data_format):
459    """Interface to native NetCDF format (.cpt)
460    """
461
462
463    def __init__(self, domain, mode = 'w'):
464        from Scientific.IO.NetCDF import NetCDFFile
465        from Numeric import Int, Float, Float
466
467        self.precision = Float #Use full precision
468
469        Data_format.__init__(self, domain, 'sww', mode)
470
471
472        # NetCDF file definition
473        fid = NetCDFFile(self.filename, mode)
474
475        if mode == 'w':
476            #Create new file
477            fid.institution = 'Geoscience Australia'
478            fid.description = 'Checkpoint data'
479            #fid.smooth = domain.smooth
480            fid.order = domain.default_order
481
482            # dimension definitions
483            fid.createDimension('number_of_volumes', self.number_of_volumes)
484            fid.createDimension('number_of_vertices', 3)
485
486            #Store info at all vertices (no smoothing)
487            fid.createDimension('number_of_points', 3*self.number_of_volumes)
488            fid.createDimension('number_of_timesteps', None) #extensible
489
490            # variable definitions
491
492            #Mesh
493            fid.createVariable('x', self.precision, ('number_of_points',))
494            fid.createVariable('y', self.precision, ('number_of_points',))
495
496
497            fid.createVariable('volumes', Int, ('number_of_volumes',
498                                                'number_of_vertices'))
499
500            fid.createVariable('time', self.precision,
501                               ('number_of_timesteps',))
502
503            #Allocate space for all quantities
504            for name in domain.quantities.keys():
505                fid.createVariable(name, self.precision,
506                                   ('number_of_timesteps',
507                                    'number_of_points'))
508
509        #Close
510        fid.close()
511
512
513    def store_checkpoint(self):
514        """
515        Write x,y coordinates of triangles.
516        Write connectivity (
517        constituting
518        the bed elevation.
519        """
520
521        from Scientific.IO.NetCDF import NetCDFFile
522
523        from Numeric import concatenate
524
525        domain = self.domain
526
527        #Get NetCDF
528        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
529
530        # Get the variables
531        x = fid.variables['x']
532        y = fid.variables['y']
533
534        volumes = fid.variables['volumes']
535
536        # Get X, Y and bed elevation Z
537        Q = domain.quantities['elevation']
538        X,Y,Z,V = Q.get_vertex_values(xy=True,
539                                      precision = self.precision)
540
541
542
543        x[:] = X.astype(self.precision)
544        y[:] = Y.astype(self.precision)
545        z[:] = Z.astype(self.precision)
546
547        volumes[:] = V
548
549        #Close
550        fid.close()
551
552
553    def store_timestep(self, name):
554        """Store time and named quantity to file
555        """
556        from Scientific.IO.NetCDF import NetCDFFile
557        from time import sleep
558
559        #Get NetCDF
560        retries = 0
561        file_open = False
562        while not file_open and retries < 10:
563            try:
564                fid = NetCDFFile(self.filename, 'a') #Open existing file
565            except IOError:
566                #This could happen if someone was reading the file.
567                #In that case, wait a while and try again
568                msg = 'Warning (store_timestep): File %s could not be opened'\
569                %self.filename
570                msg += ' - trying again'
571                print msg
572                retries += 1
573                sleep(1)
574            else:
575               file_open = True
576
577        if not file_open:
578            msg = 'File %s could not be opened for append' %self.filename
579            raise msg
580
581
582        domain = self.domain
583
584        # Get the variables
585        time = fid.variables['time']
586        stage = fid.variables['stage']
587        i = len(time)
588
589        #Store stage
590        time[i] = self.domain.time
591
592        # Get quantity
593        Q = domain.quantities[name]
594        A,V = Q.get_vertex_values(xy=False,
595                                  precision = self.precision)
596
597        stage[i,:] = A.astype(self.precision)
598
599        #Flush and close
600        fid.sync()
601        fid.close()
602
603
604
605
606
607#Function for storing xya output
608#FIXME Not done yet for this version
609class Data_format_xya(Data_format):
610    """Generic interface to data formats
611    """
612
613
614    def __init__(self, domain, mode = 'w'):
615        from Scientific.IO.NetCDF import NetCDFFile
616        from Numeric import Int, Float, Float32
617
618        self.precision = Float32 #Use single precision
619
620        Data_format.__init__(self, domain, 'xya', mode)
621
622
623
624    #FIXME -This is the old xya format
625    def store_all(self):
626        """Specialisation of store all for xya format
627
628        Writes x,y,z coordinates of triangles constituting
629        the bed elevation.
630        """
631
632        from Numeric import concatenate
633
634        domain = self.domain
635
636        fd = open(self.filename, 'w')
637
638
639        if domain.smooth is True:
640            number_of_points =  len(domain.vertexlist)
641        else:
642            number_of_points = 3*self.number_of_volumes
643
644        numVertAttrib = 3 #Three attributes is what is assumed by the xya format
645
646        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
647                 " # <vertex #> <x> <y> [attributes]" + "\n")
648
649
650        # Get X, Y, bed elevation and friction (index=0,1)
651        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
652                                           indices = (0,1), precision = self.precision)
653
654        bed_eles = A[:,0]
655        fricts = A[:,1]
656
657        # Get stage (index=0)
658        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
659                                       indices = (0,), precision = self.precision)
660
661        stages = B[:,0]
662
663        #<vertex #> <x> <y> [attributes]
664        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
665                                               stages, fricts):
666
667            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
668            fd.write(s)
669
670        #close
671        fd.close()
672
673
674    def store_timestep(self, t, V0, V1, V2):
675        """Store time, water heights (and momentums) to file
676        """
677        pass
678
679
680#Auxiliary
681def write_obj(filename,x,y,z):
682    """Store x,y,z vectors into filename (obj format)
683       Vectors are assumed to have dimension (M,3) where
684       M corresponds to the number elements.
685       triangles are assumed to be disconnected
686
687       The three numbers in each vector correspond to three vertices,
688
689       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
690
691    """
692    #print 'Writing obj to %s' % filename
693
694    import os.path
695
696    root, ext = os.path.splitext(filename)
697    if ext == '.obj':
698        FN = filename
699    else:
700        FN = filename + '.obj'
701
702
703    outfile = open(FN, 'wb')
704    outfile.write("# Triangulation as an obj file\n")
705
706    M, N = x.shape
707    assert N==3  #Assuming three vertices per element
708
709    for i in range(M):
710        for j in range(N):
711            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
712
713    for i in range(M):
714        base = i*N
715        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
716
717    outfile.close()
718
719
720
721#Conversion routines
722def sww2obj(basefilename, size):
723    """Convert netcdf based data output to obj
724    """
725    from Scientific.IO.NetCDF import NetCDFFile
726
727    from Numeric import Float, zeros
728
729    #Get NetCDF
730    FN = create_filename('.', basefilename, 'sww', size)
731    print 'Reading from ', FN
732    fid = NetCDFFile(FN, 'r')  #Open existing file for read
733
734
735    # Get the variables
736    x = fid.variables['x']
737    y = fid.variables['y']
738    z = fid.variables['elevation']
739    time = fid.variables['time']
740    stage = fid.variables['stage']
741
742    M = size  #Number of lines
743    xx = zeros((M,3), Float)
744    yy = zeros((M,3), Float)
745    zz = zeros((M,3), Float)
746
747    for i in range(M):
748        for j in range(3):
749            xx[i,j] = x[i+j*M]
750            yy[i,j] = y[i+j*M]
751            zz[i,j] = z[i+j*M]
752
753    #Write obj for bathymetry
754    FN = create_filename('.', basefilename, 'obj', size)
755    write_obj(FN,xx,yy,zz)
756
757
758    #Now read all the data with variable information, combine with
759    #x,y info and store as obj
760
761    for k in range(len(time)):
762        t = time[k]
763        print 'Processing timestep %f' %t
764
765        for i in range(M):
766            for j in range(3):
767                zz[i,j] = stage[k,i+j*M]
768
769
770        #Write obj for variable data
771        #FN = create_filename(basefilename, 'obj', size, time=t)
772        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
773        write_obj(FN,xx,yy,zz)
774
775
776def dat2obj(basefilename):
777    """Convert line based data output to obj
778    FIXME: Obsolete?
779    """
780
781    import glob, os
782    from config import data_dir
783
784
785    #Get bathymetry and x,y's
786    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
787
788    from Numeric import zeros, Float
789
790    M = len(lines)  #Number of lines
791    x = zeros((M,3), Float)
792    y = zeros((M,3), Float)
793    z = zeros((M,3), Float)
794
795    ##i = 0
796    for i, line in enumerate(lines):
797        tokens = line.split()
798        values = map(float,tokens)
799
800        for j in range(3):
801            x[i,j] = values[j*3]
802            y[i,j] = values[j*3+1]
803            z[i,j] = values[j*3+2]
804
805        ##i += 1
806
807
808    #Write obj for bathymetry
809    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
810
811
812    #Now read all the data files with variable information, combine with
813    #x,y info
814    #and store as obj
815
816    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
817
818    for filename in files:
819        print 'Processing %s' % filename
820
821        lines = open(data_dir+os.sep+filename,'r').readlines()
822        assert len(lines) == M
823        root, ext = os.path.splitext(filename)
824
825        #Get time from filename
826        i0 = filename.find('_time=')
827        if i0 == -1:
828            #Skip bathymetry file
829            continue
830
831        i0 += 6  #Position where time starts
832        i1 = filename.find('.dat')
833
834        if i1 > i0:
835            t = float(filename[i0:i1])
836        else:
837            raise 'Hmmmm'
838
839
840
841        ##i = 0
842        for i, line in enumerate(lines):
843            tokens = line.split()
844            values = map(float,tokens)
845
846            for j in range(3):
847                z[i,j] = values[j]
848
849            ##i += 1
850
851        #Write obj for variable data
852        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
853
854
855def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
856    """Read netcdf filename1, pick timesteps first:step:last and save to
857    nettcdf file filename2
858    """
859    from Scientific.IO.NetCDF import NetCDFFile
860
861    #Get NetCDF
862    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
863    outfile = NetCDFFile(filename2, 'w')  #Open new file
864
865
866    #Copy dimensions
867    for d in infile.dimensions:
868        outfile.createDimension(d, infile.dimensions[d])
869
870    for name in infile.variables:
871        var = infile.variables[name]
872        outfile.createVariable(name, var.typecode(), var.dimensions)
873
874
875    #Copy the static variables
876    for name in infile.variables:
877        if name == 'time' or name == 'stage':
878            pass
879        else:
880            #Copy
881            outfile.variables[name][:] = infile.variables[name][:]
882
883    #Copy selected timesteps
884    time = infile.variables['time']
885    stage = infile.variables['stage']
886
887    newtime = outfile.variables['time']
888    newstage = outfile.variables['stage']
889
890    if last is None:
891        last = len(time)
892
893    selection = range(first, last, step)
894    for i, j in enumerate(selection):
895        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
896        newtime[i] = time[j]
897        newstage[i,:] = stage[j,:]
898
899    #Close
900    infile.close()
901    outfile.close()
902
903
904#Get data objects
905def get_dataobject(domain, mode='w'):
906    """Return instance of class of given format using filename
907    """
908
909    cls = eval('Data_format_%s' %domain.format)
910    return cls(domain, mode)
911
912def dem2pts(basename_in, basename_out=None, verbose=False):
913    """Read Digitial Elevation model from the following NetCDF format (.dem)
914
915    Example:
916
917    ncols         3121
918    nrows         1800
919    xllcorner     722000
920    yllcorner     5893000
921    cellsize      25
922    NODATA_value  -9999
923    138.3698 137.4194 136.5062 135.5558 ..........
924
925    Convert to NetCDF pts format which is
926
927    points:  (Nx2) Float array
928    elevation: N Float array
929    """
930
931    import os
932    from Scientific.IO.NetCDF import NetCDFFile
933    from Numeric import Float, arrayrange, concatenate
934
935    root = basename_in
936
937    #Get NetCDF
938    infile = NetCDFFile(root + '.dem', 'r')  #Open existing netcdf file for read
939
940    if verbose: print 'Reading DEM from %s' %(root + '.dem')
941
942    ncols = infile.ncols[0]
943    nrows = infile.nrows[0]
944    xllcorner = infile.xllcorner[0]  #Easting of lower left corner
945    yllcorner = infile.yllcorner[0]  #Northing of lower left corner
946    cellsize = infile.cellsize[0]
947    NODATA_value = infile.NODATA_value[0]
948    dem_elevation = infile.variables['elevation']
949
950    zone = infile.zone[0]
951    false_easting = infile.false_easting[0]
952    false_northing = infile.false_northing[0]
953
954    #Text strings
955    projection = infile.projection
956    datum = infile.datum
957    units = infile.units
958
959
960    #Get output file
961    if basename_out == None:
962        ptsname = root + '.pts'
963    else:
964        ptsname = basename_out + '.pts'
965
966    if verbose: print 'Store to NetCDF file %s' %ptsname
967    # NetCDF file definition
968    outfile = NetCDFFile(ptsname, 'w')
969
970    #Create new file
971    outfile.institution = 'Geoscience Australia'
972    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
973                      'of spatial point data'
974
975    #Georeferencing
976    outfile.zone = zone
977    outfile.xllcorner = xllcorner #Easting of lower left corner
978    outfile.yllcorner = yllcorner #Northing of lower left corner
979    outfile.false_easting = false_easting
980    outfile.false_northing =false_northing
981
982    outfile.projection = projection
983    outfile.datum = datum
984    outfile.units = units
985
986
987    #Grid info (FIXME: probably not going to be used, but heck)
988    outfile.ncols = ncols
989    outfile.nrows = nrows
990
991
992    # dimension definitions
993    outfile.createDimension('number_of_points', nrows*ncols)
994    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
995
996    # variable definitions
997    outfile.createVariable('points', Float, ('number_of_points',
998                                             'number_of_dimensions'))
999    outfile.createVariable('elevation', Float, ('number_of_points',))
1000
1001    # Get handles to the variables
1002    points = outfile.variables['points']
1003    elevation = outfile.variables['elevation']
1004
1005    #Store data
1006    #FIXME: Could perhaps be faster using array operations
1007    for i in range(nrows):
1008        if verbose: print 'Processing row %d of %d' %(i, nrows)
1009
1010        y = (nrows-i)*cellsize
1011        for j in range(ncols):
1012            index = i*ncols + j
1013
1014            x = j*cellsize
1015            points[index, :] = [x,y]
1016            elevation[index] = dem_elevation[i, j]
1017
1018    infile.close()
1019    outfile.close()
1020
1021
1022def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
1023                                  verbose=False):
1024    """Read Digitial Elevation model from the following ASCII format (.asc)
1025
1026    Example:
1027
1028    ncols         3121
1029    nrows         1800
1030    xllcorner     722000
1031    yllcorner     5893000
1032    cellsize      25
1033    NODATA_value  -9999
1034    138.3698 137.4194 136.5062 135.5558 ..........
1035
1036    Convert basename_in + '.asc' to NetCDF format (.dem)
1037    mimicking the ASCII format closely.
1038
1039
1040    An accompanying file with same basename_in but extension .prj must exist
1041    and is used to fix the UTM zone, datum, false northings and eastings.
1042
1043    The prj format is assumed to be as
1044
1045    Projection    UTM
1046    Zone          56
1047    Datum         WGS84
1048    Zunits        NO
1049    Units         METERS
1050    Spheroid      WGS84
1051    Xshift        0.0000000000
1052    Yshift        10000000.0000000000
1053    Parameters
1054    """
1055
1056    import os
1057    from Scientific.IO.NetCDF import NetCDFFile
1058    from Numeric import Float, array
1059
1060    #root, ext = os.path.splitext(basename_in)
1061    root = basename_in
1062
1063    ###########################################
1064    # Read Meta data
1065    if verbose: print 'Reading METADATA from %s' %root + '.prj'
1066    metadatafile = open(root + '.prj')
1067    metalines = metadatafile.readlines()
1068    metadatafile.close()
1069
1070    L = metalines[0].strip().split()
1071    assert L[0].strip().lower() == 'projection'
1072    projection = L[1].strip()                   #TEXT
1073
1074    L = metalines[1].strip().split()
1075    assert L[0].strip().lower() == 'zone'
1076    zone = int(L[1].strip())
1077
1078    L = metalines[2].strip().split()
1079    assert L[0].strip().lower() == 'datum'
1080    datum = L[1].strip()                        #TEXT
1081
1082    L = metalines[3].strip().split()
1083    assert L[0].strip().lower() == 'zunits'     #IGNORE
1084    zunits = L[1].strip()                       #TEXT
1085
1086    L = metalines[4].strip().split()
1087    assert L[0].strip().lower() == 'units'
1088    units = L[1].strip()                        #TEXT
1089
1090    L = metalines[5].strip().split()
1091    assert L[0].strip().lower() == 'spheroid'   #IGNORE
1092    spheroid = L[1].strip()                     #TEXT
1093
1094    L = metalines[6].strip().split()
1095    assert L[0].strip().lower() == 'xshift'
1096    false_easting = float(L[1].strip())
1097
1098    L = metalines[7].strip().split()
1099    assert L[0].strip().lower() == 'yshift'
1100    false_northing = float(L[1].strip())
1101
1102    #print false_easting, false_northing, zone, datum
1103
1104
1105    ###########################################
1106    #Read DEM data
1107
1108    datafile = open(basename_in + '.asc')
1109
1110    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
1111    lines = datafile.readlines()
1112    datafile.close()
1113
1114    if verbose: print 'Got', len(lines), ' lines'
1115
1116    ncols = int(lines[0].split()[1].strip())
1117    nrows = int(lines[1].split()[1].strip())
1118    xllcorner = float(lines[2].split()[1].strip())
1119    yllcorner = float(lines[3].split()[1].strip())
1120    cellsize = float(lines[4].split()[1].strip())
1121    NODATA_value = int(lines[5].split()[1].strip())
1122
1123    assert len(lines) == nrows + 6
1124
1125
1126    ##########################################
1127
1128
1129    if basename_out == None:
1130        netcdfname = root + '.dem'
1131    else:
1132        netcdfname = basename_out + '.dem'
1133
1134    if verbose: print 'Store to NetCDF file %s' %netcdfname
1135    # NetCDF file definition
1136    fid = NetCDFFile(netcdfname, 'w')
1137
1138    #Create new file
1139    fid.institution = 'Geoscience Australia'
1140    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
1141                      'of spatial point data'
1142
1143    fid.ncols = ncols
1144    fid.nrows = nrows
1145    fid.xllcorner = xllcorner
1146    fid.yllcorner = yllcorner
1147    fid.cellsize = cellsize
1148    fid.NODATA_value = NODATA_value
1149
1150    fid.zone = zone
1151    fid.false_easting = false_easting
1152    fid.false_northing = false_northing
1153    fid.projection = projection
1154    fid.datum = datum
1155    fid.units = units
1156
1157
1158    # dimension definitions
1159    fid.createDimension('number_of_rows', nrows)
1160    fid.createDimension('number_of_columns', ncols)
1161
1162    # variable definitions
1163    fid.createVariable('elevation', Float, ('number_of_rows',
1164                                            'number_of_columns'))
1165
1166    # Get handles to the variables
1167    elevation = fid.variables['elevation']
1168
1169    #Store data
1170    for i, line in enumerate(lines[6:]):
1171        fields = line.split()
1172        if verbose: print 'Processing row %d of %d' %(i, nrows)
1173
1174        elevation[i, :] = array([float(x) for x in fields])
1175
1176    fid.close()
1177
1178
1179
1180def ferret2sww(basename_in, basename_out = None,
1181               verbose = False,
1182               minlat = None, maxlat = None,
1183               minlon = None, maxlon = None,
1184               mint = None, maxt = None, mean_stage = 0,
1185               origin = None, zscale = 1, 
1186               elevation = -100): #FIXME: Bathymetry should be obtained
1187                                  #from MOST somehow.
1188                                  #Alternatively from elsewhere
1189                                  #or, as a last resort,
1190                                  #specified here.
1191                                  #The value of -100 will work
1192                                  #for the Wollongong tsunami
1193                                  #scenario but is very hacky
1194    """Convert 'Ferret' NetCDF format for wave propagation to
1195    sww format native to pyvolution.
1196
1197    Specify only basename_in and read files of the form
1198    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1199    relative height, x-velocity and y-velocity, respectively.
1200
1201    Also convert latitude and longitude to UTM. All coordinates are
1202    assumed to be given in the GDA94 datum.
1203
1204    min's and max's: If omitted - full extend is used.
1205    To include a value min may equal it, while max must exceed it.
1206    Lat and lon are assuemd to be in decimal degrees
1207
1208    origin is a 3-tuple with geo referenced
1209    UTM coordinates (zone, easting, northing)
1210
1211    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
1212    which means that longitude is the fastest
1213    varying dimension (row major order, so to speak)
1214
1215    ferret2sww uses grid points as vertices in a triangular grid
1216    counting vertices from lower left corner upwards, then right
1217    """ 
1218
1219    import os
1220    from Scientific.IO.NetCDF import NetCDFFile
1221    from Numeric import Float, Int, Int32, searchsorted, zeros, array
1222    precision = Float
1223
1224
1225    #Get NetCDF data
1226    if verbose: print 'Reading files %s_*.nc' %basename_in
1227    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
1228    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
1229    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
1230
1231    if basename_out is None:
1232        swwname = basename_in + '.sww'
1233    else:
1234        swwname = basename_out + '.sww'       
1235
1236    times = file_h.variables['TIME']
1237    latitudes = file_h.variables['LAT']
1238    longitudes = file_h.variables['LON']
1239
1240    if mint == None:
1241        jmin = 0
1242    else:   
1243        jmin = searchsorted(times, mint)
1244   
1245    if maxt == None:
1246        jmax=len(times)
1247    else:   
1248        jmax = searchsorted(times, maxt)
1249       
1250    if minlat == None:
1251        kmin=0
1252    else:
1253        kmin = searchsorted(latitudes, minlat)
1254
1255    if maxlat == None:
1256        kmax = len(latitudes)
1257    else:
1258        kmax = searchsorted(latitudes, maxlat)
1259
1260    if minlon == None:
1261        lmin=0
1262    else:   
1263        lmin = searchsorted(longitudes, minlon)
1264   
1265    if maxlon == None:
1266        lmax = len(longitudes)
1267    else:
1268        lmax = searchsorted(longitudes, maxlon)           
1269
1270    times = times[jmin:jmax]       
1271    latitudes = latitudes[kmin:kmax]
1272    longitudes = longitudes[lmin:lmax]
1273
1274
1275    if verbose: print 'cropping'
1276    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1277    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
1278    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
1279
1280    number_of_times = times.shape[0]
1281    number_of_latitudes = latitudes.shape[0]
1282    number_of_longitudes = longitudes.shape[0]
1283
1284    assert amplitudes.shape[0] == number_of_times
1285    assert amplitudes.shape[1] == number_of_latitudes
1286    assert amplitudes.shape[2] == number_of_longitudes         
1287
1288    if verbose:
1289        print '------------------------------------------------'
1290        print 'Statistics:'
1291        print '  Extent (lat/lon):'
1292        print '    lat in [%f, %f], len(lat) == %d'\
1293              %(min(latitudes.flat), max(latitudes.flat),
1294                len(latitudes.flat))
1295        print '    lon in [%f, %f], len(lon) == %d'\
1296              %(min(longitudes.flat), max(longitudes.flat),
1297                len(longitudes.flat))
1298        print '    t in [%f, %f], len(t) == %d'\
1299              %(min(times.flat), max(times.flat), len(times.flat))
1300       
1301        q = amplitudes.flat       
1302        name = 'Amplitudes (ha) [cm]'
1303        print '  %s in [%f, %f]' %(name, min(q), max(q))
1304       
1305        q = uspeed.flat
1306        name = 'Speeds (ua) [cm/s]'       
1307        print '  %s in [%f, %f]' %(name, min(q), max(q))
1308       
1309        q = vspeed.flat
1310        name = 'Speeds (va) [cm/s]'               
1311        print '  %s in [%f, %f]' %(name, min(q), max(q))               
1312
1313
1314    #print number_of_latitudes, number_of_longitudes
1315    number_of_points = number_of_latitudes*number_of_longitudes
1316    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
1317   
1318
1319    file_h.close()
1320    file_u.close()
1321    file_v.close()   
1322
1323
1324    if verbose: print 'Store to SWW file %s' %swwname
1325    # NetCDF file definition
1326    outfile = NetCDFFile(swwname, 'w')
1327       
1328    #Create new file
1329    outfile.institution = 'Geoscience Australia'
1330    outfile.description = 'Converted from Ferret files: %s, %s, %s'\
1331                          %(basename_in + '_ha.nc',
1332                            basename_in + '_ua.nc',
1333                            basename_in + '_va.nc')
1334
1335
1336    #For sww compatibility
1337    outfile.smoothing = 'Yes' 
1338    outfile.order = 1
1339           
1340    #Start time in seconds since the epoch (midnight 1/1/1970)
1341    outfile.starttime = times[0]
1342   
1343    # dimension definitions
1344    outfile.createDimension('number_of_volumes', number_of_volumes)
1345
1346    outfile.createDimension('number_of_vertices', 3)
1347    outfile.createDimension('number_of_points', number_of_points)
1348                           
1349               
1350    #outfile.createDimension('number_of_timesteps', len(times))
1351    outfile.createDimension('number_of_timesteps', len(times))
1352
1353    # variable definitions
1354    outfile.createVariable('x', precision, ('number_of_points',))
1355    outfile.createVariable('y', precision, ('number_of_points',))
1356    outfile.createVariable('elevation', precision, ('number_of_points',))
1357           
1358    #FIXME: Backwards compatibility
1359    outfile.createVariable('z', precision, ('number_of_points',))
1360    #################################
1361       
1362    outfile.createVariable('volumes', Int, ('number_of_volumes',
1363                                            'number_of_vertices'))
1364   
1365    outfile.createVariable('time', precision,
1366                           ('number_of_timesteps',))
1367       
1368    outfile.createVariable('stage', precision,
1369                           ('number_of_timesteps',
1370                            'number_of_points'))
1371
1372    outfile.createVariable('xmomentum', precision,
1373                           ('number_of_timesteps',
1374                            'number_of_points'))
1375   
1376    outfile.createVariable('ymomentum', precision,
1377                           ('number_of_timesteps',
1378                            'number_of_points'))
1379
1380
1381    #Store
1382    from coordinate_transforms.redfearn import redfearn
1383    x = zeros(number_of_points, Float)  #Easting
1384    y = zeros(number_of_points, Float)  #Northing
1385
1386
1387    #Check zone boundaries
1388    refzone, _, _ = redfearn(latitudes[0],longitudes[0])   
1389
1390    vertices = {}
1391    i = 0   
1392    for k, lat in enumerate(latitudes):       #Y direction
1393        for l, lon in enumerate(longitudes):  #X direction
1394
1395            vertices[l,k] = i
1396
1397            zone, easting, northing = redfearn(lat,lon)               
1398
1399            msg = 'Zone boundary crossed at longitude =', lon
1400            assert zone == refzone, msg
1401            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
1402            x[i] = easting
1403            y[i] = northing
1404            i += 1
1405
1406
1407    #Construct 2 triangles per 'rectangular' element
1408    volumes = []
1409    for l in range(number_of_longitudes-1):    #X direction
1410        for k in range(number_of_latitudes-1): #Y direction
1411            v1 = vertices[l,k+1]
1412            v2 = vertices[l,k]           
1413            v3 = vertices[l+1,k+1]           
1414            v4 = vertices[l+1,k]           
1415
1416            volumes.append([v1,v2,v3]) #Upper element
1417            volumes.append([v4,v3,v2]) #Lower element           
1418
1419    volumes = array(volumes)       
1420
1421    if origin == None:
1422        zone = refzone       
1423        xllcorner = min(x)
1424        yllcorner = min(y)
1425    else:
1426        zone = origin[0]
1427        xllcorner = origin[1]
1428        yllcorner = origin[2]               
1429
1430   
1431    outfile.xllcorner = xllcorner
1432    outfile.yllcorner = yllcorner
1433    outfile.zone = zone
1434
1435
1436    if elevation is not None:
1437        z = elevation
1438    else:
1439        pass
1440        #FIXME: z should be obtained from MOST and passed in here
1441       
1442    outfile.variables['x'][:] = x - xllcorner
1443    outfile.variables['y'][:] = y - yllcorner
1444    outfile.variables['z'][:] = z
1445    outfile.variables['elevation'][:] = z  #FIXME HACK
1446    outfile.variables['time'][:] = times
1447    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
1448   
1449
1450
1451    #Time stepping
1452    stage = outfile.variables['stage']
1453    xmomentum = outfile.variables['xmomentum']
1454    ymomentum = outfile.variables['ymomentum']       
1455
1456
1457   
1458    for j in range(len(times)):
1459        i = 0
1460        for k in range(number_of_latitudes):      #Y direction
1461            for l in range(number_of_longitudes): #X direction
1462                w = zscale*amplitudes[j,k,l]/100 + mean_stage
1463                stage[j,i] = w
1464                h = w - z
1465                xmomentum[j,i] = uspeed[j,k,l]/100*h
1466                ymomentum[j,i] = vspeed[j,k,l]/100*h
1467                i += 1
1468
1469
1470    if verbose:
1471        x = outfile.variables['x'][:]
1472        y = outfile.variables['y'][:]       
1473       
1474        print '------------------------------------------------'
1475        print 'Statistics of output file:'
1476        print '  Name: %s' %swwname
1477        print '  Reference:'
1478        print '    Lower left corner: [%f, %f]'\
1479              %(xllcorner, yllcorner)
1480        print '  Extent:'
1481        print '    x [m] in [%f, %f], len(x) == %d'\
1482              %(min(x.flat), max(x.flat), len(x.flat))
1483        print '    y [m] in [%f, %f], len(y) == %d'\
1484              %(min(y.flat), max(y.flat), len(y.flat))
1485        print '    t [s] in [%f, %f], len(t) == %d'\
1486              %(min(times), max(times), len(times))
1487        print '  Quantities [SI units]:'
1488        for name in ['stage', 'xmomentum', 'ymomentum']:
1489            q = outfile.variables[name][:].flat
1490            print '    %s in [%f, %f]' %(name, min(q), max(q))
1491           
1492       
1493
1494       
1495    outfile.close()               
1496   
1497
1498
1499def extent_sww(file_name):
1500    """
1501    Read in an sww file.
1502
1503    Input;
1504    file_name - the sww file
1505
1506    Output;
1507    z - Vector of bed elevation
1508    volumes - Array.  Each row has 3 values, representing
1509    the vertices that define the volume
1510    time - Vector of the times where there is stage information
1511    stage - array with respect to time and vertices (x,y)
1512    """
1513
1514
1515    from Scientific.IO.NetCDF import NetCDFFile
1516
1517    #Check contents
1518    #Get NetCDF
1519    fid = NetCDFFile(file_name, 'r')
1520
1521    # Get the variables
1522    x = fid.variables['x'][:]
1523    y = fid.variables['y'][:]
1524    stage = fid.variables['stage'][:]
1525    #print "stage",stage
1526    #print "stage.shap",stage.shape
1527    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
1528    #print "min(stage)",min(stage)
1529
1530    fid.close()
1531
1532    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
1533
1534
1535def sww2domain(filename,t=None,fail_if_NaN=True,NaN_filler=0):
1536    """Read sww Net CDF file containing Shallow Water Wave simulation
1537
1538    Quantities stage, elevation, xmomentum and ymomentum.
1539
1540    The momentum is not always stored.
1541
1542    """
1543    NaN=9.969209968386869e+036
1544    from Scientific.IO.NetCDF import NetCDFFile
1545    from domain import Domain
1546    from Numeric import asarray, transpose
1547    #print 'Reading from ', filename
1548    fid = NetCDFFile(filename, 'r')    #Open existing file for read
1549    time = fid.variables['time']       #Timesteps
1550    if t is None:
1551        t = time[-1]
1552    time_interp = get_time_interp(time,t)
1553################################
1554########################################   
1555    # Get the variables as Numeric arrays
1556    x = fid.variables['x'][:]             #x-coordinates of vertices
1557    y = fid.variables['y'][:]             #y-coordinates of vertices
1558    elevation = fid.variables['elevation']     #Elevation
1559    stage = fid.variables['stage']     #Water level
1560    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
1561    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction   
1562#################################
1563    xllcorner = fid.xllcorner[0]
1564    yllcorner = fid.yllcorner[0]
1565    starttime = fid.starttime[0]
1566    zone = fid.zone
1567    volumes = fid.variables['volumes'][:] #Connectivity
1568    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
1569#
1570    conserved_quantities = []
1571    interpolated_quantities = {}
1572    other_quantities = []
1573#
1574    #print '    interpolating quantities'
1575    for quantity in fid.variables.keys():
1576        dimensions = fid.variables[quantity].dimensions
1577        if 'number_of_timesteps' in dimensions:
1578            conserved_quantities.append(quantity)
1579            interpolated_quantities[quantity]=\
1580                  interpolated_quantity(fid.variables[quantity][:],time_interp)
1581        else: other_quantities.append(quantity)
1582#
1583    other_quantities.remove('x')
1584    other_quantities.remove('y')
1585    other_quantities.remove('z')
1586    other_quantities.remove('volumes')
1587#
1588    conserved_quantities.remove('time')
1589#
1590    #print other_quantities
1591    #print conserved_quantities
1592    #print '    building domain'
1593    domain = Domain(coordinates, volumes,\
1594                    conserved_quantities = conserved_quantities,\
1595                    other_quantities = other_quantities,zone=zone,\
1596                    xllcorner=xllcorner, yllcorner=yllcorner)
1597    domain.starttime=starttime
1598    domain.time=t
1599    for quantity in other_quantities:
1600        X = fid.variables[quantity][:]
1601        #print quantity
1602        #print 'max(X)'
1603        #print max(X)
1604        #print 'max(X)==NaN'
1605        #print max(X)==NaN
1606        if (max(X)==NaN) or (min(X)==NaN):
1607            if fail_if_NaN:
1608                msg = 'quantity %s contains no_data entry'%quantity
1609                raise msg
1610            else:
1611                data = (X<>NaN)
1612                X = (X*data)+(data==0)*NaN_filler
1613        domain.set_quantity(quantity,X)
1614#
1615    for quantity in conserved_quantities:
1616        X = interpolated_quantities[quantity]
1617        #print quantity
1618        #print 'max(X)'
1619        #print max(X)
1620        #print 'max(X)==NaN'
1621        #print max(X)==NaN
1622        if (max(X)==NaN) or (min(X)==NaN):
1623            if fail_if_NaN:
1624                msg = 'quantity %s contains no_data entry'%quantity
1625                raise msg
1626            else:
1627                data = (X<>NaN)
1628                X = (X*data)+(data==0)*NaN_filler
1629        domain.set_quantity(quantity,X)
1630    return domain
1631
1632def interpolated_quantity(saved_quantity,time_interp):
1633       
1634    #given an index and ratio, interpolate quantity with respect to time.
1635    index,ratio = time_interp
1636    Q = saved_quantity 
1637    if ratio > 0:
1638        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
1639    else:
1640         q = Q[index]
1641    #Return vector of interpolated values
1642    return q
1643
1644def get_time_interp(time,t=None):
1645    #Finds the ratio and index for time interpolation.
1646    #It is borrowed from previous pyvolution code.
1647    if t is None:
1648        t=time[-1]
1649        index = -1
1650        ratio = 0.
1651    else:
1652        T = time
1653        tau = t
1654        index=0
1655        msg = 'Time interval derived from file %s [%s:%s]'\
1656            %('FIXMEfilename', T[0], T[-1])
1657        msg += ' does not match model time: %s' %tau
1658        if tau < time[0]: raise msg
1659        if tau > time[-1]: raise msg
1660        while tau > time[index]: index += 1
1661        while tau < time[index]: index -= 1
1662        if tau == time[index]:
1663            #Protect against case where tau == time[-1] (last time)
1664            # - also works in general when tau == time[i]
1665            ratio = 0
1666        else:
1667            #t is now between index and index+1           
1668            ratio = (tau - time[index])/(time[index+1] - time[index])
1669    return (index,ratio)
1670
1671
1672
1673#OBSOLETE STUFF
1674#Native checkpoint format.
1675#Information needed to recreate a state is preserved
1676#FIXME: Rethink and maybe use netcdf format
1677def cpt_variable_writer(filename, t, v0, v1, v2):
1678    """Store all conserved quantities to file
1679    """
1680
1681    M, N = v0.shape
1682
1683    FN = create_filename(filename, 'cpt', M, t)
1684    #print 'Writing to %s' %FN
1685
1686    fid = open(FN, 'w')
1687    for i in range(M):
1688        for j in range(N):
1689            fid.write('%.16e ' %v0[i,j])
1690        for j in range(N):
1691            fid.write('%.16e ' %v1[i,j])
1692        for j in range(N):
1693            fid.write('%.16e ' %v2[i,j])
1694
1695        fid.write('\n')
1696    fid.close()
1697
1698
1699def cpt_variable_reader(filename, t, v0, v1, v2):
1700    """Store all conserved quantities to file
1701    """
1702
1703    M, N = v0.shape
1704
1705    FN = create_filename(filename, 'cpt', M, t)
1706    #print 'Reading from %s' %FN
1707
1708    fid = open(FN)
1709
1710
1711    for i in range(M):
1712        values = fid.readline().split() #Get one line
1713
1714        for j in range(N):
1715            v0[i,j] = float(values[j])
1716            v1[i,j] = float(values[3+j])
1717            v2[i,j] = float(values[6+j])
1718
1719    fid.close()
1720
1721def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1722    """Writes x,y,z,z,z coordinates of triangles constituting the bed
1723    elevation.
1724    Not in use pt
1725    """
1726
1727    M, N = v0.shape
1728
1729    print X0
1730    import sys; sys.exit()
1731    FN = create_filename(filename, 'cpt', M)
1732    print 'Writing to %s' %FN
1733
1734    fid = open(FN, 'w')
1735    for i in range(M):
1736        for j in range(2):
1737            fid.write('%.16e ' %X0[i,j])   #x, y
1738        for j in range(N):
1739            fid.write('%.16e ' %v0[i,j])       #z,z,z,
1740
1741        for j in range(2):
1742            fid.write('%.16e ' %X1[i,j])   #x, y
1743        for j in range(N):
1744            fid.write('%.16e ' %v1[i,j])
1745
1746        for j in range(2):
1747            fid.write('%.16e ' %X2[i,j])   #x, y
1748        for j in range(N):
1749            fid.write('%.16e ' %v2[i,j])
1750
1751        fid.write('\n')
1752    fid.close()
1753
1754
1755
1756#Function for storing out to e.g. visualisation
1757#FIXME: Do we want this?
1758#FIXME: Not done yet for this version
1759def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1760    """Writes x,y,z coordinates of triangles constituting the bed elevation.
1761    """
1762
1763    M, N = v0.shape
1764
1765    FN = create_filename(filename, 'dat', M)
1766    #print 'Writing to %s' %FN
1767
1768    fid = open(FN, 'w')
1769    for i in range(M):
1770        for j in range(2):
1771            fid.write('%f ' %X0[i,j])   #x, y
1772        fid.write('%f ' %v0[i,0])       #z
1773
1774        for j in range(2):
1775            fid.write('%f ' %X1[i,j])   #x, y
1776        fid.write('%f ' %v1[i,0])       #z
1777
1778        for j in range(2):
1779            fid.write('%f ' %X2[i,j])   #x, y
1780        fid.write('%f ' %v2[i,0])       #z
1781
1782        fid.write('\n')
1783    fid.close()
1784
1785
1786
1787def dat_variable_writer(filename, t, v0, v1, v2):
1788    """Store water height to file
1789    """
1790
1791    M, N = v0.shape
1792
1793    FN = create_filename(filename, 'dat', M, t)
1794    #print 'Writing to %s' %FN
1795
1796    fid = open(FN, 'w')
1797    for i in range(M):
1798        fid.write('%.4f ' %v0[i,0])
1799        fid.write('%.4f ' %v1[i,0])
1800        fid.write('%.4f ' %v2[i,0])
1801
1802        fid.write('\n')
1803    fid.close()
1804
1805
1806def read_sww(filename):
1807    """Read sww Net CDF file containing Shallow Water Wave simulation
1808
1809    The integer array volumes is of shape Nx3 where N is the number of
1810    triangles in the mesh.
1811
1812    Each entry in volumes is an index into the x,y arrays (the location).
1813
1814    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
1815    number_of_timesteps, number_of_points.
1816
1817    The momentum is not always stored.
1818   
1819    """
1820    from Scientific.IO.NetCDF import NetCDFFile
1821    print 'Reading from ', filename
1822    fid = NetCDFFile(filename, 'r')    #Open existing file for read
1823
1824    # Get the variables as Numeric arrays
1825    x = fid.variables['x']             #x-coordinates of vertices
1826    y = fid.variables['y']             #y-coordinates of vertices
1827    z = fid.variables['elevation']     #Elevation
1828    time = fid.variables['time']       #Timesteps
1829    stage = fid.variables['stage']     #Water level
1830    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
1831    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction   
1832
1833    volumes = fid.variables['volumes'] #Connectivity
Note: See TracBrowser for help on using the repository browser.