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

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

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