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

Last change on this file since 1103 was 1103, checked in by ole, 19 years ago

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