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

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

OKOKOKOK
File size: 37.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.dem: NetCDF representation of regular DEM data
19
20.tsh: ASCII format for storing meshes and associated boundary and region info
21
22
23#Not yet implemented
24.nc: Native ferret NetCDF format
25
26
27FIXME: What else
28
29"""
30
31from Numeric import concatenate       
32
33
34def make_filename(s):
35    """Transform argument string into a suitable filename
36    """
37
38    s = s.strip()
39    s = s.replace(' ', '_')
40    s = s.replace('(', '')
41    s = s.replace(')', '')
42    s = s.replace('__', '_')               
43       
44    return s
45
46
47def check_dir(path, verbose=None):
48    """Check that specified path exists.
49    If path does not exist it will be created if possible   
50
51    USAGE:
52       checkdir(path, verbose):
53
54    ARGUMENTS:
55        path -- Directory
56        verbose -- Flag verbose output (default: None)
57
58    RETURN VALUE:
59        Verified path including trailing separator
60       
61    """
62
63    import os, sys
64    import os.path
65
66    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
67        unix = 0
68    else:
69        unix = 1
70
71
72    if path[-1] != os.sep: 
73        path = path + os.sep  # Add separator for directories
74
75    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
76    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
77        try:
78            exitcode=os.mkdir(path)
79
80            # Change access rights if possible
81            #
82            if unix:
83                exitcode=os.system('chmod 775 '+path)
84            else:
85                pass  # FIXME: What about acces rights under Windows?
86           
87            if verbose: print 'MESSAGE: Directory', path, 'created.'
88         
89        except:
90            print 'WARNING: Directory', path, 'could not be created.'
91            if unix:
92                path = '/tmp/'
93            else:
94                path = 'C:'
95               
96            print 'Using directory %s instead' %path
97
98    return(path)
99
100   
101
102def del_dir(path):
103    """Recursively delete directory path and all its contents
104    """
105
106    import os
107
108    if os.path.isdir(path):
109        for file in os.listdir(path):
110            X = os.path.join(path, file)
111
112
113            if os.path.isdir(X) and not os.path.islink(X):
114                del_dir(X) 
115            else:
116                try:
117                    os.remove(X)
118                except:
119                    print "Could not remove file %s" %X
120
121        os.rmdir(path)
122
123
124       
125def create_filename(datadir, filename, format, size=None, time=None):
126
127    import os
128    #from config import data_dir
129           
130    FN = check_dir(datadir) + filename
131   
132    if size is not None:           
133        FN += '_size%d' %size
134
135    if time is not None:
136        FN += '_time%.2f' %time
137       
138    FN += '.' + format
139    return FN
140   
141
142def get_files(datadir, filename, format, size):
143    """Get all file (names) with gven name, size and format
144    """
145
146    import glob
147
148    import os
149    #from config import data_dir
150           
151    dir = check_dir(datadir)
152
153    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
154    return glob.glob(pattern)
155
156
157
158#Generic class for storing output to e.g. visualisation or checkpointing
159class Data_format:
160    """Generic interface to data formats
161    """
162
163       
164    def __init__(self, domain, extension, mode = 'w'):
165        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
166                                        '''   'w' (write)'''+\
167                                        '''   'r' (read)''' +\
168                                        '''   'a' (append)'''   
169
170        #Create filename
171        #self.filename = create_filename(domain.get_datadir(),
172        #       domain.get_name(), extension, len(domain))
173       
174        self.filename = create_filename(domain.get_datadir(), 
175                        domain.get_name(), extension)   
176               
177        self.timestep = 0
178        self.number_of_volumes = len(domain)
179        self.domain = domain
180   
181   
182    #FIXME: Should we have a general set_precision function?
183   
184
185
186#Class for storing output to e.g. visualisation
187class Data_format_sww(Data_format):
188    """Interface to native NetCDF format (.sww)
189    """
190
191       
192    def __init__(self, domain, mode = 'w'):
193        from Scientific.IO.NetCDF import NetCDFFile
194        from Numeric import Int, Float, Float32
195
196        self.precision = Float32 #Use single precision
197       
198        Data_format.__init__(self, domain, 'sww', mode)
199
200
201        # NetCDF file definition
202        fid = NetCDFFile(self.filename, mode)
203       
204        if mode == 'w':
205           
206            #Create new file
207            fid.institution = 'Geoscience Australia'
208            fid.description = 'Output from pyvolution suitable for plotting'
209           
210            if domain.smooth:
211                fid.smoothing = 'Yes'
212            else:   
213                fid.smoothing = 'No'
214
215            fid.order = domain.default_order   
216           
217            #Start time in seconds since the epoch (midnight 1/1/1970)
218            fid.starttime = domain.starttime
219
220            # dimension definitions
221            fid.createDimension('number_of_volumes', self.number_of_volumes) 
222            fid.createDimension('number_of_vertices', 3)
223
224            if domain.smooth is True:
225                fid.createDimension('number_of_points', len(domain.vertexlist))
226            else:   
227                fid.createDimension('number_of_points', 3*self.number_of_volumes)
228               
229            fid.createDimension('number_of_timesteps', None) #extensible
230
231            # variable definitions
232            fid.createVariable('x', self.precision, ('number_of_points',))
233            fid.createVariable('y', self.precision, ('number_of_points',))
234            fid.createVariable('elevation', self.precision, ('number_of_points',))
235           
236            #FIXME: Backwards compatibility
237            fid.createVariable('z', self.precision, ('number_of_points',))
238            #################################
239       
240            fid.createVariable('volumes', Int, ('number_of_volumes',
241                                                'number_of_vertices'))
242
243            fid.createVariable('time', self.precision,
244                               ('number_of_timesteps',))
245       
246            fid.createVariable('stage', self.precision,
247                               ('number_of_timesteps',
248                                'number_of_points'))
249
250            fid.createVariable('xmomentum', self.precision,
251                               ('number_of_timesteps',
252                                'number_of_points'))
253
254            fid.createVariable('ymomentum', self.precision,
255                               ('number_of_timesteps',
256                                'number_of_points'))
257
258        #Close
259        fid.close()               
260
261
262    def store_connectivity(self):
263        """Specialisation of store_connectivity for net CDF format
264
265        Writes x,y,z coordinates of triangles constituting
266        the bed elevation.
267        """
268
269        from Scientific.IO.NetCDF import NetCDFFile
270
271        from Numeric import concatenate, Int
272
273        domain = self.domain
274
275        #Get NetCDF
276        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
277       
278        # Get the variables
279        x = fid.variables['x']
280        y = fid.variables['y']
281        z = fid.variables['elevation']
282       
283        volumes = fid.variables['volumes']
284
285        # Get X, Y and bed elevation Z
286        Q = domain.quantities['elevation']
287        X,Y,Z,V = Q.get_vertex_values(xy=True,
288                                      precision = self.precision)
289                                     
290
291
292        x[:] = X.astype(self.precision)
293        y[:] = Y.astype(self.precision)
294        z[:] = Z.astype(self.precision)
295       
296        #FIXME: Backwards compatibility
297        z = fid.variables['z']
298        z[:] = Z.astype(self.precision) 
299        ################################
300
301        volumes[:] = V.astype(volumes.typecode())
302       
303        #Close
304        fid.close()
305
306
307    def store_timestep(self, names):
308        """Store time and named quantities to file
309        """
310        from Scientific.IO.NetCDF import NetCDFFile
311        import types
312        from time import sleep
313
314               
315        #Get NetCDF
316        retries = 0
317        file_open = False
318        while not file_open and retries < 10:
319            try:       
320                fid = NetCDFFile(self.filename, 'a') #Open existing file
321            except IOError:
322                #This could happen if someone was reading the file.
323                #In that case, wait a while and try again
324                msg = 'Warning (store_timestep): File %s could not be opened'\
325                %self.filename
326                msg += ' - trying again'
327                print msg
328                retries += 1
329                sleep(1)
330            else:
331               file_open = True
332               
333        if not file_open:
334            msg = 'File %s could not be opened for append' %self.filename
335            raise msg
336               
337           
338        domain = self.domain
339       
340        # Get the variables
341        time = fid.variables['time']
342        stage = fid.variables['stage']
343        xmomentum = fid.variables['xmomentum']
344        ymomentum = fid.variables['ymomentum']                       
345        i = len(time)
346
347        #Store time
348        time[i] = self.domain.time
349
350
351        if type(names) not in [types.ListType, types.TupleType]:
352            names = [names]
353
354        for name in names:   
355            # Get quantity
356            Q = domain.quantities[name]
357            A,V = Q.get_vertex_values(xy=False,
358                                      precision = self.precision)
359
360            #FIXME: Make this more general and rethink naming
361            if name == 'stage':
362                stage[i,:] = A.astype(self.precision)
363            elif name == 'xmomentum':
364                xmomentum[i,:] = A.astype(self.precision)
365            elif name == 'ymomentum':
366                ymomentum[i,:] = A.astype(self.precision)   
367                                     
368            #As in....
369            #eval( name + '[i,:] = A.astype(self.precision)' )
370            #FIXME: But we need a UNIT test for that before refactoring
371           
372           
373
374        #Flush and close
375        fid.sync()
376        fid.close()
377
378
379#Class for handling checkpoints data
380class Data_format_cpt(Data_format):
381    """Interface to native NetCDF format (.cpt)
382    """
383
384       
385    def __init__(self, domain, mode = 'w'):
386        from Scientific.IO.NetCDF import NetCDFFile
387        from Numeric import Int, Float, Float
388
389        self.precision = Float #Use full precision
390       
391        Data_format.__init__(self, domain, 'sww', mode)
392
393
394        # NetCDF file definition
395        fid = NetCDFFile(self.filename, mode)
396       
397        if mode == 'w':
398            #Create new file
399            fid.institution = 'Geoscience Australia'
400            fid.description = 'Checkpoint data'
401            #fid.smooth = domain.smooth               
402            fid.order = domain.default_order   
403
404            # dimension definitions
405            fid.createDimension('number_of_volumes', self.number_of_volumes) 
406            fid.createDimension('number_of_vertices', 3)
407
408            #Store info at all vertices (no smoothing)
409            fid.createDimension('number_of_points', 3*self.number_of_volumes)
410            fid.createDimension('number_of_timesteps', None) #extensible
411
412            # variable definitions
413
414            #Mesh
415            fid.createVariable('x', self.precision, ('number_of_points',))
416            fid.createVariable('y', self.precision, ('number_of_points',))
417
418       
419            fid.createVariable('volumes', Int, ('number_of_volumes',
420                                                'number_of_vertices'))
421
422            fid.createVariable('time', self.precision,
423                               ('number_of_timesteps',))
424
425            #Allocate space for all quantities
426            for name in domain.quantities.keys():
427                fid.createVariable(name, self.precision,
428                                   ('number_of_timesteps',
429                                    'number_of_points'))
430
431        #Close
432        fid.close()               
433
434
435    def store_checkpoint(self):
436        """
437        Write x,y coordinates of triangles.
438        Write connectivity (
439        constituting
440        the bed elevation.
441        """
442
443        from Scientific.IO.NetCDF import NetCDFFile
444
445        from Numeric import concatenate
446
447        domain = self.domain
448
449        #Get NetCDF
450        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
451       
452        # Get the variables
453        x = fid.variables['x']
454        y = fid.variables['y']
455       
456        volumes = fid.variables['volumes']
457
458        # Get X, Y and bed elevation Z
459        Q = domain.quantities['elevation']
460        X,Y,Z,V = Q.get_vertex_values(xy=True,
461                                      precision = self.precision)
462                                     
463
464
465        x[:] = X.astype(self.precision)
466        y[:] = Y.astype(self.precision)
467        z[:] = Z.astype(self.precision)
468
469        volumes[:] = V
470       
471        #Close
472        fid.close()
473
474
475    def store_timestep(self, name):
476        """Store time and named quantity to file
477        """
478        from Scientific.IO.NetCDF import NetCDFFile
479        from time import sleep
480               
481        #Get NetCDF
482        retries = 0
483        file_open = False
484        while not file_open and retries < 10:
485            try:       
486                fid = NetCDFFile(self.filename, 'a') #Open existing file
487            except IOError:
488                #This could happen if someone was reading the file.
489                #In that case, wait a while and try again
490                msg = 'Warning (store_timestep): File %s could not be opened'\
491                %self.filename
492                msg += ' - trying again'
493                print msg
494                retries += 1
495                sleep(1)
496            else:
497               file_open = True
498               
499        if not file_open:
500            msg = 'File %s could not be opened for append' %self.filename
501            raise msg
502               
503           
504        domain = self.domain
505       
506        # Get the variables
507        time = fid.variables['time']
508        stage = fid.variables['stage']       
509        i = len(time)
510
511        #Store stage
512        time[i] = self.domain.time
513
514        # Get quantity
515        Q = domain.quantities[name]
516        A,V = Q.get_vertex_values(xy=False,
517                                  precision = self.precision)
518       
519        stage[i,:] = A.astype(self.precision)
520
521        #Flush and close
522        fid.sync()
523        fid.close()
524
525
526
527
528
529
530
531#Function for storing xya output
532#FIXME Not done yet for this version
533class Data_format_xya(Data_format):
534    """Generic interface to data formats
535    """
536
537
538    def __init__(self, domain, mode = 'w'):
539        from Scientific.IO.NetCDF import NetCDFFile
540        from Numeric import Int, Float, Float32
541
542        self.precision = Float32 #Use single precision
543       
544        Data_format.__init__(self, domain, 'xya', mode)
545       
546   
547   
548    #FIXME -This is the old xya format 
549    def store_all(self):
550        """Specialisation of store all for xya format
551
552        Writes x,y,z coordinates of triangles constituting
553        the bed elevation.
554        """
555
556        from Numeric import concatenate
557
558        domain = self.domain
559       
560        fd = open(self.filename, 'w')
561
562
563        if domain.smooth is True:
564            number_of_points =  len(domain.vertexlist)
565        else:   
566            number_of_points = 3*self.number_of_volumes
567
568        numVertAttrib = 3 #Three attributes is what is assumed by the xya format
569
570        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
571                 " # <vertex #> <x> <y> [attributes]" + "\n")
572
573
574        # Get X, Y, bed elevation and friction (index=0,1)
575        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
576                                           indices = (0,1), precision = self.precision)
577
578        bed_eles = A[:,0]
579        fricts = A[:,1]
580
581        # Get stage (index=0)
582        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
583                                       indices = (0,), precision = self.precision)
584
585        stages = B[:,0]
586
587        #<vertex #> <x> <y> [attributes]
588        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
589                                               stages, fricts):
590
591            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
592            fd.write(s)
593           
594        #close
595        fd.close()
596
597
598    def store_timestep(self, t, V0, V1, V2):
599        """Store time, water heights (and momentums) to file
600        """
601        pass
602   
603
604#Auxiliary
605def write_obj(filename,x,y,z):
606    """Store x,y,z vectors into filename (obj format)
607       Vectors are assumed to have dimension (M,3) where
608       M corresponds to the number elements.
609       triangles are assumed to be disconnected
610
611       The three numbers in each vector correspond to three vertices,
612
613       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
614       
615    """
616    #print 'Writing obj to %s' % filename
617
618    import os.path
619
620    root, ext = os.path.splitext(filename)
621    if ext == '.obj':
622        FN = filename
623    else:
624        FN = filename + '.obj'
625       
626
627    outfile = open(FN, 'wb')
628    outfile.write("# Triangulation as an obj file\n")
629
630    M, N = x.shape
631    assert N==3  #Assuming three vertices per element
632
633    for i in range(M):
634        for j in range(N):
635            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
636
637    for i in range(M):
638        base = i*N
639        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
640
641    outfile.close()
642
643
644
645#Conversion routines
646def sww2obj(basefilename, size):
647    """Convert netcdf based data output to obj
648    """
649    from Scientific.IO.NetCDF import NetCDFFile       
650
651    from Numeric import Float, zeros
652
653    #Get NetCDF
654    FN = create_filename('.', basefilename, 'sww', size)
655    print 'Reading from ', FN
656    fid = NetCDFFile(FN, 'r')  #Open existing file for read
657
658       
659    # Get the variables
660    x = fid.variables['x']
661    y = fid.variables['y']
662    z = fid.variables['elevation']       
663    time = fid.variables['time']
664    stage = fid.variables['stage']
665       
666    M = size  #Number of lines
667    xx = zeros((M,3), Float)
668    yy = zeros((M,3), Float)
669    zz = zeros((M,3), Float)
670
671    for i in range(M):
672        for j in range(3):
673            xx[i,j] = x[i+j*M]
674            yy[i,j] = y[i+j*M]
675            zz[i,j] = z[i+j*M]
676       
677    #Write obj for bathymetry
678    FN = create_filename('.', basefilename, 'obj', size)       
679    write_obj(FN,xx,yy,zz)
680
681
682    #Now read all the data with variable information, combine with
683    #x,y info and store as obj
684
685    for k in range(len(time)):
686        t = time[k]
687        print 'Processing timestep %f' %t
688
689        for i in range(M):
690            for j in range(3):
691                zz[i,j] = stage[k,i+j*M]
692
693
694        #Write obj for variable data
695        #FN = create_filename(basefilename, 'obj', size, time=t)
696        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)   
697        write_obj(FN,xx,yy,zz)
698
699
700def dat2obj(basefilename):
701    """Convert line based data output to obj
702    FIXME: Obsolete?
703    """
704
705    import glob, os
706    from config import data_dir
707   
708   
709    #Get bathymetry and x,y's
710    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
711   
712    from Numeric import zeros, Float
713
714    M = len(lines)  #Number of lines
715    x = zeros((M,3), Float)
716    y = zeros((M,3), Float)
717    z = zeros((M,3), Float)
718
719    ##i = 0
720    for i, line in enumerate(lines):
721        tokens = line.split()
722        values = map(float,tokens)       
723
724        for j in range(3):
725            x[i,j] = values[j*3]
726            y[i,j] = values[j*3+1]
727            z[i,j] = values[j*3+2]       
728       
729        ##i += 1
730
731
732    #Write obj for bathymetry
733    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
734
735
736    #Now read all the data files with variable information, combine with
737    #x,y info
738    #and store as obj
739
740    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
741
742    for filename in files:
743        print 'Processing %s' % filename
744
745        lines = open(data_dir+os.sep+filename,'r').readlines()
746        assert len(lines) == M
747        root, ext = os.path.splitext(filename)
748
749        #Get time from filename
750        i0 = filename.find('_time=')
751        if i0 == -1:
752            #Skip bathymetry file
753            continue
754       
755        i0 += 6  #Position where time starts           
756        i1 = filename.find('.dat')         
757
758        if i1 > i0:
759            t = float(filename[i0:i1])
760        else:
761            raise 'Hmmmm'
762               
763       
764       
765        ##i = 0
766        for i, line in enumerate(lines):
767            tokens = line.split()
768            values = map(float,tokens)       
769
770            for j in range(3):
771                z[i,j] = values[j]       
772       
773            ##i += 1
774
775        #Write obj for variable data
776        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
777
778
779def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
780    """Read netcdf filename1, pick timesteps first:step:last and save to
781    nettcdf file filename2
782    """
783    from Scientific.IO.NetCDF import NetCDFFile       
784   
785    #Get NetCDF
786    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
787    outfile = NetCDFFile(filename2, 'w')  #Open new file
788   
789
790    #Copy dimensions
791    for d in infile.dimensions:
792        outfile.createDimension(d, infile.dimensions[d])
793
794    for name in infile.variables:
795        var = infile.variables[name]
796        outfile.createVariable(name, var.typecode(), var.dimensions)
797       
798   
799    #Copy the static variables
800    for name in infile.variables:
801        if name == 'time' or name == 'stage':
802            pass
803        else:
804            #Copy
805            outfile.variables[name][:] = infile.variables[name][:]           
806
807    #Copy selected timesteps
808    time = infile.variables['time']
809    stage = infile.variables['stage']
810
811    newtime = outfile.variables['time']
812    newstage = outfile.variables['stage']   
813
814    if last is None:
815        last = len(time)
816
817    selection = range(first, last, step)
818    for i, j in enumerate(selection):
819        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
820        newtime[i] = time[j]
821        newstage[i,:] = stage[j,:]       
822       
823    #Close
824    infile.close()
825    outfile.close()
826
827
828#Get data objects
829def get_dataobject(domain, mode='w'):
830    """Return instance of class of given format using filename
831    """
832
833    cls = eval('Data_format_%s' %domain.format)
834    return cls(domain, mode)
835
836def dem2pts(filename, verbose=False):
837    """Read Digitial Elevation model from the following NetCDF format (.dem)
838
839    Example:
840
841    ncols         3121
842    nrows         1800
843    xllcorner     722000
844    yllcorner     5893000
845    cellsize      25
846    NODATA_value  -9999
847    138.3698 137.4194 136.5062 135.5558 ..........
848
849    Convert to NetCDF pts format which is
850
851    points:  (Nx2) Float array
852    elevation: N Float array
853    """ 
854
855    import os
856    from Scientific.IO.NetCDF import NetCDFFile
857    from Numeric import Float, arrayrange, concatenate   
858
859    root, ext = os.path.splitext(filename)
860
861    #Get NetCDF
862    infile = NetCDFFile(filename, 'r')  #Open existing netcdf file for read
863
864    if verbose: print 'Reading DEM from %s' %filename
865   
866    ncols = infile.ncols[0]
867    nrows = infile.nrows[0]
868    xllcorner = infile.xllcorner[0]  #Easting of lower left corner
869    yllcorner = infile.yllcorner[0]  #Northing of lower left corner
870    cellsize = infile.cellsize[0]
871    NODATA_value = infile.NODATA_value[0]
872    dem_elevation = infile.variables['elevation']
873
874    #Get output file                               
875    xyaname = root + '.pts'
876    if verbose: print 'Store to NetCDF file %s' %xyaname
877    # NetCDF file definition
878    outfile = NetCDFFile(xyaname, 'w')
879       
880    #Create new file
881    outfile.institution = 'Geoscience Australia'
882    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
883                      'of spatial point data'
884
885    #Georeferencing
886    outfile.zone = 56   #FIXME: Must be read from somewhere. Talk to Don.
887    outfile.xllcorner = xllcorner #Easting of lower left corner
888    outfile.yllcorner = yllcorner #Northing of lower left corner
889   
890
891    #Grid info (FIXME: probably not going to be used, but heck)
892    outfile.ncols = ncols
893    outfile.nrows = nrows
894
895   
896    # dimension definitions
897    outfile.createDimension('number_of_points', nrows*ncols)   
898    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
899   
900    # variable definitions
901    outfile.createVariable('points', Float, ('number_of_points',
902                                             'number_of_dimensions'))
903    outfile.createVariable('elevation', Float, ('number_of_points',))
904
905    # Get handles to the variables
906    points = outfile.variables['points']
907    elevation = outfile.variables['elevation']
908
909    #Store data
910    #FIXME: Could perhaps be faster using array operations
911    for i in range(nrows):
912        if verbose: print 'Processing row %d of %d' %(i, nrows)
913       
914        y = (nrows-i)*cellsize           
915        for j in range(ncols):
916            index = i*ncols + j
917           
918            x = j*cellsize
919            points[index, :] = [x,y]
920            elevation[index] = dem_elevation[i, j]           
921
922    infile.close()           
923    outfile.close()
924
925                                     
926def convert_dem_from_ascii2netcdf(filename, verbose=False):
927    """Read Digitial Elevation model from the following ASCII format (.asc)
928
929    Example:
930
931    ncols         3121
932    nrows         1800
933    xllcorner     722000
934    yllcorner     5893000
935    cellsize      25
936    NODATA_value  -9999
937    138.3698 137.4194 136.5062 135.5558 ..........
938
939    Convert to NetCDF format (.dem) mimcing the ASCII format closely.
940    """ 
941
942    import os
943    from Scientific.IO.NetCDF import NetCDFFile
944    from Numeric import Float, array
945
946    root, ext = os.path.splitext(filename)
947    fid = open(filename)
948
949    if verbose: print 'Reading DEM from %s' %filename
950    lines = fid.readlines()
951    fid.close()
952
953    if verbose: print 'Got', len(lines), ' lines'
954   
955    ncols = int(lines[0].split()[1].strip())
956    nrows = int(lines[1].split()[1].strip())
957    xllcorner = float(lines[2].split()[1].strip())
958    yllcorner = float(lines[3].split()[1].strip())
959    cellsize = float(lines[4].split()[1].strip())
960    NODATA_value = int(lines[5].split()[1].strip())
961
962    assert len(lines) == nrows + 6
963
964    netcdfname = root + '.dem'
965    if verbose: print 'Store to NetCDF file %s' %netcdfname
966    # NetCDF file definition
967    fid = NetCDFFile(netcdfname, 'w')
968       
969    #Create new file
970    fid.institution = 'Geoscience Australia'
971    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
972                      'of spatial point data'
973
974    fid.ncols = ncols
975    fid.nrows = nrows
976    fid.xllcorner = xllcorner
977    fid.yllcorner = yllcorner
978    fid.cellsize = cellsize
979    fid.NODATA_value = NODATA_value
980
981    # dimension definitions
982    fid.createDimension('number_of_rows', nrows)
983    fid.createDimension('number_of_columns', ncols)           
984
985    # variable definitions
986    fid.createVariable('elevation', Float, ('number_of_rows',
987                                            'number_of_columns'))
988
989    # Get handles to the variables
990    elevation = fid.variables['elevation']
991
992    #Store data
993    for i, line in enumerate(lines[6:]):
994        fields = line.split()
995        if verbose: print 'Processing row %d of %d' %(i, nrows)
996
997        elevation[i, :] = array([float(x) for x in fields])
998
999    fid.close()
1000
1001   
1002
1003
1004
1005
1006
1007
1008
1009def ferret2sww(basefilename, verbose=False,
1010               minlat = None, maxlat =None,
1011               minlon = None, maxlon =None,
1012               mint = None, maxt = None, mean_stage = 0):
1013    """Convert 'Ferret' NetCDF format for wave propagation to
1014    sww format native to pyvolution.
1015
1016    Specify only basefilename and read files of the form
1017    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1018    relative height, x-velocity and y-velocity, respectively.
1019
1020    Also convert latitude and longitude to UTM. All coordinates are
1021    assumed to be given in the GDA94 datum.
1022
1023
1024    min's and max's: If omitted - full extend is used.
1025    To include a value min may equal it, while max must exceed it.
1026    """ 
1027
1028    import os
1029    from Scientific.IO.NetCDF import NetCDFFile
1030    from Numeric import Float, Int, searchsorted, zeros, array
1031    precision = Float
1032
1033
1034    #Get NetCDF data
1035    file_h = NetCDFFile(basefilename + '_ha.nc', 'r') #Wave amplitude (cm)
1036    file_u = NetCDFFile(basefilename + '_ua.nc', 'r') #Velocity (x) (cm/s)
1037    file_v = NetCDFFile(basefilename + '_va.nc', 'r') #Velocity (y) (cm/s)   
1038
1039    swwname = basefilename + '.sww'
1040
1041    times = file_h.variables['TIME']
1042    latitudes = file_h.variables['LAT']
1043    longitudes = file_h.variables['LON']
1044
1045    if mint == None:
1046        jmin = 0
1047    else:   
1048        jmin = searchsorted(times, mint)
1049   
1050    if maxt == None:
1051        jmax=len(times)
1052    else:   
1053        jmax = searchsorted(times, maxt)
1054       
1055    if minlat == None:
1056        kmin=0
1057    else:
1058        kmin = searchsorted(latitudes, minlat)
1059
1060    if maxlat == None:
1061        kmax = len(latitudes)
1062    else:
1063        kmax = searchsorted(latitudes, maxlat)
1064
1065    if minlon == None:
1066        lmin=0
1067    else:   
1068        lmin = searchsorted(longitudes, minlon)
1069   
1070    if maxlon == None:
1071        lmax = len(longitudes)
1072    else:
1073        lmax = searchsorted(longitudes, maxlon)           
1074       
1075    latitudes = latitudes[kmin:kmax]
1076    longitudes = longitudes[lmin:lmax]
1077    times = times[jmin:jmax]
1078   
1079    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1080    xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax]
1081    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]   
1082
1083    number_of_latitudes = latitudes.shape[0]
1084    number_of_longitudes = longitudes.shape[0]
1085
1086
1087    #print times
1088    #print latitudes
1089    #print longitudes
1090
1091    #print 'MIN', min(min(min(amplitudes)))
1092    #print 'MAX', max(max(max(amplitudes)))   
1093
1094    #print number_of_latitudes, number_of_longitudes
1095    number_of_points = number_of_latitudes*number_of_longitudes
1096    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
1097   
1098    #print file_h.dimensions.keys()
1099    #print file_h.variables.keys()   
1100
1101    if verbose: print 'Store to SWW file %s' %swwname
1102    # NetCDF file definition
1103    outfile = NetCDFFile(swwname, 'w')
1104       
1105    #Create new file
1106    outfile.institution = 'Geoscience Australia'
1107    outfile.description = 'Converted from Ferret files: %s, %s, %s'\
1108                          %(basefilename + '_ha.nc',
1109                            basefilename + '_ua.nc',
1110                            basefilename + '_va.nc')
1111
1112
1113    #For sww compatibility
1114    outfile.smoothing = 'Yes' 
1115    outfile.order = 1
1116           
1117    #Start time in seconds since the epoch (midnight 1/1/1970)
1118    outfile.starttime = times[0]
1119   
1120    # dimension definitions
1121    outfile.createDimension('number_of_volumes', number_of_volumes)
1122
1123    outfile.createDimension('number_of_vertices', 3)
1124    outfile.createDimension('number_of_points', number_of_points)
1125                           
1126               
1127    #outfile.createDimension('number_of_timesteps', len(times))
1128    outfile.createDimension('number_of_timesteps', len(times))
1129
1130    # variable definitions
1131    outfile.createVariable('x', precision, ('number_of_points',))
1132    outfile.createVariable('y', precision, ('number_of_points',))
1133    outfile.createVariable('elevation', precision, ('number_of_points',))
1134           
1135    #FIXME: Backwards compatibility
1136    outfile.createVariable('z', precision, ('number_of_points',))
1137    #################################
1138       
1139    outfile.createVariable('volumes', Int, ('number_of_volumes',
1140                                            'number_of_vertices'))
1141   
1142    outfile.createVariable('time', precision,
1143                           ('number_of_timesteps',))
1144       
1145    outfile.createVariable('stage', precision,
1146                           ('number_of_timesteps',
1147                            'number_of_points'))
1148
1149    outfile.createVariable('xmomentum', precision,
1150                           ('number_of_timesteps',
1151                            'number_of_points'))
1152   
1153    outfile.createVariable('ymomentum', precision,
1154                           ('number_of_timesteps',
1155                            'number_of_points'))
1156
1157
1158    #Store
1159    from coordinate_transforms.redfearn import redfearn
1160    x = zeros(number_of_points, Float)  #Easting
1161    y = zeros(number_of_points, Float)  #Northing
1162    #volumes = zeros(number_of_volumes, Int)
1163    i = 0
1164
1165    #Check zone boundaries
1166    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
1167
1168    vertices = {}
1169    for l, lon in enumerate(longitudes):   
1170        for k, lat in enumerate(latitudes):
1171            vertices[l,k] = i
1172           
1173            zone, easting, northing = redfearn(lat,lon)
1174
1175            msg = 'Zone boundary crossed at longitude =', lon
1176            assert zone == refzone, msg
1177            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
1178            x[i] = easting
1179            y[i] = northing
1180            i += 1
1181
1182
1183    #Construct 2 triangles per 'rectangular' element
1184    volumes = []
1185    for l in range(number_of_longitudes-1):
1186        for k in range(number_of_latitudes-1):
1187            v1 = vertices[l,k+1]
1188            v2 = vertices[l,k]           
1189            v3 = vertices[l+1,k+1]           
1190            v4 = vertices[l+1,k]           
1191
1192            volumes.append([v1,v2,v3]) #Upper element
1193            volumes.append([v4,v3,v2]) #Lower element           
1194
1195    volumes = array(volumes)       
1196
1197    origin = (min(x), min(y))
1198    outfile.minx = origin[0]
1199    outfile.miny = origin[1]   
1200
1201    #print x - origin[0]
1202    #print y - origin[1]
1203
1204    #print len(x)
1205    #print number_of_points
1206   
1207    outfile.variables['x'][:] = x - origin[0]
1208    outfile.variables['y'][:] = y - origin[1]
1209    outfile.variables['z'][:] = 0.0
1210    outfile.variables['elevation'][:] = 0.0   
1211    outfile.variables['volumes'][:] = volumes
1212    outfile.variables['time'][:] = times   
1213
1214    #Time stepping
1215    stage = outfile.variables['stage']
1216    xmomentum = outfile.variables['xmomentum']
1217    ymomentum = outfile.variables['ymomentum']       
1218
1219    for j in range(len(times)):
1220        i = 0
1221        for l in range(number_of_longitudes):   
1222            for k in range(number_of_latitudes):
1223                h = amplitudes[j,k,l]/100 + mean_stage
1224                stage[j,i] = h
1225                xmomentum[j,i] = xspeed[j,k,l]/100*h
1226                ymomentum[j,i] = yspeed[j,k,l]/100*h               
1227                i += 1
1228       
1229    outfile.close()               
1230   
1231   
1232
1233
1234
1235#OBSOLETE STUFF
1236#Native checkpoint format.
1237#Information needed to recreate a state is preserved
1238#FIXME: Rethink and maybe use netcdf format
1239def cpt_variable_writer(filename, t, v0, v1, v2):
1240    """Store all conserved quantities to file
1241    """
1242
1243    M, N = v0.shape
1244       
1245    FN = create_filename(filename, 'cpt', M, t)
1246    #print 'Writing to %s' %FN
1247   
1248    fid = open(FN, 'w')
1249    for i in range(M):
1250        for j in range(N):       
1251            fid.write('%.16e ' %v0[i,j])           
1252        for j in range(N):
1253            fid.write('%.16e ' %v1[i,j])                       
1254        for j in range(N):
1255            fid.write('%.16e ' %v2[i,j])                       
1256           
1257        fid.write('\n')
1258    fid.close()
1259
1260
1261def cpt_variable_reader(filename, t, v0, v1, v2):
1262    """Store all conserved quantities to file
1263    """
1264
1265    M, N = v0.shape
1266       
1267    FN = create_filename(filename, 'cpt', M, t)
1268    #print 'Reading from %s' %FN
1269   
1270    fid = open(FN)
1271
1272
1273    for i in range(M):
1274        values = fid.readline().split() #Get one line
1275
1276        for j in range(N):
1277            v0[i,j] = float(values[j])
1278            v1[i,j] = float(values[3+j])
1279            v2[i,j] = float(values[6+j])
1280   
1281    fid.close()
1282
1283def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1284    """Writes x,y,z,z,z coordinates of triangles constituting the bed
1285    elevation.
1286    Not in use pt
1287    """
1288
1289    M, N = v0.shape
1290
1291    print X0
1292    import sys; sys.exit() 
1293    FN = create_filename(filename, 'cpt', M)
1294    print 'Writing to %s' %FN
1295   
1296    fid = open(FN, 'w')
1297    for i in range(M):
1298        for j in range(2):       
1299            fid.write('%.16e ' %X0[i,j])   #x, y
1300        for j in range(N):                   
1301            fid.write('%.16e ' %v0[i,j])       #z,z,z,
1302           
1303        for j in range(2):                   
1304            fid.write('%.16e ' %X1[i,j])   #x, y
1305        for j in range(N):                               
1306            fid.write('%.16e ' %v1[i,j]) 
1307       
1308        for j in range(2):                               
1309            fid.write('%.16e ' %X2[i,j])   #x, y
1310        for j in range(N):                                           
1311            fid.write('%.16e ' %v2[i,j]) 
1312           
1313        fid.write('\n')
1314    fid.close()
1315
1316
1317
1318#Function for storing out to e.g. visualisation
1319#FIXME: Do we want this?
1320#FIXME: Not done yet for this version
1321def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
1322    """Writes x,y,z coordinates of triangles constituting the bed elevation.
1323    """
1324
1325    M, N = v0.shape
1326       
1327    FN = create_filename(filename, 'dat', M)
1328    #print 'Writing to %s' %FN
1329   
1330    fid = open(FN, 'w')
1331    for i in range(M):
1332        for j in range(2):       
1333            fid.write('%f ' %X0[i,j])   #x, y
1334        fid.write('%f ' %v0[i,0])       #z
1335           
1336        for j in range(2):                   
1337            fid.write('%f ' %X1[i,j])   #x, y
1338        fid.write('%f ' %v1[i,0])       #z
1339       
1340        for j in range(2):                               
1341            fid.write('%f ' %X2[i,j])   #x, y
1342        fid.write('%f ' %v2[i,0])       #z
1343           
1344        fid.write('\n')
1345    fid.close()
1346
1347
1348
1349def dat_variable_writer(filename, t, v0, v1, v2):
1350    """Store water height to file
1351    """
1352
1353    M, N = v0.shape
1354       
1355    FN = create_filename(filename, 'dat', M, t)
1356    #print 'Writing to %s' %FN
1357   
1358    fid = open(FN, 'w')
1359    for i in range(M):
1360        fid.write('%.4f ' %v0[i,0])
1361        fid.write('%.4f ' %v1[i,0])
1362        fid.write('%.4f ' %v2[i,0])       
1363
1364        fid.write('\n')
1365    fid.close()
1366
Note: See TracBrowser for help on using the repository browser.