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

Last change on this file since 1393 was 1393, checked in by steve, 20 years ago

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