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

Last change on this file since 1142 was 1142, checked in by ole, 20 years ago

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