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

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

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