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

Last change on this file since 1114 was 1114, checked in by prow, 20 years ago
File size: 56.4 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               fail_on_NaN = True,
1187               NaN_filler = 0,
1188               elevation = -100
1189               ): #FIXME: Bathymetry should be obtained
1190                                  #from MOST somehow.
1191                                  #Alternatively from elsewhere
1192                                  #or, as a last resort,
1193                                  #specified here.
1194                                  #The value of -100 will work
1195                                  #for the Wollongong tsunami
1196                                  #scenario but is very hacky
1197    """Convert 'Ferret' NetCDF format for wave propagation to
1198    sww format native to pyvolution.
1199
1200    Specify only basename_in and read files of the form
1201    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1202    relative height, x-velocity and y-velocity, respectively.
1203
1204    Also convert latitude and longitude to UTM. All coordinates are
1205    assumed to be given in the GDA94 datum.
1206
1207    min's and max's: If omitted - full extend is used.
1208    To include a value min may equal it, while max must exceed it.
1209    Lat and lon are assuemd to be in decimal degrees
1210
1211    origin is a 3-tuple with geo referenced
1212    UTM coordinates (zone, easting, northing)
1213
1214    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
1215    which means that longitude is the fastest
1216    varying dimension (row major order, so to speak)
1217
1218    ferret2sww uses grid points as vertices in a triangular grid
1219    counting vertices from lower left corner upwards, then right
1220    """ 
1221
1222    import os
1223    from Scientific.IO.NetCDF import NetCDFFile
1224    from Numeric import Float, Int, Int32, searchsorted, zeros, array
1225    precision = Float
1226
1227
1228    #Get NetCDF data
1229    if verbose: print 'Reading files %s_*.nc' %basename_in
1230    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
1231    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
1232    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
1233    file_e = NetCDFFile(basename_in + '_va.nc', 'r') #Elevation (z) (m)   
1234
1235    if basename_out is None:
1236        swwname = basename_in + '.sww'
1237    else:
1238        swwname = basename_out + '.sww'       
1239
1240    times = file_h.variables['TIME']
1241    latitudes = file_h.variables['LAT']
1242    longitudes = file_h.variables['LON']
1243
1244    if mint == None:
1245        jmin = 0
1246    else:   
1247        jmin = searchsorted(times, mint)
1248   
1249    if maxt == None:
1250        jmax=len(times)
1251    else:   
1252        jmax = searchsorted(times, maxt)
1253       
1254    if minlat == None:
1255        kmin=0
1256    else:
1257        kmin = searchsorted(latitudes, minlat)
1258
1259    if maxlat == None:
1260        kmax = len(latitudes)
1261    else:
1262        kmax = searchsorted(latitudes, maxlat)
1263
1264    if minlon == None:
1265        lmin=0
1266    else:   
1267        lmin = searchsorted(longitudes, minlon)
1268   
1269    if maxlon == None:
1270        lmax = len(longitudes)
1271    else:
1272        lmax = searchsorted(longitudes, maxlon)           
1273
1274    times = times[jmin:jmax]       
1275    latitudes = latitudes[kmin:kmax]
1276    longitudes = longitudes[lmin:lmax]
1277
1278
1279    if verbose: print 'cropping'
1280    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1281    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
1282    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
1283    elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
1284    #Get missing values
1285    nan_ha = file_h.variables['HA'].missing_value[0]
1286    nan_ua = file_u.variables['UA'].missing_value[0]
1287    nan_va = file_v.variables['VA'].missing_value[0]
1288    nan_e  = file_e.variables['ELEVATION'].missing_value[0]
1289
1290   
1291    #Cleanup
1292    from Numeric import sometrue
1293   
1294    missing = (amplitudes == nan_ha)
1295    if sometrue (missing):
1296        if fail_on_NaN:
1297            msg = 'NetCDFFile %s contains missing values'\
1298                  %(basename_in+'_ha.nc')
1299            raise msg
1300        else:
1301            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
1302
1303    missing = (uspeed == nan_ua)
1304    if sometrue (missing):
1305        if fail_on_NaN:
1306            msg = 'NetCDFFile %s contains missing values'\
1307                  %(basename_in+'_ua.nc')
1308            raise msg
1309        else:
1310            uspeed = uspeed*(missing==0) + missing*NaN_filler
1311
1312    missing = (vspeed == nan_va)
1313    if sometrue (missing):
1314        if fail_on_NaN:
1315            msg = 'NetCDFFile %s contains missing values'\
1316                  %(basename_in+'_va.nc')
1317            raise msg
1318        else:
1319            vspeed = vspeed*(missing==0) + missing*NaN_filler
1320
1321
1322    missing = (elevations == nan_e)
1323    if sometrue (missing):
1324        if fail_on_NaN:
1325            msg = 'NetCDFFile %s contains missing values'\
1326                  %(basename_in+'_e.nc')
1327            raise msg
1328        else:
1329            elevations = elevations*(missing==0) + missing*NaN_filler
1330
1331    #######
1332
1333               
1334
1335    number_of_times = times.shape[0]
1336    number_of_latitudes = latitudes.shape[0]
1337    number_of_longitudes = longitudes.shape[0]
1338
1339    assert amplitudes.shape[0] == number_of_times
1340    assert amplitudes.shape[1] == number_of_latitudes
1341    assert amplitudes.shape[2] == number_of_longitudes         
1342
1343    if verbose:
1344        print '------------------------------------------------'
1345        print 'Statistics:'
1346        print '  Extent (lat/lon):'
1347        print '    lat in [%f, %f], len(lat) == %d'\
1348              %(min(latitudes.flat), max(latitudes.flat),
1349                len(latitudes.flat))
1350        print '    lon in [%f, %f], len(lon) == %d'\
1351              %(min(longitudes.flat), max(longitudes.flat),
1352                len(longitudes.flat))
1353        print '    t in [%f, %f], len(t) == %d'\
1354              %(min(times.flat), max(times.flat), len(times.flat))
1355       
1356        q = amplitudes.flat       
1357        name = 'Amplitudes (ha) [cm]'
1358        print '  %s in [%f, %f]' %(name, min(q), max(q))
1359       
1360        q = uspeed.flat
1361        name = 'Speeds (ua) [cm/s]'       
1362        print '  %s in [%f, %f]' %(name, min(q), max(q))
1363       
1364        q = vspeed.flat
1365        name = 'Speeds (va) [cm/s]'               
1366        print '  %s in [%f, %f]' %(name, min(q), max(q))               
1367
1368        q = elevations.flat
1369        name = 'Elevations (e) [m]'               
1370        print '  %s in [%f, %f]' %(name, min(q), max(q))               
1371
1372
1373    #print number_of_latitudes, number_of_longitudes
1374    number_of_points = number_of_latitudes*number_of_longitudes
1375    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
1376   
1377
1378    file_h.close()
1379    file_u.close()
1380    file_v.close()   
1381    file_e.close()   
1382
1383
1384    # NetCDF file definition
1385    outfile = NetCDFFile(swwname, 'w')
1386       
1387    #Create new file
1388    outfile.institution = 'Geoscience Australia'
1389    outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
1390                          %(basename_in + '_ha.nc',
1391                            basename_in + '_ua.nc',
1392                            basename_in + '_va.nc'
1393                            basename_in + '_e.nc')
1394
1395
1396    #For sww compatibility
1397    outfile.smoothing = 'Yes' 
1398    outfile.order = 1
1399           
1400    #Start time in seconds since the epoch (midnight 1/1/1970)
1401    outfile.starttime = starttime = times[0]
1402    times = times - starttime  #Store relative times
1403   
1404    # dimension definitions
1405    outfile.createDimension('number_of_volumes', number_of_volumes)
1406
1407    outfile.createDimension('number_of_vertices', 3)
1408    outfile.createDimension('number_of_points', number_of_points)
1409                           
1410               
1411    #outfile.createDimension('number_of_timesteps', len(times))
1412    outfile.createDimension('number_of_timesteps', len(times))
1413
1414    # variable definitions
1415    outfile.createVariable('x', precision, ('number_of_points',))
1416    outfile.createVariable('y', precision, ('number_of_points',))
1417    outfile.createVariable('elevation', precision, ('number_of_points',))
1418           
1419    #FIXME: Backwards compatibility
1420    outfile.createVariable('z', precision, ('number_of_points',))
1421    #################################
1422       
1423    outfile.createVariable('volumes', Int, ('number_of_volumes',
1424                                            'number_of_vertices'))
1425   
1426    outfile.createVariable('time', precision,
1427                           ('number_of_timesteps',))
1428       
1429    outfile.createVariable('stage', precision,
1430                           ('number_of_timesteps',
1431                            'number_of_points'))
1432
1433    outfile.createVariable('xmomentum', precision,
1434                           ('number_of_timesteps',
1435                            'number_of_points'))
1436   
1437    outfile.createVariable('ymomentum', precision,
1438                           ('number_of_timesteps',
1439                            'number_of_points'))
1440
1441
1442    #Store
1443    from coordinate_transforms.redfearn import redfearn
1444    x = zeros(number_of_points, Float)  #Easting
1445    y = zeros(number_of_points, Float)  #Northing
1446
1447
1448    if verbose: print 'Making triangular grid'
1449    #Check zone boundaries
1450    refzone, _, _ = redfearn(latitudes[0],longitudes[0])   
1451
1452    vertices = {}
1453    i = 0   
1454    for k, lat in enumerate(latitudes):       #Y direction
1455        for l, lon in enumerate(longitudes):  #X direction
1456
1457            vertices[l,k] = i
1458
1459            zone, easting, northing = redfearn(lat,lon)               
1460
1461            msg = 'Zone boundary crossed at longitude =', lon
1462            assert zone == refzone, msg
1463            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
1464            x[i] = easting
1465            y[i] = northing
1466            i += 1
1467
1468
1469    #Construct 2 triangles per 'rectangular' element
1470    volumes = []
1471    for l in range(number_of_longitudes-1):    #X direction
1472        for k in range(number_of_latitudes-1): #Y direction
1473            v1 = vertices[l,k+1]
1474            v2 = vertices[l,k]           
1475            v3 = vertices[l+1,k+1]           
1476            v4 = vertices[l+1,k]           
1477
1478            volumes.append([v1,v2,v3]) #Upper element
1479            volumes.append([v4,v3,v2]) #Lower element           
1480
1481    volumes = array(volumes)       
1482
1483    if origin == None:
1484        zone = refzone       
1485        xllcorner = min(x)
1486        yllcorner = min(y)
1487    else:
1488        zone = origin[0]
1489        xllcorner = origin[1]
1490        yllcorner = origin[2]               
1491
1492   
1493    outfile.xllcorner = xllcorner
1494    outfile.yllcorner = yllcorner
1495    outfile.zone = zone
1496
1497
1498    if elevation is not None:
1499        z = elevation
1500    else:
1501        if inverted_bathymetry:
1502            z = -1*elevations
1503        else:
1504            z = elevations
1505        #FIXME: z should be obtained from MOST and passed in here
1506       
1507    outfile.variables['x'][:] = x - xllcorner
1508    outfile.variables['y'][:] = y - yllcorner
1509    outfile.variables['z'][:] = z
1510    outfile.variables['elevation'][:] = z  #FIXME HACK
1511    outfile.variables['time'][:] = times   #Store time relative
1512    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
1513   
1514
1515
1516    #Time stepping
1517    stage = outfile.variables['stage']
1518    xmomentum = outfile.variables['xmomentum']
1519    ymomentum = outfile.variables['ymomentum'] 
1520    z = outfile.variables['elevation']         
1521
1522
1523    if verbose: print 'Converting quantities'
1524    n = len(times)
1525    for j in range(n):
1526        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
1527        i = 0
1528        for k in range(number_of_latitudes):      #Y direction
1529            for l in range(number_of_longitudes): #X direction
1530                w = zscale*amplitudes[j,k,l]/100 + mean_stage
1531                stage[j,i] = w
1532                h = w - z[j,i]
1533                xmomentum[j,i] = uspeed[j,k,l]/100*h
1534                ymomentum[j,i] = vspeed[j,k,l]/100*h
1535                i += 1
1536
1537
1538    if verbose:
1539        x = outfile.variables['x'][:]
1540        y = outfile.variables['y'][:]       
1541        print '------------------------------------------------'
1542        print 'Statistics of output file:'
1543        print '  Name: %s' %swwname
1544        print '  Reference:'
1545        print '    Lower left corner: [%f, %f]'\
1546              %(xllcorner, yllcorner)
1547        print '    Start time: %f' %starttime
1548        print '  Extent:'
1549        print '    x [m] in [%f, %f], len(x) == %d'\
1550              %(min(x.flat), max(x.flat), len(x.flat))
1551        print '    y [m] in [%f, %f], len(y) == %d'\
1552              %(min(y.flat), max(y.flat), len(y.flat))
1553        print '    t [s] in [%f, %f], len(t) == %d'\
1554              %(min(times), max(times), len(times))
1555        print '  Quantities [SI units]:'
1556        for name in ['stage', 'xmomentum', 'ymomentum']:
1557            q = outfile.variables[name][:].flat
1558            print '    %s in [%f, %f]' %(name, min(q), max(q))
1559           
1560       
1561
1562       
1563    outfile.close()               
1564   
1565
1566
1567def extent_sww(file_name):
1568    """
1569    Read in an sww file.
1570
1571    Input;
1572    file_name - the sww file
1573
1574    Output;
1575    z - Vector of bed elevation
1576    volumes - Array.  Each row has 3 values, representing
1577    the vertices that define the volume
1578    time - Vector of the times where there is stage information
1579    stage - array with respect to time and vertices (x,y)
1580    """
1581
1582
1583    from Scientific.IO.NetCDF import NetCDFFile
1584
1585    #Check contents
1586    #Get NetCDF
1587    fid = NetCDFFile(file_name, 'r')
1588
1589    # Get the variables
1590    x = fid.variables['x'][:]
1591    y = fid.variables['y'][:]
1592    stage = fid.variables['stage'][:]
1593    #print "stage",stage
1594    #print "stage.shap",stage.shape
1595    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
1596    #print "min(stage)",min(stage)
1597
1598    fid.close()
1599
1600    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
1601
1602
1603def sww2domain(filename,t=None,fail_if_NaN=True,NaN_filler=0):
1604    """Read sww Net CDF file containing Shallow Water Wave simulation
1605
1606    Quantities stage, elevation, xmomentum and ymomentum.
1607
1608    The momentum is not always stored.
1609
1610    """
1611    NaN=9.969209968386869e+036
1612    from Scientific.IO.NetCDF import NetCDFFile
1613    from domain import Domain
1614    from Numeric import asarray, transpose
1615    #print 'Reading from ', filename
1616    fid = NetCDFFile(filename, 'r')    #Open existing file for read
1617    time = fid.variables['time']       #Timesteps
1618    if t is None:
1619        t = time[-1]
1620    time_interp = get_time_interp(time,t)
1621################################
1622########################################   
1623    # Get the variables as Numeric arrays
1624    x = fid.variables['x'][:]             #x-coordinates of vertices
1625    y = fid.variables['y'][:]             #y-coordinates of vertices
1626    elevation = fid.variables['elevation']     #Elevation
1627    stage = fid.variables['stage']     #Water level
1628    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
1629    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction   
1630#################################
1631    xllcorner = fid.xllcorner[0]
1632    yllcorner = fid.yllcorner[0]
1633    starttime = fid.starttime[0]
1634    zone = fid.zone
1635    volumes = fid.variables['volumes'][:] #Connectivity
1636    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
1637#
1638    conserved_quantities = []
1639    interpolated_quantities = {}
1640    other_quantities = []
1641#
1642    #print '    interpolating quantities'
1643    for quantity in fid.variables.keys():
1644        dimensions = fid.variables[quantity].dimensions
1645        if 'number_of_timesteps' in dimensions:
1646            conserved_quantities.append(quantity)
1647            interpolated_quantities[quantity]=\
1648                  interpolated_quantity(fid.variables[quantity][:],time_interp)
1649        else: other_quantities.append(quantity)
1650#
1651    other_quantities.remove('x')
1652    other_quantities.remove('y')
1653    other_quantities.remove('z')
1654    other_quantities.remove('volumes')
1655#
1656    conserved_quantities.remove('time')
1657#
1658    #print other_quantities
1659    #print conserved_quantities
1660    #print '    building domain'
1661    domain = Domain(coordinates, volumes,\
1662                    conserved_quantities = conserved_quantities,\
1663                    other_quantities = other_quantities,zone=zone,\
1664                    xllcorner=xllcorner, yllcorner=yllcorner)
1665    domain.starttime=starttime
1666    domain.time=t
1667    for quantity in other_quantities:
1668        X = fid.variables[quantity][:]
1669        #print quantity
1670        #print 'max(X)'
1671        #print max(X)
1672        #print 'max(X)==NaN'
1673        #print max(X)==NaN
1674        if (max(X)==NaN) or (min(X)==NaN):
1675            if fail_if_NaN:
1676                msg = 'quantity %s contains no_data entry'%quantity
1677                raise msg
1678            else:
1679                data = (X<>NaN)
1680                X = (X*data)+(data==0)*NaN_filler
1681        domain.set_quantity(quantity,X)
1682#
1683    for quantity in conserved_quantities:
1684        X = interpolated_quantities[quantity]
1685        #print quantity
1686        #print 'max(X)'
1687        #print max(X)
1688        #print 'max(X)==NaN'
1689        #print max(X)==NaN
1690        if (max(X)==NaN) or (min(X)==NaN):
1691            if fail_if_NaN:
1692                msg = 'quantity %s contains no_data entry'%quantity
1693                raise msg
1694            else:
1695                data = (X<>NaN)
1696                X = (X*data)+(data==0)*NaN_filler
1697        domain.set_quantity(quantity,X)
1698    return domain
1699
1700def interpolated_quantity(saved_quantity,time_interp):
1701       
1702    #given an index and ratio, interpolate quantity with respect to time.
1703    index,ratio = time_interp
1704    Q = saved_quantity 
1705    if ratio > 0:
1706        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
1707    else:
1708         q = Q[index]
1709    #Return vector of interpolated values
1710    return q
1711
1712def get_time_interp(time,t=None):
1713    #Finds the ratio and index for time interpolation.
1714    #It is borrowed from previous pyvolution code.
1715    if t is None:
1716        t=time[-1]
1717        index = -1
1718        ratio = 0.
1719    else:
1720        T = time
1721        tau = t
1722        index=0
1723        msg = 'Time interval derived from file %s [%s:%s]'\
1724            %('FIXMEfilename', T[0], T[-1])
1725        msg += ' does not match model time: %s' %tau
1726        if tau < time[0]: raise msg
1727        if tau > time[-1]: raise msg
1728        while tau > time[index]: index += 1
1729        while tau < time[index]: index -= 1
1730        if tau == time[index]:
1731            #Protect against case where tau == time[-1] (last time)
1732            # - also works in general when tau == time[i]
1733            ratio = 0
1734        else:
1735            #t is now between index and index+1           
1736            ratio = (tau - time[index])/(time[index+1] - time[index])
1737    return (index,ratio)
1738
1739
1740
1741#OBSOLETE STUFF
1742#Native checkpoint format.
1743#Information needed to recreate a state is preserved
1744#FIXME: Rethink and maybe use netcdf format
1745def cpt_variable_writer(filename, t, v0, v1, v2):
1746    """Store all conserved quantities to file
1747    """
1748
1749    M, N = v0.shape
1750
1751    FN = create_filename(filename, 'cpt', M, t)
1752    #print 'Writing to %s' %FN
1753
1754    fid = open(FN, 'w')
1755    for i in range(M):
1756        for j in range(N):
1757            fid.write('%.16e ' %v0[i,j])
1758        for j in range(N):
1759            fid.write('%.16e ' %v1[i,j])
1760        for j in range(N):
1761            fid.write('%.16e ' %v2[i,j])
1762
1763        fid.write('\n')
1764    fid.close()
1765
1766
1767def cpt_variable_reader(filename, t, v0, v1, v2):
1768    """Store all conserved quantities to file
1769    """
1770
1771    M, N = v0.shape
1772
1773    FN = create_filename(filename, 'cpt', M, t)
1774    #print 'Reading from %s' %FN
1775
1776    fid = open(FN)
1777
1778
1779    for i in range(M):
1780        values = fid.readline().split() #Get one line
1781
1782        for j in range(N):
1783            v0[i,j] = float(values[j])
1784            v1[i,j] = float(values[3+j])
1785            v2[i,j] = float(values[6+j])
1786
1787    fid.close()
1788
1789def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1790    """Writes x,y,z,z,z coordinates of triangles constituting the bed
1791    elevation.
1792    Not in use pt
1793    """
1794
1795    M, N = v0.shape
1796
1797    print X0
1798    import sys; sys.exit()
1799    FN = create_filename(filename, 'cpt', M)
1800    print 'Writing to %s' %FN
1801
1802    fid = open(FN, 'w')
1803    for i in range(M):
1804        for j in range(2):
1805            fid.write('%.16e ' %X0[i,j])   #x, y
1806        for j in range(N):
1807            fid.write('%.16e ' %v0[i,j])       #z,z,z,
1808
1809        for j in range(2):
1810            fid.write('%.16e ' %X1[i,j])   #x, y
1811        for j in range(N):
1812            fid.write('%.16e ' %v1[i,j])
1813
1814        for j in range(2):
1815            fid.write('%.16e ' %X2[i,j])   #x, y
1816        for j in range(N):
1817            fid.write('%.16e ' %v2[i,j])
1818
1819        fid.write('\n')
1820    fid.close()
1821
1822
1823
1824#Function for storing out to e.g. visualisation
1825#FIXME: Do we want this?
1826#FIXME: Not done yet for this version
1827def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1828    """Writes x,y,z coordinates of triangles constituting the bed elevation.
1829    """
1830
1831    M, N = v0.shape
1832
1833    FN = create_filename(filename, 'dat', M)
1834    #print 'Writing to %s' %FN
1835
1836    fid = open(FN, 'w')
1837    for i in range(M):
1838        for j in range(2):
1839            fid.write('%f ' %X0[i,j])   #x, y
1840        fid.write('%f ' %v0[i,0])       #z
1841
1842        for j in range(2):
1843            fid.write('%f ' %X1[i,j])   #x, y
1844        fid.write('%f ' %v1[i,0])       #z
1845
1846        for j in range(2):
1847            fid.write('%f ' %X2[i,j])   #x, y
1848        fid.write('%f ' %v2[i,0])       #z
1849
1850        fid.write('\n')
1851    fid.close()
1852
1853
1854
1855def dat_variable_writer(filename, t, v0, v1, v2):
1856    """Store water height to file
1857    """
1858
1859    M, N = v0.shape
1860
1861    FN = create_filename(filename, 'dat', M, t)
1862    #print 'Writing to %s' %FN
1863
1864    fid = open(FN, 'w')
1865    for i in range(M):
1866        fid.write('%.4f ' %v0[i,0])
1867        fid.write('%.4f ' %v1[i,0])
1868        fid.write('%.4f ' %v2[i,0])
1869
1870        fid.write('\n')
1871    fid.close()
1872
1873
1874def read_sww(filename):
1875    """Read sww Net CDF file containing Shallow Water Wave simulation
1876
1877    The integer array volumes is of shape Nx3 where N is the number of
1878    triangles in the mesh.
1879
1880    Each entry in volumes is an index into the x,y arrays (the location).
1881
1882    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
1883    number_of_timesteps, number_of_points.
1884
1885    The momentum is not always stored.
1886   
1887    """
1888    from Scientific.IO.NetCDF import NetCDFFile
1889    print 'Reading from ', filename
1890    fid = NetCDFFile(filename, 'r')    #Open existing file for read
1891
1892    # Get the variables as Numeric arrays
1893    x = fid.variables['x']             #x-coordinates of vertices
1894    y = fid.variables['y']             #y-coordinates of vertices
1895    z = fid.variables['elevation']     #Elevation
1896    time = fid.variables['time']       #Timesteps
1897    stage = fid.variables['stage']     #Water level
1898    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
1899    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction   
1900
1901    volumes = fid.variables['volumes'] #Connectivity
Note: See TracBrowser for help on using the repository browser.