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

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

Improved optimisation in sww2asc (happy Easter everyone!)
File size: 63.7 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         sometrue
1072   
1073
1074    #FIXME: Should be variable
1075    datum = 'WGS84'
1076    false_easting = 500000
1077    false_northing = 10000000
1078    NODATA_value = -9999
1079    #
1080
1081   
1082    if quantity is None:
1083        quantity = 'elevation'
1084
1085    if reduction is None:
1086        reduction = max
1087
1088    if basename_out is None:
1089        basename_out = basename_in
1090       
1091    swwfile = basename_in + '.sww'
1092    ascfile = basename_out + '.asc'
1093    prjfile = basename_out + '.prj'   
1094
1095
1096    if verbose: print 'Reading from %s' %swwfile
1097    #Read sww file
1098    from Scientific.IO.NetCDF import NetCDFFile
1099    fid = NetCDFFile(swwfile)
1100
1101    #Get extent and reference
1102    x = fid.variables['x'][:]
1103    y = fid.variables['y'][:]
1104    volumes = fid.variables['volumes'][:]
1105
1106    ymin = min(y); ymax = max(y)
1107    xmin = min(x); xmax = max(x)
1108   
1109    number_of_timesteps = fid.dimensions['number_of_timesteps']
1110    number_of_points = fid.dimensions['number_of_points']
1111    if origin is None:
1112        xllcorner = fid.xllcorner[0]
1113        yllcorner = fid.yllcorner[0]
1114        zone = fid.zone[0]
1115    else:
1116        zone = origin[0]
1117        xllcorner = origin[1]
1118        yllcorner = origin[2]
1119       
1120
1121    #Get quantity and reduce if applicable
1122    if verbose: print 'Reading quantity %s' %quantity   
1123    q = fid.variables[quantity][:]
1124
1125   
1126    if len(q.shape) == 2:
1127        if verbose: print 'Reducing quantity %s' %quantity           
1128        q_reduced = zeros( number_of_points, Float ) 
1129
1130        for k in range(number_of_points):
1131            q_reduced[k] = reduction( q[:,k] )
1132
1133        q = q_reduced
1134
1135    #Now q has dimension: number_of_points
1136
1137   
1138    #Write prj file
1139    if verbose: print 'Writing %s' %prjfile
1140    prjid = open(prjfile, 'w')
1141    prjid.write('Projection    %s\n' %'UTM')
1142    prjid.write('Zone          %d\n' %zone)               
1143    prjid.write('Datum         %s\n' %datum)
1144    prjid.write('Zunits        NO\n')
1145    prjid.write('Units         METERS\n')
1146    prjid.write('Spheroid      %s\n' %datum)
1147    prjid.write('Xshift        %d\n' %false_easting)
1148    prjid.write('Yshift        %d\n' %false_northing)
1149    prjid.write('Parameters\n')
1150    prjid.close()
1151
1152    #Create grid and update xll/yll corner and x,y
1153    if verbose: print 'Creating grid'
1154    ncols = int((xmax-xmin)/cellsize)+1
1155    nrows = int((ymax-ymin)/cellsize)+1
1156
1157    newxllcorner = xmin+xllcorner
1158    newyllcorner = ymin+yllcorner
1159   
1160    x = x+xllcorner-newxllcorner
1161    y = y+yllcorner-newyllcorner
1162
1163    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
1164    assert len(vertex_points.shape) == 2
1165   
1166
1167    from Numeric import zeros, Float
1168    grid_points = zeros ( (ncols*nrows, 2), Float )
1169
1170
1171    for i in xrange(nrows):
1172        yg = i*cellsize       
1173        for j in xrange(ncols):
1174            xg = j*cellsize
1175            k = i*ncols + j
1176
1177            #print k, xg, yg
1178            grid_points[k,0] = xg
1179            grid_points[k,1] = yg
1180           
1181
1182    #Interpolate
1183   
1184    from least_squares import Interpolation
1185    from util import inside_polygon
1186
1187    #FIXME: This should be done with precrop = True, otherwise it'll
1188    #take forever. With expand_search  set to False, some grid points might
1189    #miss out....
1190   
1191    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
1192                           precrop = False, expand_search = False,
1193                           verbose = verbose)
1194
1195    #Interpolate using quantity values
1196    if verbose: print 'Interpolating'
1197    grid_values = interp.interpolate(q).flat
1198
1199    #Write
1200    if verbose: print 'Writing %s' %ascfile   
1201    ascid = open(ascfile, 'w')
1202
1203    ascid.write('ncols         %d\n' %ncols)
1204    ascid.write('nrows         %d\n' %nrows)
1205    ascid.write('xllcorner     %d\n' %newxllcorner)
1206    ascid.write('yllcorner     %d\n' %newyllcorner)
1207    ascid.write('cellsize      %f\n' %cellsize)
1208    ascid.write('NODATA_value  %d\n' %NODATA_value)                   
1209
1210
1211    #Get bounding polygon from mesh
1212    P = interp.mesh.get_boundary_polygon()
1213    inside_indices = inside_polygon(grid_points, P)
1214   
1215    for i in range(nrows):
1216        if verbose and i%((nrows+1)/10)==0:
1217            print 'Doing row %d of %d' %(i, nrows)
1218           
1219        for j in range(ncols):
1220            index = (nrows-i-1)*ncols+j
1221           
1222            if sometrue(inside_indices == index):
1223                ascid.write('%f ' %grid_values[index])               
1224            else:
1225                ascid.write('%d ' %NODATA_value)
1226               
1227        ascid.write('\n')
1228       
1229    #Close
1230    ascid.close()
1231    fid.close()
1232   
1233
1234def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
1235                                  verbose=False):
1236    """Read Digitial Elevation model from the following ASCII format (.asc)
1237
1238    Example:
1239
1240    ncols         3121
1241    nrows         1800
1242    xllcorner     722000
1243    yllcorner     5893000
1244    cellsize      25
1245    NODATA_value  -9999
1246    138.3698 137.4194 136.5062 135.5558 ..........
1247
1248    Convert basename_in + '.asc' to NetCDF format (.dem)
1249    mimicking the ASCII format closely.
1250
1251
1252    An accompanying file with same basename_in but extension .prj must exist
1253    and is used to fix the UTM zone, datum, false northings and eastings.
1254
1255    The prj format is assumed to be as
1256
1257    Projection    UTM
1258    Zone          56
1259    Datum         WGS84
1260    Zunits        NO
1261    Units         METERS
1262    Spheroid      WGS84
1263    Xshift        0.0000000000
1264    Yshift        10000000.0000000000
1265    Parameters
1266    """
1267
1268    import os
1269    from Scientific.IO.NetCDF import NetCDFFile
1270    from Numeric import Float, array
1271
1272    #root, ext = os.path.splitext(basename_in)
1273    root = basename_in
1274
1275    ###########################################
1276    # Read Meta data
1277    if verbose: print 'Reading METADATA from %s' %root + '.prj'
1278    metadatafile = open(root + '.prj')
1279    metalines = metadatafile.readlines()
1280    metadatafile.close()
1281
1282    L = metalines[0].strip().split()
1283    assert L[0].strip().lower() == 'projection'
1284    projection = L[1].strip()                   #TEXT
1285
1286    L = metalines[1].strip().split()
1287    assert L[0].strip().lower() == 'zone'
1288    zone = int(L[1].strip())
1289
1290    L = metalines[2].strip().split()
1291    assert L[0].strip().lower() == 'datum'
1292    datum = L[1].strip()                        #TEXT
1293
1294    L = metalines[3].strip().split()
1295    assert L[0].strip().lower() == 'zunits'     #IGNORE
1296    zunits = L[1].strip()                       #TEXT
1297
1298    L = metalines[4].strip().split()
1299    assert L[0].strip().lower() == 'units'
1300    units = L[1].strip()                        #TEXT
1301
1302    L = metalines[5].strip().split()
1303    assert L[0].strip().lower() == 'spheroid'   #IGNORE
1304    spheroid = L[1].strip()                     #TEXT
1305
1306    L = metalines[6].strip().split()
1307    assert L[0].strip().lower() == 'xshift'
1308    false_easting = float(L[1].strip())
1309
1310    L = metalines[7].strip().split()
1311    assert L[0].strip().lower() == 'yshift'
1312    false_northing = float(L[1].strip())
1313
1314    #print false_easting, false_northing, zone, datum
1315
1316
1317    ###########################################
1318    #Read DEM data
1319
1320    datafile = open(basename_in + '.asc')
1321
1322    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
1323    lines = datafile.readlines()
1324    datafile.close()
1325
1326    if verbose: print 'Got', len(lines), ' lines'
1327
1328    ncols = int(lines[0].split()[1].strip())
1329    nrows = int(lines[1].split()[1].strip())
1330    xllcorner = float(lines[2].split()[1].strip())
1331    yllcorner = float(lines[3].split()[1].strip())
1332    cellsize = float(lines[4].split()[1].strip())
1333    NODATA_value = int(lines[5].split()[1].strip())
1334
1335    assert len(lines) == nrows + 6
1336
1337
1338    ##########################################
1339
1340
1341    if basename_out == None:
1342        netcdfname = root + '.dem'
1343    else:
1344        netcdfname = basename_out + '.dem'
1345
1346    if verbose: print 'Store to NetCDF file %s' %netcdfname
1347    # NetCDF file definition
1348    fid = NetCDFFile(netcdfname, 'w')
1349
1350    #Create new file
1351    fid.institution = 'Geoscience Australia'
1352    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
1353                      'of spatial point data'
1354
1355    fid.ncols = ncols
1356    fid.nrows = nrows
1357    fid.xllcorner = xllcorner
1358    fid.yllcorner = yllcorner
1359    fid.cellsize = cellsize
1360    fid.NODATA_value = NODATA_value
1361
1362    fid.zone = zone
1363    fid.false_easting = false_easting
1364    fid.false_northing = false_northing
1365    fid.projection = projection
1366    fid.datum = datum
1367    fid.units = units
1368
1369
1370    # dimension definitions
1371    fid.createDimension('number_of_rows', nrows)
1372    fid.createDimension('number_of_columns', ncols)
1373
1374    # variable definitions
1375    fid.createVariable('elevation', Float, ('number_of_rows',
1376                                            'number_of_columns'))
1377
1378    # Get handles to the variables
1379    elevation = fid.variables['elevation']
1380
1381    #Store data
1382    for i, line in enumerate(lines[6:]):
1383        fields = line.split()
1384        if verbose: print 'Processing row %d of %d' %(i, nrows)
1385
1386        elevation[i, :] = array([float(x) for x in fields])
1387
1388    fid.close()
1389
1390
1391
1392def ferret2sww(basename_in, basename_out = None,
1393               verbose = False,
1394               minlat = None, maxlat = None,
1395               minlon = None, maxlon = None,
1396               mint = None, maxt = None, mean_stage = 0,
1397               origin = None, zscale = 1,
1398               fail_on_NaN = True,
1399               NaN_filler = 0,
1400               elevation = None,
1401               inverted_bathymetry = False
1402               ): #FIXME: Bathymetry should be obtained
1403                                  #from MOST somehow.
1404                                  #Alternatively from elsewhere
1405                                  #or, as a last resort,
1406                                  #specified here.
1407                                  #The value of -100 will work
1408                                  #for the Wollongong tsunami
1409                                  #scenario but is very hacky
1410    """Convert 'Ferret' NetCDF format for wave propagation to
1411    sww format native to pyvolution.
1412
1413    Specify only basename_in and read files of the form
1414    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1415    relative height, x-velocity and y-velocity, respectively.
1416
1417    Also convert latitude and longitude to UTM. All coordinates are
1418    assumed to be given in the GDA94 datum.
1419
1420    min's and max's: If omitted - full extend is used.
1421    To include a value min may equal it, while max must exceed it.
1422    Lat and lon are assuemd to be in decimal degrees
1423
1424    origin is a 3-tuple with geo referenced
1425    UTM coordinates (zone, easting, northing)
1426
1427    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
1428    which means that longitude is the fastest
1429    varying dimension (row major order, so to speak)
1430
1431    ferret2sww uses grid points as vertices in a triangular grid
1432    counting vertices from lower left corner upwards, then right
1433    """ 
1434
1435    import os
1436    from Scientific.IO.NetCDF import NetCDFFile
1437    from Numeric import Float, Int, Int32, searchsorted, zeros, array
1438    precision = Float
1439
1440
1441    #Get NetCDF data
1442    if verbose: print 'Reading files %s_*.nc' %basename_in
1443    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
1444    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
1445    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
1446    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)   
1447
1448    if basename_out is None:
1449        swwname = basename_in + '.sww'
1450    else:
1451        swwname = basename_out + '.sww'       
1452
1453    times = file_h.variables['TIME']
1454    latitudes = file_h.variables['LAT']
1455    longitudes = file_h.variables['LON']
1456
1457    if mint == None:
1458        jmin = 0
1459    else:   
1460        jmin = searchsorted(times, mint)
1461   
1462    if maxt == None:
1463        jmax=len(times)
1464    else:   
1465        jmax = searchsorted(times, maxt)
1466       
1467    if minlat == None:
1468        kmin=0
1469    else:
1470        kmin = searchsorted(latitudes, minlat)
1471
1472    if maxlat == None:
1473        kmax = len(latitudes)
1474    else:
1475        kmax = searchsorted(latitudes, maxlat)
1476
1477    if minlon == None:
1478        lmin=0
1479    else:   
1480        lmin = searchsorted(longitudes, minlon)
1481   
1482    if maxlon == None:
1483        lmax = len(longitudes)
1484    else:
1485        lmax = searchsorted(longitudes, maxlon)           
1486
1487
1488
1489    times = times[jmin:jmax]       
1490    latitudes = latitudes[kmin:kmax]
1491    longitudes = longitudes[lmin:lmax]
1492
1493
1494    if verbose: print 'cropping'
1495    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1496    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
1497    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
1498    elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
1499
1500#    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
1501#        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
1502#    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
1503#        from Numeric import asarray
1504#        elevations=elevations.tolist()
1505#        elevations.reverse()
1506#        elevations=asarray(elevations)
1507#    else:
1508#        from Numeric import asarray
1509#        elevations=elevations.tolist()
1510#        elevations.reverse()
1511#        elevations=asarray(elevations)
1512#        'print hmmm'
1513
1514
1515
1516    #Get missing values
1517    nan_ha = file_h.variables['HA'].missing_value[0]
1518    nan_ua = file_u.variables['UA'].missing_value[0]
1519    nan_va = file_v.variables['VA'].missing_value[0]
1520    if hasattr(file_e.variables['ELEVATION'],'missing_value'):
1521        nan_e  = file_e.variables['ELEVATION'].missing_value[0]
1522    else:
1523        nan_e = None
1524   
1525    #Cleanup
1526    from Numeric import sometrue
1527   
1528    missing = (amplitudes == nan_ha)
1529    if sometrue (missing):
1530        if fail_on_NaN:
1531            msg = 'NetCDFFile %s contains missing values'\
1532                  %(basename_in+'_ha.nc')
1533            raise msg
1534        else:
1535            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
1536
1537    missing = (uspeed == nan_ua)
1538    if sometrue (missing):
1539        if fail_on_NaN:
1540            msg = 'NetCDFFile %s contains missing values'\
1541                  %(basename_in+'_ua.nc')
1542            raise msg
1543        else:
1544            uspeed = uspeed*(missing==0) + missing*NaN_filler
1545
1546    missing = (vspeed == nan_va)
1547    if sometrue (missing):
1548        if fail_on_NaN:
1549            msg = 'NetCDFFile %s contains missing values'\
1550                  %(basename_in+'_va.nc')
1551            raise msg
1552        else:
1553            vspeed = vspeed*(missing==0) + missing*NaN_filler
1554
1555
1556    missing = (elevations == nan_e)
1557    if sometrue (missing):
1558        if fail_on_NaN:
1559            msg = 'NetCDFFile %s contains missing values'\
1560                  %(basename_in+'_e.nc')
1561            raise msg
1562        else:
1563            elevations = elevations*(missing==0) + missing*NaN_filler
1564
1565    #######
1566
1567               
1568
1569    number_of_times = times.shape[0]
1570    number_of_latitudes = latitudes.shape[0]
1571    number_of_longitudes = longitudes.shape[0]
1572
1573    assert amplitudes.shape[0] == number_of_times
1574    assert amplitudes.shape[1] == number_of_latitudes
1575    assert amplitudes.shape[2] == number_of_longitudes         
1576
1577    if verbose:
1578        print '------------------------------------------------'
1579        print 'Statistics:'
1580        print '  Extent (lat/lon):'
1581        print '    lat in [%f, %f], len(lat) == %d'\
1582              %(min(latitudes.flat), max(latitudes.flat),
1583                len(latitudes.flat))
1584        print '    lon in [%f, %f], len(lon) == %d'\
1585              %(min(longitudes.flat), max(longitudes.flat),
1586                len(longitudes.flat))
1587        print '    t in [%f, %f], len(t) == %d'\
1588              %(min(times.flat), max(times.flat), len(times.flat))
1589       
1590        q = amplitudes.flat       
1591        name = 'Amplitudes (ha) [cm]'
1592        print '  %s in [%f, %f]' %(name, min(q), max(q))
1593       
1594        q = uspeed.flat
1595        name = 'Speeds (ua) [cm/s]'       
1596        print '  %s in [%f, %f]' %(name, min(q), max(q))
1597       
1598        q = vspeed.flat
1599        name = 'Speeds (va) [cm/s]'               
1600        print '  %s in [%f, %f]' %(name, min(q), max(q))               
1601
1602        q = elevations.flat
1603        name = 'Elevations (e) [m]'               
1604        print '  %s in [%f, %f]' %(name, min(q), max(q))               
1605
1606
1607    #print number_of_latitudes, number_of_longitudes
1608    number_of_points = number_of_latitudes*number_of_longitudes
1609    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
1610   
1611
1612    file_h.close()
1613    file_u.close()
1614    file_v.close()   
1615    file_e.close()   
1616
1617
1618    # NetCDF file definition
1619    outfile = NetCDFFile(swwname, 'w')
1620       
1621    #Create new file
1622    outfile.institution = 'Geoscience Australia'
1623    outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
1624                          %(basename_in + '_ha.nc',
1625                            basename_in + '_ua.nc',
1626                            basename_in + '_va.nc',
1627                            basename_in + '_e.nc')
1628
1629
1630    #For sww compatibility
1631    outfile.smoothing = 'Yes' 
1632    outfile.order = 1
1633           
1634    #Start time in seconds since the epoch (midnight 1/1/1970)
1635    outfile.starttime = starttime = times[0]
1636    times = times - starttime  #Store relative times
1637   
1638    # dimension definitions
1639    outfile.createDimension('number_of_volumes', number_of_volumes)
1640
1641    outfile.createDimension('number_of_vertices', 3)
1642    outfile.createDimension('number_of_points', number_of_points)
1643                           
1644               
1645    #outfile.createDimension('number_of_timesteps', len(times))
1646    outfile.createDimension('number_of_timesteps', len(times))
1647
1648    # variable definitions
1649    outfile.createVariable('x', precision, ('number_of_points',))
1650    outfile.createVariable('y', precision, ('number_of_points',))
1651    outfile.createVariable('elevation', precision, ('number_of_points',))
1652           
1653    #FIXME: Backwards compatibility
1654    outfile.createVariable('z', precision, ('number_of_points',))
1655    #################################
1656       
1657    outfile.createVariable('volumes', Int, ('number_of_volumes',
1658                                            'number_of_vertices'))
1659   
1660    outfile.createVariable('time', precision,
1661                           ('number_of_timesteps',))
1662       
1663    outfile.createVariable('stage', precision,
1664                           ('number_of_timesteps',
1665                            'number_of_points'))
1666
1667    outfile.createVariable('xmomentum', precision,
1668                           ('number_of_timesteps',
1669                            'number_of_points'))
1670   
1671    outfile.createVariable('ymomentum', precision,
1672                           ('number_of_timesteps',
1673                            'number_of_points'))
1674
1675
1676    #Store
1677    from coordinate_transforms.redfearn import redfearn
1678    x = zeros(number_of_points, Float)  #Easting
1679    y = zeros(number_of_points, Float)  #Northing
1680
1681
1682    if verbose: print 'Making triangular grid'
1683    #Check zone boundaries
1684    refzone, _, _ = redfearn(latitudes[0],longitudes[0])   
1685
1686    vertices = {}
1687    i = 0   
1688    for k, lat in enumerate(latitudes):       #Y direction
1689        for l, lon in enumerate(longitudes):  #X direction
1690
1691            vertices[l,k] = i
1692
1693            zone, easting, northing = redfearn(lat,lon)               
1694
1695            msg = 'Zone boundary crossed at longitude =', lon
1696            #assert zone == refzone, msg
1697            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
1698            x[i] = easting
1699            y[i] = northing
1700            i += 1
1701
1702
1703    #Construct 2 triangles per 'rectangular' element
1704    volumes = []
1705    for l in range(number_of_longitudes-1):    #X direction
1706        for k in range(number_of_latitudes-1): #Y direction
1707            v1 = vertices[l,k+1]
1708            v2 = vertices[l,k]           
1709            v3 = vertices[l+1,k+1]           
1710            v4 = vertices[l+1,k]           
1711
1712            volumes.append([v1,v2,v3]) #Upper element
1713            volumes.append([v4,v3,v2]) #Lower element           
1714
1715    volumes = array(volumes)       
1716
1717    if origin == None:
1718        zone = refzone       
1719        xllcorner = min(x)
1720        yllcorner = min(y)
1721    else:
1722        zone = origin[0]
1723        xllcorner = origin[1]
1724        yllcorner = origin[2]               
1725
1726   
1727    outfile.xllcorner = xllcorner
1728    outfile.yllcorner = yllcorner
1729    outfile.zone = zone
1730
1731
1732    if elevation is not None:
1733        z = elevation
1734    else:
1735        if inverted_bathymetry:
1736            z = -1*elevations
1737        else:
1738            z = elevations
1739        #FIXME: z should be obtained from MOST and passed in here
1740 
1741    from Numeric import resize
1742    z = resize(z,outfile.variables['z'][:].shape)
1743    outfile.variables['x'][:] = x - xllcorner
1744    outfile.variables['y'][:] = y - yllcorner
1745    outfile.variables['z'][:] = z
1746    outfile.variables['elevation'][:] = z  #FIXME HACK
1747    outfile.variables['time'][:] = times   #Store time relative
1748    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
1749   
1750
1751
1752    #Time stepping
1753    stage = outfile.variables['stage']
1754    xmomentum = outfile.variables['xmomentum']
1755    ymomentum = outfile.variables['ymomentum'] 
1756
1757    if verbose: print 'Converting quantities'
1758    n = len(times)
1759    for j in range(n):
1760        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
1761        i = 0
1762        for k in range(number_of_latitudes):      #Y direction
1763            for l in range(number_of_longitudes): #X direction
1764                w = zscale*amplitudes[j,k,l]/100 + mean_stage
1765                stage[j,i] = w
1766                h = w - z[i]
1767                xmomentum[j,i] = uspeed[j,k,l]/100*h
1768                ymomentum[j,i] = vspeed[j,k,l]/100*h
1769                i += 1
1770
1771
1772    if verbose:
1773        x = outfile.variables['x'][:]
1774        y = outfile.variables['y'][:]       
1775        print '------------------------------------------------'
1776        print 'Statistics of output file:'
1777        print '  Name: %s' %swwname
1778        print '  Reference:'
1779        print '    Lower left corner: [%f, %f]'\
1780              %(xllcorner, yllcorner)
1781        print '    Start time: %f' %starttime
1782        print '  Extent:'
1783        print '    x [m] in [%f, %f], len(x) == %d'\
1784              %(min(x.flat), max(x.flat), len(x.flat))
1785        print '    y [m] in [%f, %f], len(y) == %d'\
1786              %(min(y.flat), max(y.flat), len(y.flat))
1787        print '    t [s] in [%f, %f], len(t) == %d'\
1788              %(min(times), max(times), len(times))
1789        print '  Quantities [SI units]:'
1790        for name in ['stage', 'xmomentum', 'ymomentum']:
1791            q = outfile.variables[name][:].flferret2swwat
1792            print '    %s in [%f, %f]' %(name, min(q), max(q))
1793           
1794       
1795
1796       
1797    outfile.close()               
1798   
1799
1800
1801def extent_sww(file_name):
1802    """
1803    Read in an sww file.
1804
1805    Input;
1806    file_name - the sww file
1807
1808    Output;
1809    z - Vector of bed elevation
1810    volumes - Array.  Each row has 3 values, representing
1811    the vertices that define the volume
1812    time - Vector of the times where there is stage information
1813    stage - array with respect to time and vertices (x,y)
1814    """
1815
1816
1817    from Scientific.IO.NetCDF import NetCDFFile
1818
1819    #Check contents
1820    #Get NetCDF
1821    fid = NetCDFFile(file_name, 'r')
1822
1823    # Get the variables
1824    x = fid.variables['x'][:]
1825    y = fid.variables['y'][:]
1826    stage = fid.variables['stage'][:]
1827    #print "stage",stage
1828    #print "stage.shap",stage.shape
1829    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
1830    #print "min(stage)",min(stage)
1831
1832    fid.close()
1833
1834    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
1835
1836
1837def sww2domain(filename,t=None,fail_if_NaN=True,NaN_filler=0,verbose = True):
1838    """
1839    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
1840
1841    If the sww has stages, but not
1842    """
1843    NaN=9.969209968386869e+036
1844    #initialise NaN.
1845
1846    from Scientific.IO.NetCDF import NetCDFFile
1847    from domain import Domain
1848    from Numeric import asarray, transpose
1849
1850    if verbose: print 'Reading from ', filename
1851    fid = NetCDFFile(filename, 'r')    #Open existing file for read
1852    time = fid.variables['time']       #Timesteps
1853    if t is None:
1854        t = time[-1]
1855    time_interp = get_time_interp(time,t)
1856
1857    # Get the variables as Numeric arrays
1858    x = fid.variables['x'][:]             #x-coordinates of vertices
1859    y = fid.variables['y'][:]             #y-coordinates of vertices
1860    elevation = fid.variables['elevation']     #Elevation
1861    stage = fid.variables['stage']     #Water level
1862    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
1863    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
1864
1865    xllcorner = fid.xllcorner[0]
1866    yllcorner = fid.yllcorner[0]
1867    starttime = fid.starttime[0]
1868    zone = fid.zone
1869    volumes = fid.variables['volumes'][:] #Connectivity
1870    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
1871
1872    conserved_quantities = []
1873    interpolated_quantities = {}
1874    other_quantities = []
1875
1876    if verbose: print '    interpolating quantities'
1877    for quantity in fid.variables.keys():
1878        dimensions = fid.variables[quantity].dimensions
1879        if 'number_of_timesteps' in dimensions:
1880            conserved_quantities.append(quantity)
1881            interpolated_quantities[quantity]=\
1882                  interpolated_quantity(fid.variables[quantity][:],time_interp)
1883        else: other_quantities.append(quantity)
1884
1885    other_quantities.remove('x')
1886    other_quantities.remove('y')
1887    other_quantities.remove('z')
1888    other_quantities.remove('volumes')
1889
1890    conserved_quantities.remove('time')
1891
1892    if verbose: print '    building domain'
1893    domain = Domain(coordinates, volumes,\
1894                    conserved_quantities = conserved_quantities,\
1895                    other_quantities = other_quantities,zone=zone,\
1896                    xllcorner=xllcorner, yllcorner=yllcorner)
1897    domain.starttime=starttime
1898    domain.time=t
1899    for quantity in other_quantities:
1900        try:
1901            NaN = fid.variables[quantity].missing_value
1902        except:
1903            pass #quantity has no missing_value number
1904        X = fid.variables[quantity][:]
1905        if verbose:
1906            print '       ',quantity
1907            print '        NaN =',NaN
1908            print '        max(X)'
1909            print '       ',max(X)
1910            print '        max(X)==NaN'
1911            print '       ',max(X)==NaN
1912            print ''
1913        if (max(X)==NaN) or (min(X)==NaN):
1914            if fail_if_NaN:
1915                msg = 'quantity "%s" contains no_data entry'%quantity
1916                raise msg
1917            else:
1918                data = (X<>NaN)
1919                X = (X*data)+(data==0)*NaN_filler
1920        domain.set_quantity(quantity,X)
1921#
1922    for quantity in conserved_quantities:
1923        try:
1924            NaN = fid.variables[quantity].missing_value
1925        except:
1926            pass #quantity has no missing_value number
1927        X = interpolated_quantities[quantity]
1928        if verbose:
1929            print '       ',quantity
1930            print '        NaN =',NaN
1931            print '        max(X)'
1932            print '       ',max(X)
1933            print '        max(X)==NaN'
1934            print '       ',max(X)==NaN
1935            print ''
1936        if (max(X)==NaN) or (min(X)==NaN):
1937            if fail_if_NaN:
1938                msg = 'quantity "%s" contains no_data entry'%quantity
1939                raise msg
1940            else:
1941                data = (X<>NaN)
1942                X = (X*data)+(data==0)*NaN_filler
1943        domain.set_quantity(quantity,X)
1944    fid.close()
1945    return domain
1946
1947def interpolated_quantity(saved_quantity,time_interp):
1948       
1949    #given an index and ratio, interpolate quantity with respect to time.
1950    index,ratio = time_interp
1951    Q = saved_quantity 
1952    if ratio > 0:
1953        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
1954    else:
1955         q = Q[index]
1956    #Return vector of interpolated values
1957    return q
1958
1959def get_time_interp(time,t=None):
1960    #Finds the ratio and index for time interpolation.
1961    #It is borrowed from previous pyvolution code.
1962    if t is None:
1963        t=time[-1]
1964        index = -1
1965        ratio = 0.
1966    else:
1967        T = time
1968        tau = t
1969        index=0
1970        msg = 'Time interval derived from file %s [%s:%s]'\
1971            %('FIXMEfilename', T[0], T[-1])
1972        msg += ' does not match model time: %s' %tau
1973        if tau < time[0]: raise msg
1974        if tau > time[-1]: raise msg
1975        while tau > time[index]: index += 1
1976        while tau < time[index]: index -= 1
1977        if tau == time[index]:
1978            #Protect against case where tau == time[-1] (last time)
1979            # - also works in general when tau == time[i]
1980            ratio = 0
1981        else:
1982            #t is now between index and index+1           
1983            ratio = (tau - time[index])/(time[index+1] - time[index])
1984    return (index,ratio)
1985
1986
1987
1988#OBSOLETE STUFF
1989#Native checkpoint format.
1990#Information needed to recreate a state is preserved
1991#FIXME: Rethink and maybe use netcdf format
1992def cpt_variable_writer(filename, t, v0, v1, v2):
1993    """Store all conserved quantities to file
1994    """
1995
1996    M, N = v0.shape
1997
1998    FN = create_filename(filename, 'cpt', M, t)
1999    #print 'Writing to %s' %FN
2000
2001    fid = open(FN, 'w')
2002    for i in range(M):
2003        for j in range(N):
2004            fid.write('%.16e ' %v0[i,j])
2005        for j in range(N):
2006            fid.write('%.16e ' %v1[i,j])
2007        for j in range(N):
2008            fid.write('%.16e ' %v2[i,j])
2009
2010        fid.write('\n')
2011    fid.close()
2012
2013
2014def cpt_variable_reader(filename, t, v0, v1, v2):
2015    """Store all conserved quantities to file
2016    """
2017
2018    M, N = v0.shape
2019
2020    FN = create_filename(filename, 'cpt', M, t)
2021    #print 'Reading from %s' %FN
2022
2023    fid = open(FN)
2024
2025
2026    for i in range(M):
2027        values = fid.readline().split() #Get one line
2028
2029        for j in range(N):
2030            v0[i,j] = float(values[j])
2031            v1[i,j] = float(values[3+j])
2032            v2[i,j] = float(values[6+j])
2033
2034    fid.close()
2035
2036def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
2037    """Writes x,y,z,z,z coordinates of triangles constituting the bed
2038    elevation.
2039    Not in use pt
2040    """
2041
2042    M, N = v0.shape
2043
2044    print X0
2045    import sys; sys.exit()
2046    FN = create_filename(filename, 'cpt', M)
2047    print 'Writing to %s' %FN
2048
2049    fid = open(FN, 'w')
2050    for i in range(M):
2051        for j in range(2):
2052            fid.write('%.16e ' %X0[i,j])   #x, y
2053        for j in range(N):
2054            fid.write('%.16e ' %v0[i,j])       #z,z,z,
2055
2056        for j in range(2):
2057            fid.write('%.16e ' %X1[i,j])   #x, y
2058        for j in range(N):
2059            fid.write('%.16e ' %v1[i,j])
2060
2061        for j in range(2):
2062            fid.write('%.16e ' %X2[i,j])   #x, y
2063        for j in range(N):
2064            fid.write('%.16e ' %v2[i,j])
2065
2066        fid.write('\n')
2067    fid.close()
2068
2069
2070
2071#Function for storing out to e.g. visualisation
2072#FIXME: Do we want this?
2073#FIXME: Not done yet for this version
2074def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
2075    """Writes x,y,z coordinates of triangles constituting the bed elevation.
2076    """
2077
2078    M, N = v0.shape
2079
2080    FN = create_filename(filename, 'dat', M)
2081    #print 'Writing to %s' %FN
2082
2083    fid = open(FN, 'w')
2084    for i in range(M):
2085        for j in range(2):
2086            fid.write('%f ' %X0[i,j])   #x, y
2087        fid.write('%f ' %v0[i,0])       #z
2088
2089        for j in range(2):
2090            fid.write('%f ' %X1[i,j])   #x, y
2091        fid.write('%f ' %v1[i,0])       #z
2092
2093        for j in range(2):
2094            fid.write('%f ' %X2[i,j])   #x, y
2095        fid.write('%f ' %v2[i,0])       #z
2096
2097        fid.write('\n')
2098    fid.close()
2099
2100
2101
2102def dat_variable_writer(filename, t, v0, v1, v2):
2103    """Store water height to file
2104    """
2105
2106    M, N = v0.shape
2107
2108    FN = create_filename(filename, 'dat', M, t)
2109    #print 'Writing to %s' %FN
2110
2111    fid = open(FN, 'w')
2112    for i in range(M):
2113        fid.write('%.4f ' %v0[i,0])
2114        fid.write('%.4f ' %v1[i,0])
2115        fid.write('%.4f ' %v2[i,0])
2116
2117        fid.write('\n')
2118    fid.close()
2119
2120
2121def read_sww(filename):
2122    """Read sww Net CDF file containing Shallow Water Wave simulation
2123
2124    The integer array volumes is of shape Nx3 where N is the number of
2125    triangles in the mesh.
2126
2127    Each entry in volumes is an index into the x,y arrays (the location).
2128
2129    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
2130    number_of_timesteps, number_of_points.
2131
2132    The momentum is not always stored.
2133   
2134    """
2135    from Scientific.IO.NetCDF import NetCDFFile
2136    print 'Reading from ', filename
2137    fid = NetCDFFile(filename, 'r')    #Open existing file for read
2138#latitude, longitude
2139    # Get the variables as Numeric arrays
2140    x = fid.variables['x']             #x-coordinates of vertices
2141    y = fid.variables['y']             #y-coordinates of vertices
2142    z = fid.variables['elevation']     #Elevation
2143    time = fid.variables['time']       #Timesteps
2144    stage = fid.variables['stage']     #Water level
2145    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
2146    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction   
2147
2148    volumes = fid.variables['volumes'] #Connectivity
Note: See TracBrowser for help on using the repository browser.