source: inundation/pyvolution/data_manager.py @ 1875

Last change on this file since 1875 was 1875, checked in by ole, 18 years ago

Bounding box for sww2dem

File size: 93.7 KB
RevLine 
[1665]1"""datamanager.py - input output for AnuGA
[1502]2
3
[1665]4This module takes care of reading and writing datafiles such as topograhies,
5model output, etc
[1502]6
7
[1665]8Formats used within AnuGA:
[1502]9
[1680]10.sww: Netcdf format for storing model output f(t,x,y)
11.tms: Netcdf format for storing time series f(t)
[1502]12
13.xya: ASCII format for storing arbitrary points and associated attributes
14.pts: NetCDF format for storing arbitrary points and associated attributes
15
[1866]16.asc: ASCII format of regular DEMs as output from ArcView
[1502]17.prj: Associated ArcView file giving more meta data for asc format
[1866]18.ers: ERMapper header format of regular DEMs for ArcView
[1502]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.geo: Houdinis ascii geometry format (?)
27
28
29A typical dataflow can be described as follows
30
31Manually created files:
32ASC, PRJ:     Digital elevation models (gridded)
33TSH:          Triangular meshes (e.g. created from pmesh)
34NC            Model outputs for use as boundary conditions (e.g from MOST)
35
36
37AUTOMATICALLY CREATED FILES:
38
39ASC, PRJ  ->  DEM  ->  PTS: Conversion of DEM's to native pts file
40
41NC -> SWW: Conversion of MOST bundary files to boundary sww
42
43PTS + TSH -> TSH with elevation: Least squares fit
44
45TSH -> SWW:  Conversion of TSH to sww viewable using Swollen
46
47TSH + Boundary SWW -> SWW: Simluation using pyvolution
48
49"""
50
51from Numeric import concatenate
52
[1661]53from coordinate_transforms.geo_reference import Geo_reference
[1502]54
55def make_filename(s):
56    """Transform argument string into a suitable filename
57    """
58
59    s = s.strip()
60    s = s.replace(' ', '_')
61    s = s.replace('(', '')
62    s = s.replace(')', '')
63    s = s.replace('__', '_')
64
65    return s
66
67
68def check_dir(path, verbose=None):
69    """Check that specified path exists.
70    If path does not exist it will be created if possible
71
72    USAGE:
73       checkdir(path, verbose):
74
75    ARGUMENTS:
76        path -- Directory
77        verbose -- Flag verbose output (default: None)
78
79    RETURN VALUE:
80        Verified path including trailing separator
81
82    """
83
84    import os, sys
85    import os.path
86
87    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
88        unix = 0
89    else:
90        unix = 1
91
92
93    if path[-1] != os.sep:
94        path = path + os.sep  # Add separator for directories
95
96    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
97    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
98        try:
99            exitcode=os.mkdir(path)
100
101            # Change access rights if possible
102            #
103            if unix:
104                exitcode=os.system('chmod 775 '+path)
105            else:
106                pass  # FIXME: What about acces rights under Windows?
107
108            if verbose: print 'MESSAGE: Directory', path, 'created.'
109
110        except:
111            print 'WARNING: Directory', path, 'could not be created.'
112            if unix:
113                path = '/tmp/'
114            else:
115                path = 'C:'
116
117            print 'Using directory %s instead' %path
118
119    return(path)
120
121
122
123def del_dir(path):
124    """Recursively delete directory path and all its contents
125    """
126
127    import os
128
129    if os.path.isdir(path):
130        for file in os.listdir(path):
131            X = os.path.join(path, file)
132
133
134            if os.path.isdir(X) and not os.path.islink(X):
135                del_dir(X)
136            else:
137                try:
138                    os.remove(X)
139                except:
140                    print "Could not remove file %s" %X
141
142        os.rmdir(path)
143
144
145
146def create_filename(datadir, filename, format, size=None, time=None):
147
148    import os
149    #from config import data_dir
150
151    FN = check_dir(datadir) + filename
152
153    if size is not None:
154        FN += '_size%d' %size
155
156    if time is not None:
157        FN += '_time%.2f' %time
158
159    FN += '.' + format
160    return FN
161
162
163def get_files(datadir, filename, format, size):
164    """Get all file (names) with gven name, size and format
165    """
166
167    import glob
168
169    import os
170    #from config import data_dir
171
172    dir = check_dir(datadir)
173
174    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
175    return glob.glob(pattern)
176
177
178
179#Generic class for storing output to e.g. visualisation or checkpointing
180class Data_format:
181    """Generic interface to data formats
182    """
183
184
185    def __init__(self, domain, extension, mode = 'w'):
186        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
187                                        '''   'w' (write)'''+\
188                                        '''   'r' (read)''' +\
189                                        '''   'a' (append)'''
190
191        #Create filename
192        #self.filename = create_filename(domain.get_datadir(),
193        #domain.get_name(), extension, len(domain))
194
195
196        self.filename = create_filename(domain.get_datadir(),
197                        domain.get_name(), extension)
198
199        #print 'F', self.filename
200        self.timestep = 0
201        self.number_of_volumes = len(domain)
202        self.domain = domain
203
204
205        #FIXME: Should we have a general set_precision function?
206
207
208
209#Class for storing output to e.g. visualisation
210class Data_format_sww(Data_format):
[1665]211    """Interface to native NetCDF format (.sww) for storing model output
212
213    There are two kinds of data
214
215    1: Constant data: Vertex coordinates and field values. Stored once
216    2: Variable data: Conserved quantities. Stored once per timestep.
217
218    All data is assumed to reside at vertex locations.
[1502]219    """
220
221
222    def __init__(self, domain, mode = 'w',\
223                 max_size = 2000000000,
[1854]224                 recursion = False):
[1502]225        from Scientific.IO.NetCDF import NetCDFFile
226        from Numeric import Int, Float, Float32
227
228        self.precision = Float32 #Use single precision
229        if hasattr(domain, 'max_size'):
230            self.max_size = domain.max_size #file size max is 2Gig
231        else:
232            self.max_size = max_size
233        self.recursion = recursion
234        self.mode = mode
235
236        Data_format.__init__(self, domain, 'sww', mode)
237
238
239        # NetCDF file definition
240        fid = NetCDFFile(self.filename, mode)
241
242        if mode == 'w':
243
244            #Create new file
245            fid.institution = 'Geoscience Australia'
246            fid.description = 'Output from pyvolution suitable for plotting'
247
248            if domain.smooth:
249                fid.smoothing = 'Yes'
250            else:
251                fid.smoothing = 'No'
252
253            fid.order = domain.default_order
254           
255            if hasattr(domain, 'texture'):                 
256                fid.texture = domain.texture
[1672]257            #else:                         
258            #    fid.texture = 'None'       
[1502]259
260            #Reference point
261            #Start time in seconds since the epoch (midnight 1/1/1970)
262            #FIXME: Use Georef
263            fid.starttime = domain.starttime
264
265            # dimension definitions
266            fid.createDimension('number_of_volumes', self.number_of_volumes)
267            fid.createDimension('number_of_vertices', 3)
268
269            if domain.smooth is True:
270                fid.createDimension('number_of_points', len(domain.vertexlist))
271            else:
272                fid.createDimension('number_of_points', 3*self.number_of_volumes)
273
274            fid.createDimension('number_of_timesteps', None) #extensible
275
276            # variable definitions
277            fid.createVariable('x', self.precision, ('number_of_points',))
278            fid.createVariable('y', self.precision, ('number_of_points',))
279            fid.createVariable('elevation', self.precision, ('number_of_points',))
280            if domain.geo_reference is not None:
281                domain.geo_reference.write_NetCDF(fid)
282               
283            #FIXME: Backwards compatibility
284            fid.createVariable('z', self.precision, ('number_of_points',))
285            #################################
286
287            fid.createVariable('volumes', Int, ('number_of_volumes',
288                                                'number_of_vertices'))
289
290            fid.createVariable('time', self.precision,
291                               ('number_of_timesteps',))
292
293            fid.createVariable('stage', self.precision,
294                               ('number_of_timesteps',
295                                'number_of_points'))
296
297            fid.createVariable('xmomentum', self.precision,
298                               ('number_of_timesteps',
299                                'number_of_points'))
300
301            fid.createVariable('ymomentum', self.precision,
302                               ('number_of_timesteps',
303                                'number_of_points'))
304
305        #Close
306        fid.close()
307
308
309    def store_connectivity(self):
310        """Specialisation of store_connectivity for net CDF format
311
312        Writes x,y,z coordinates of triangles constituting
313        the bed elevation.
314        """
315
316        from Scientific.IO.NetCDF import NetCDFFile
317
318        from Numeric import concatenate, Int
319
320        domain = self.domain
321
322        #Get NetCDF
323        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
324
325        # Get the variables
326        x = fid.variables['x']
327        y = fid.variables['y']
328        z = fid.variables['elevation']
329
330        volumes = fid.variables['volumes']
331
332        # Get X, Y and bed elevation Z
333        Q = domain.quantities['elevation']
334        X,Y,Z,V = Q.get_vertex_values(xy=True,
335                                      precision = self.precision)
336
337
338
339        x[:] = X.astype(self.precision)
340        y[:] = Y.astype(self.precision)
341        z[:] = Z.astype(self.precision)
342
343        #FIXME: Backwards compatibility
344        z = fid.variables['z']
345        z[:] = Z.astype(self.precision)
346        ################################
347
348        volumes[:] = V.astype(volumes.typecode())
349
350        #Close
351        fid.close()
352
353    def store_timestep(self, names):
354        """Store time and named quantities to file
355        """
356        from Scientific.IO.NetCDF import NetCDFFile
357        import types
358        from time import sleep
359        from os import stat
360
[1740]361        minimum_allowed_depth = 0.001
362        #minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
363        from Numeric import choose
[1502]364
365        #Get NetCDF
366        retries = 0
367        file_open = False
368        while not file_open and retries < 10:
369            try:
370                fid = NetCDFFile(self.filename, 'a') #Open existing file
371            except IOError:
372                #This could happen if someone was reading the file.
373                #In that case, wait a while and try again
374                msg = 'Warning (store_timestep): File %s could not be opened'\
375                      %self.filename
[1690]376                msg += ' - trying step %s again' %self.domain.time
[1502]377                print msg
378                retries += 1
379                sleep(1)
380            else:
381                file_open = True
382
383        if not file_open:
384            msg = 'File %s could not be opened for append' %self.filename
385            raise msg
386
387
388
389        #Check to see if the file is already too big:
390        time = fid.variables['time']
391        i = len(time)+1
392        file_size = stat(self.filename)[6]
393        file_size_increase =  file_size/i
394        if file_size + file_size_increase > self.max_size*(2**self.recursion):
395            #in order to get the file name and start time correct,
396            #I change the domian.filename and domain.starttime.
397            #This is the only way to do this without changing
398            #other modules (I think).
399
400            #write a filename addon that won't break swollens reader
401            #(10.sww is bad)
402            filename_ext = '_time_%s'%self.domain.time
403            filename_ext = filename_ext.replace('.', '_')
404            #remember the old filename, then give domain a
405            #name with the extension
406            old_domain_filename = self.domain.filename
407            if not self.recursion:
408                self.domain.filename = self.domain.filename+filename_ext
409
410            #change the domain starttime to the current time
411            old_domain_starttime = self.domain.starttime
412            self.domain.starttime = self.domain.time
413
414            #build a new data_structure.
415            next_data_structure=\
416                Data_format_sww(self.domain, mode=self.mode,\
417                                max_size = self.max_size,\
418                                recursion = self.recursion+1)
419            if not self.recursion:
420                print '    file_size = %s'%file_size
421                print '    saving file to %s'%next_data_structure.filename
422            #set up the new data_structure
423            self.domain.writer = next_data_structure
424
425            #FIXME - could be cleaner to use domain.store_timestep etc.
426            next_data_structure.store_connectivity()
427            next_data_structure.store_timestep(names)
428            fid.sync()
429            fid.close()
430
431            #restore the old starttime and filename
432            self.domain.starttime = old_domain_starttime
433            self.domain.filename = old_domain_filename
434        else:
435            self.recursion = False
436            domain = self.domain
437
438            # Get the variables
439            time = fid.variables['time']
440            stage = fid.variables['stage']
441            xmomentum = fid.variables['xmomentum']
442            ymomentum = fid.variables['ymomentum']
443            i = len(time)
444
445            #Store time
446            time[i] = self.domain.time
447
448
449            if type(names) not in [types.ListType, types.TupleType]:
450                names = [names]
451
452            for name in names:
453                # Get quantity
454                Q = domain.quantities[name]
[1740]455                A,V = Q.get_vertex_values(xy = False,
456                                          precision = self.precision)
[1502]457
458                #FIXME: Make this general (see below)
459                if name == 'stage':
[1740]460                    z = fid.variables['elevation']
461                    #print z[:]
462                    #print A-z[:]
463                    A = choose( A-z[:] >= minimum_allowed_depth, (z[:], A))
[1502]464                    stage[i,:] = A.astype(self.precision)
465                elif name == 'xmomentum':
466                    xmomentum[i,:] = A.astype(self.precision)
467                elif name == 'ymomentum':
468                    ymomentum[i,:] = A.astype(self.precision)
469
470                #As in....
471                #eval( name + '[i,:] = A.astype(self.precision)' )
472                #FIXME: But we need a UNIT test for that before refactoring
473
474
475
476            #Flush and close
477            fid.sync()
478            fid.close()
479
480
481
482#Class for handling checkpoints data
483class Data_format_cpt(Data_format):
484    """Interface to native NetCDF format (.cpt)
485    """
486
487
488    def __init__(self, domain, mode = 'w'):
489        from Scientific.IO.NetCDF import NetCDFFile
490        from Numeric import Int, Float, Float
491
492        self.precision = Float #Use full precision
493
494        Data_format.__init__(self, domain, 'sww', mode)
495
496
497        # NetCDF file definition
498        fid = NetCDFFile(self.filename, mode)
499
500        if mode == 'w':
501            #Create new file
502            fid.institution = 'Geoscience Australia'
503            fid.description = 'Checkpoint data'
504            #fid.smooth = domain.smooth
505            fid.order = domain.default_order
506
507            # dimension definitions
508            fid.createDimension('number_of_volumes', self.number_of_volumes)
509            fid.createDimension('number_of_vertices', 3)
510
511            #Store info at all vertices (no smoothing)
512            fid.createDimension('number_of_points', 3*self.number_of_volumes)
513            fid.createDimension('number_of_timesteps', None) #extensible
514
515            # variable definitions
516
517            #Mesh
518            fid.createVariable('x', self.precision, ('number_of_points',))
519            fid.createVariable('y', self.precision, ('number_of_points',))
520
521
522            fid.createVariable('volumes', Int, ('number_of_volumes',
523                                                'number_of_vertices'))
524
525            fid.createVariable('time', self.precision,
526                               ('number_of_timesteps',))
527
528            #Allocate space for all quantities
529            for name in domain.quantities.keys():
530                fid.createVariable(name, self.precision,
531                                   ('number_of_timesteps',
532                                    'number_of_points'))
533
534        #Close
535        fid.close()
536
537
538    def store_checkpoint(self):
539        """
540        Write x,y coordinates of triangles.
541        Write connectivity (
542        constituting
543        the bed elevation.
544        """
545
546        from Scientific.IO.NetCDF import NetCDFFile
547
548        from Numeric import concatenate
549
550        domain = self.domain
551
552        #Get NetCDF
553        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
554
555        # Get the variables
556        x = fid.variables['x']
557        y = fid.variables['y']
558
559        volumes = fid.variables['volumes']
560
561        # Get X, Y and bed elevation Z
562        Q = domain.quantities['elevation']
563        X,Y,Z,V = Q.get_vertex_values(xy=True,
564                                      precision = self.precision)
565
566
567
568        x[:] = X.astype(self.precision)
569        y[:] = Y.astype(self.precision)
570        z[:] = Z.astype(self.precision)
571
572        volumes[:] = V
573
574        #Close
575        fid.close()
576
577
578    def store_timestep(self, name):
579        """Store time and named quantity to file
580        """
581        from Scientific.IO.NetCDF import NetCDFFile
582        from time import sleep
583
584        #Get NetCDF
585        retries = 0
586        file_open = False
587        while not file_open and retries < 10:
588            try:
589                fid = NetCDFFile(self.filename, 'a') #Open existing file
590            except IOError:
591                #This could happen if someone was reading the file.
592                #In that case, wait a while and try again
593                msg = 'Warning (store_timestep): File %s could not be opened'\
594                          %self.filename
595                msg += ' - trying again'
596                print msg
597                retries += 1
598                sleep(1)
599            else:
600                file_open = True
601
602            if not file_open:
603                msg = 'File %s could not be opened for append' %self.filename
604            raise msg
605
606
607        domain = self.domain
608
609        # Get the variables
610        time = fid.variables['time']
611        stage = fid.variables['stage']
612        i = len(time)
613
614        #Store stage
615        time[i] = self.domain.time
616
617        # Get quantity
618        Q = domain.quantities[name]
619        A,V = Q.get_vertex_values(xy=False,
620                                  precision = self.precision)
621
622        stage[i,:] = A.astype(self.precision)
623
624        #Flush and close
625        fid.sync()
626        fid.close()
627
628
629
630
631
632#Function for storing xya output
633#FIXME Not done yet for this version
634class Data_format_xya(Data_format):
635    """Generic interface to data formats
636    """
637
638
639    def __init__(self, domain, mode = 'w'):
640        from Scientific.IO.NetCDF import NetCDFFile
641        from Numeric import Int, Float, Float32
642
643        self.precision = Float32 #Use single precision
644
645        Data_format.__init__(self, domain, 'xya', mode)
646
647
648
649    #FIXME -This is the old xya format
650    def store_all(self):
651        """Specialisation of store all for xya format
652
653        Writes x,y,z coordinates of triangles constituting
654        the bed elevation.
655        """
656
657        from Numeric import concatenate
658
659        domain = self.domain
660
661        fd = open(self.filename, 'w')
662
663
664        if domain.smooth is True:
665            number_of_points =  len(domain.vertexlist)
666        else:
667            number_of_points = 3*self.number_of_volumes
668
669        numVertAttrib = 3 #Three attributes is what is assumed by the xya format
670
671        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
672                 " # <vertex #> <x> <y> [attributes]" + "\n")
673
674
675        # Get X, Y, bed elevation and friction (index=0,1)
676        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
677                                           indices = (0,1), precision = self.precision)
678
679        bed_eles = A[:,0]
680        fricts = A[:,1]
681
682        # Get stage (index=0)
683        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
684                                       indices = (0,), precision = self.precision)
685
686        stages = B[:,0]
687
688        #<vertex #> <x> <y> [attributes]
689        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
690                                               stages, fricts):
691
692            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
693            fd.write(s)
694
695        #close
696        fd.close()
697
698
699    def store_timestep(self, t, V0, V1, V2):
700        """Store time, water heights (and momentums) to file
701        """
702        pass
703
704
705#Auxiliary
706def write_obj(filename,x,y,z):
707    """Store x,y,z vectors into filename (obj format)
708       Vectors are assumed to have dimension (M,3) where
709       M corresponds to the number elements.
710       triangles are assumed to be disconnected
711
712       The three numbers in each vector correspond to three vertices,
713
714       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
715
716    """
717    #print 'Writing obj to %s' % filename
718
719    import os.path
720
721    root, ext = os.path.splitext(filename)
722    if ext == '.obj':
723        FN = filename
724    else:
725        FN = filename + '.obj'
726
727
728    outfile = open(FN, 'wb')
729    outfile.write("# Triangulation as an obj file\n")
730
731    M, N = x.shape
732    assert N==3  #Assuming three vertices per element
733
734    for i in range(M):
735        for j in range(N):
736            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
737
738    for i in range(M):
739        base = i*N
740        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
741
742    outfile.close()
743
744
[1671]745#########################################################
746#Conversion routines
747########################################################
[1502]748
749def sww2obj(basefilename, size):
750    """Convert netcdf based data output to obj
751    """
752    from Scientific.IO.NetCDF import NetCDFFile
753
754    from Numeric import Float, zeros
755
756    #Get NetCDF
757    FN = create_filename('.', basefilename, 'sww', size)
758    print 'Reading from ', FN
759    fid = NetCDFFile(FN, 'r')  #Open existing file for read
760
761
762    # Get the variables
763    x = fid.variables['x']
764    y = fid.variables['y']
765    z = fid.variables['elevation']
766    time = fid.variables['time']
767    stage = fid.variables['stage']
768
769    M = size  #Number of lines
770    xx = zeros((M,3), Float)
771    yy = zeros((M,3), Float)
772    zz = zeros((M,3), Float)
773
774    for i in range(M):
775        for j in range(3):
776            xx[i,j] = x[i+j*M]
777            yy[i,j] = y[i+j*M]
778            zz[i,j] = z[i+j*M]
779
780    #Write obj for bathymetry
781    FN = create_filename('.', basefilename, 'obj', size)
782    write_obj(FN,xx,yy,zz)
783
784
785    #Now read all the data with variable information, combine with
786    #x,y info and store as obj
787
788    for k in range(len(time)):
789        t = time[k]
790        print 'Processing timestep %f' %t
791
792        for i in range(M):
793            for j in range(3):
794                zz[i,j] = stage[k,i+j*M]
795
796
797        #Write obj for variable data
798        #FN = create_filename(basefilename, 'obj', size, time=t)
799        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
800        write_obj(FN,xx,yy,zz)
801
802
803def dat2obj(basefilename):
804    """Convert line based data output to obj
805    FIXME: Obsolete?
806    """
807
808    import glob, os
809    from config import data_dir
810
811
812    #Get bathymetry and x,y's
813    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
814
815    from Numeric import zeros, Float
816
817    M = len(lines)  #Number of lines
818    x = zeros((M,3), Float)
819    y = zeros((M,3), Float)
820    z = zeros((M,3), Float)
821
822    ##i = 0
823    for i, line in enumerate(lines):
824        tokens = line.split()
825        values = map(float,tokens)
826
827        for j in range(3):
828            x[i,j] = values[j*3]
829            y[i,j] = values[j*3+1]
830            z[i,j] = values[j*3+2]
831
832        ##i += 1
833
834
835    #Write obj for bathymetry
836    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
837
838
839    #Now read all the data files with variable information, combine with
840    #x,y info
841    #and store as obj
842
843    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
844
845    for filename in files:
846        print 'Processing %s' % filename
847
848        lines = open(data_dir+os.sep+filename,'r').readlines()
849        assert len(lines) == M
850        root, ext = os.path.splitext(filename)
851
852        #Get time from filename
853        i0 = filename.find('_time=')
854        if i0 == -1:
855            #Skip bathymetry file
856            continue
857
858        i0 += 6  #Position where time starts
859        i1 = filename.find('.dat')
860
861        if i1 > i0:
862            t = float(filename[i0:i1])
863        else:
864            raise 'Hmmmm'
865
866
867
868        ##i = 0
869        for i, line in enumerate(lines):
870            tokens = line.split()
871            values = map(float,tokens)
872
873            for j in range(3):
874                z[i,j] = values[j]
875
876            ##i += 1
877
878        #Write obj for variable data
879        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
880
881
882def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
883    """Read netcdf filename1, pick timesteps first:step:last and save to
884    nettcdf file filename2
885    """
886    from Scientific.IO.NetCDF import NetCDFFile
887
888    #Get NetCDF
889    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
890    outfile = NetCDFFile(filename2, 'w')  #Open new file
891
892
893    #Copy dimensions
894    for d in infile.dimensions:
895        outfile.createDimension(d, infile.dimensions[d])
896
897    for name in infile.variables:
898        var = infile.variables[name]
899        outfile.createVariable(name, var.typecode(), var.dimensions)
900
901
902    #Copy the static variables
903    for name in infile.variables:
904        if name == 'time' or name == 'stage':
905            pass
906        else:
907            #Copy
908            outfile.variables[name][:] = infile.variables[name][:]
909
910    #Copy selected timesteps
911    time = infile.variables['time']
912    stage = infile.variables['stage']
913
914    newtime = outfile.variables['time']
915    newstage = outfile.variables['stage']
916
917    if last is None:
918        last = len(time)
919
920    selection = range(first, last, step)
921    for i, j in enumerate(selection):
922        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
923        newtime[i] = time[j]
924        newstage[i,:] = stage[j,:]
925
926    #Close
927    infile.close()
928    outfile.close()
929
930
931#Get data objects
932def get_dataobject(domain, mode='w'):
933    """Return instance of class of given format using filename
934    """
935
936    cls = eval('Data_format_%s' %domain.format)
937    return cls(domain, mode)
938
[1661]939def xya2pts(basename_in, basename_out=None, verbose=False,
940            #easting_min=None, easting_max=None,
941            #northing_min=None, northing_max=None,
[1740]942            stride = 1,           
[1680]943            attribute_name = 'elevation',
944            z_func = None):
[1661]945    """Read points data from ascii (.xya)
946
947    Example:
948
949              x(m)        y(m)        z(m)
950         0.00000e+00  0.00000e+00  1.3535000e-01
951         0.00000e+00  1.40000e-02  1.3535000e-01
952
953   
954
955    Convert to NetCDF pts format which is
956
957    points:  (Nx2) Float array
958    elevation: N Float array
959
960    Only lines that contain three numeric values are processed
[1680]961
962    If z_func is specified, it will be applied to the third field
[1661]963    """
964
965    import os 
[1753]966    #from Scientific.IO.NetCDF import NetCDFFile
[1661]967    from Numeric import Float, arrayrange, concatenate
968
969    root, ext = os.path.splitext(basename_in)
970
971    if ext == '': ext = '.xya'
972
973    #Get NetCDF
974    infile = open(root + ext, 'r')  #Open existing xya file for read
975
976    if verbose: print 'Reading xya points from %s' %(root + ext)
977
978    points = []
979    attribute = []
980    for i, line in enumerate(infile.readlines()):
[1740]981
982        if i % stride != 0: continue
983       
[1661]984        fields = line.split()
985
986        try:
987            assert len(fields) == 3
988        except:   
989            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) 
990
991        try:
992            x = float( fields[0] )
993            y = float( fields[1] )
994            z = float( fields[2] )           
995        except:
996            continue
997
998        points.append( [x, y] )
999
[1680]1000        if callable(z_func):
1001            attribute.append(z_func(z))           
1002        else:   
1003            attribute.append(z)
[1661]1004
[1680]1005
[1661]1006    #Get output file
1007    if basename_out == None:
1008        ptsname = root + '.pts'
1009    else:
1010        ptsname = basename_out + '.pts'
1011
1012    if verbose: print 'Store to NetCDF file %s' %ptsname
[1753]1013    write_ptsfile(ptsname, points, attribute, attribute_name)
1014
1015
1016def write_ptsfile(ptsname, points, attribute, attribute_name = None):
1017    """Write points and associated attribute to pts (NetCDF) format
1018    """
1019
1020    from Numeric import Float
1021   
1022    if attribute_name is None:
1023        attribute_name = 'attribute'       
1024   
1025
1026    from Scientific.IO.NetCDF import NetCDFFile
1027   
[1661]1028    # NetCDF file definition
1029    outfile = NetCDFFile(ptsname, 'w')
1030
[1753]1031
[1661]1032    #Create new file
1033    outfile.institution = 'Geoscience Australia'
1034    outfile.description = 'NetCDF pts format for compact and '\
1035                          'portable storage of spatial point data'
1036   
1037
1038    #Georeferencing
1039    from coordinate_transforms.geo_reference import Geo_reference
1040    Geo_reference().write_NetCDF(outfile)
1041   
1042
1043    outfile.createDimension('number_of_points', len(points))
1044    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1045
1046    # variable definitions
1047    outfile.createVariable('points', Float, ('number_of_points',
1048                                             'number_of_dimensions'))
1049    outfile.createVariable(attribute_name, Float, ('number_of_points',))
1050
1051    # Get handles to the variables
1052    nc_points = outfile.variables['points']
1053    nc_attribute = outfile.variables[attribute_name]
1054
1055    #Store data
1056    nc_points[:, :] = points
1057    nc_attribute[:] = attribute
1058
1059    outfile.close()
1060   
[1753]1061
[1577]1062def dem2pts(basename_in, basename_out=None, verbose=False,
1063            easting_min=None, easting_max=None,
1064            northing_min=None, northing_max=None):
[1502]1065    """Read Digitial Elevation model from the following NetCDF format (.dem)
1066
1067    Example:
1068
1069    ncols         3121
1070    nrows         1800
1071    xllcorner     722000
1072    yllcorner     5893000
1073    cellsize      25
1074    NODATA_value  -9999
1075    138.3698 137.4194 136.5062 135.5558 ..........
1076
1077    Convert to NetCDF pts format which is
1078
1079    points:  (Nx2) Float array
1080    elevation: N Float array
1081    """
1082
[1753]1083    #FIXME: Can this be written feasibly using write_pts?
1084
[1502]1085    import os
1086    from Scientific.IO.NetCDF import NetCDFFile
[1740]1087    from Numeric import Float, zeros, reshape
[1502]1088
1089    root = basename_in
1090
1091    #Get NetCDF
1092    infile = NetCDFFile(root + '.dem', 'r')  #Open existing netcdf file for read
1093
1094    if verbose: print 'Reading DEM from %s' %(root + '.dem')
1095
1096    ncols = infile.ncols[0]
1097    nrows = infile.nrows[0]
1098    xllcorner = infile.xllcorner[0]  #Easting of lower left corner
1099    yllcorner = infile.yllcorner[0]  #Northing of lower left corner
1100    cellsize = infile.cellsize[0]
1101    NODATA_value = infile.NODATA_value[0]
1102    dem_elevation = infile.variables['elevation']
1103
1104    zone = infile.zone[0]
1105    false_easting = infile.false_easting[0]
1106    false_northing = infile.false_northing[0]
1107
1108    #Text strings
1109    projection = infile.projection
1110    datum = infile.datum
1111    units = infile.units
1112
1113
1114    #Get output file
1115    if basename_out == None:
1116        ptsname = root + '.pts'
1117    else:
1118        ptsname = basename_out + '.pts'
1119
1120    if verbose: print 'Store to NetCDF file %s' %ptsname
1121    # NetCDF file definition
1122    outfile = NetCDFFile(ptsname, 'w')
1123
1124    #Create new file
1125    outfile.institution = 'Geoscience Australia'
1126    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
1127                      'of spatial point data'
[1577]1128    #assign default values
1129    if easting_min is None: easting_min = xllcorner
1130    if easting_max is None: easting_max = xllcorner + ncols*cellsize
1131    if northing_min is None: northing_min = yllcorner
1132    if northing_max is None: northing_max = yllcorner + nrows*cellsize
[1502]1133
[1577]1134    #compute offsets to update georeferencing
1135    easting_offset = xllcorner - easting_min
1136    northing_offset = yllcorner - northing_min
1137
[1502]1138    #Georeferencing
1139    outfile.zone = zone
[1577]1140    outfile.xllcorner = easting_min #Easting of lower left corner
1141    outfile.yllcorner = northing_min #Northing of lower left corner
[1502]1142    outfile.false_easting = false_easting
1143    outfile.false_northing = false_northing
1144
1145    outfile.projection = projection
1146    outfile.datum = datum
1147    outfile.units = units
1148
1149
1150    #Grid info (FIXME: probably not going to be used, but heck)
1151    outfile.ncols = ncols
1152    outfile.nrows = nrows
1153
1154
1155    # dimension definitions
[1660]1156    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
1157    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
[1577]1158    outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box)
[1502]1159    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1160
1161    # variable definitions
1162    outfile.createVariable('points', Float, ('number_of_points',
1163                                             'number_of_dimensions'))
1164    outfile.createVariable('elevation', Float, ('number_of_points',))
1165
1166    # Get handles to the variables
1167    points = outfile.variables['points']
1168    elevation = outfile.variables['elevation']
1169
[1740]1170    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
1171
[1502]1172    #Store data
[1660]1173    #FIXME: Could perhaps be faster using array operations (Fixed 27/7/05)
1174    global_index = 0
[1502]1175    for i in range(nrows):
[1787]1176        if verbose and i%((nrows+10)/10)==0:
1177            print 'Processing row %d of %d' %(i, nrows)       
1178       
[1660]1179        lower_index = global_index
1180        tpoints = zeros((ncols_in_bounding_box, 2), Float)
1181        telev =  zeros(ncols_in_bounding_box, Float)
1182        local_index = 0
[1502]1183
[1577]1184        y = (nrows-i)*cellsize + yllcorner
[1502]1185        for j in range(ncols):
1186
[1577]1187            x = j*cellsize + xllcorner
1188            if easting_min <= x <= easting_max and \
1189               northing_min <= y <= northing_max:
[1660]1190                tpoints[local_index, :] = [x-easting_min,y-northing_min]
[1740]1191                telev[local_index] = dem_elevation_r[i, j]
[1660]1192                global_index += 1
1193                local_index += 1
[1502]1194
[1660]1195        upper_index = global_index
[1577]1196
[1660]1197        if upper_index == lower_index + ncols_in_bounding_box:
1198            points[lower_index:upper_index, :] = tpoints
1199            elevation[lower_index:upper_index] = telev
1200
1201    assert global_index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points'
1202
[1502]1203    infile.close()
1204    outfile.close()
1205
1206
1207
[1865]1208def sww2dem(basename_in, basename_out = None,
1209            quantity = None,
1210            timestep = None,
1211            reduction = None,
1212            cellsize = 10,
[1875]1213            easting_min = None,
1214            easting_max = None,
1215            northing_min = None,
1216            northing_max = None,         
[1865]1217            verbose = False,
1218            origin = None,
1219            datum = 'WGS84',
1220            format = 'ers'):
1221   
1222    """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
1223
1224    Example (ASC):
1225
1226    ncols         3121
1227    nrows         1800
1228    xllcorner     722000
1229    yllcorner     5893000
1230    cellsize      25
1231    NODATA_value  -9999
1232    138.3698 137.4194 136.5062 135.5558 ..........
1233
1234    Also write accompanying file with same basename_in but extension .prj
1235    used to fix the UTM zone, datum, false northings and eastings.
1236
1237    The prj format is assumed to be as
1238
1239    Projection    UTM
1240    Zone          56
1241    Datum         WGS84
1242    Zunits        NO
1243    Units         METERS
1244    Spheroid      WGS84
1245    Xshift        0.0000000000
1246    Yshift        10000000.0000000000
1247    Parameters
1248
1249
1250    if quantity is given, out values from quantity otherwise default to
1251    elevation
1252
1253    if timestep (an index) is given, output quantity at that timestep
1254
1255    if reduction is given use that to reduce quantity over all timesteps.
1256
1257    datum
1258   
1259    format can be either 'asc' or 'ers'
1260    """
1261
1262    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
1263
1264    msg = 'Format must be either asc or ers'
1265    assert format.lower() in ['asc', 'ers'], msg
1266   
1267    false_easting = 500000
1268    false_northing = 10000000
1269
1270    if quantity is None:
1271        quantity = 'elevation'
1272
1273    if reduction is None:
1274        reduction = max
1275
1276    if basename_out is None:
1277        basename_out = basename_in + '_%s' %quantity
1278 
1279    swwfile = basename_in + '.sww'
1280    demfile = basename_out + '.' + format
1281    # Note the use of a .ers extension is optional (write_ermapper_grid will
1282    # deal with either option
1283
[1875]1284    #Read sww file
[1865]1285    if verbose: print 'Reading from %s' %swwfile
1286    from Scientific.IO.NetCDF import NetCDFFile
1287    fid = NetCDFFile(swwfile)
1288
1289    #Get extent and reference
1290    x = fid.variables['x'][:]
1291    y = fid.variables['y'][:]
1292    volumes = fid.variables['volumes'][:]
1293
1294    number_of_timesteps = fid.dimensions['number_of_timesteps']
1295    number_of_points = fid.dimensions['number_of_points']
1296    if origin is None:
1297
1298        #Get geo_reference
1299        #sww files don't have to have a geo_ref
1300        try:
1301            geo_reference = Geo_reference(NetCDFObject=fid)
1302        except AttributeError, e:
1303            geo_reference = Geo_reference() #Default georef object
1304
1305        xllcorner = geo_reference.get_xllcorner()
1306        yllcorner = geo_reference.get_yllcorner()
1307        zone = geo_reference.get_zone()
1308    else:
1309        zone = origin[0]
1310        xllcorner = origin[1]
1311        yllcorner = origin[2]
1312
1313
[1875]1314
[1865]1315    #Get quantity and reduce if applicable
1316    if verbose: print 'Reading quantity %s' %quantity
1317
1318    if quantity.lower() == 'depth':
1319        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
1320    else:
1321        q = fid.variables[quantity][:]
1322
1323
1324    if len(q.shape) == 2:
[1875]1325        #q has a time component and needs to be reduced along
1326        #the temporal dimension
[1865]1327        if verbose: print 'Reducing quantity %s' %quantity
1328        q_reduced = zeros( number_of_points, Float )
1329
1330        for k in range(number_of_points):
1331            q_reduced[k] = reduction( q[:,k] )
1332
1333        q = q_reduced
1334
[1875]1335    #Post condition: Now q has dimension: number_of_points
1336    assert len(q.shape) == 1
1337    assert q.shape[0] == number_of_points   
1338   
[1865]1339
1340    #Create grid and update xll/yll corner and x,y
[1875]1341
1342    #Relative extent
1343    if easting_min is None:
1344        xmin = min(x)
1345    else:
1346        xmin = easting_min - xllcorner
1347
1348    if easting_max is None:
1349        xmax = max(x)
1350    else:
1351        xmax = easting_max - xllcorner
1352       
1353    if northing_min is None:       
1354        ymin = min(y)
1355    else:
1356        ymin = northing_min - yllcorner
1357       
1358    if northing_max is None:       
1359        ymax = max(y)
1360    else:
1361        ymax = northing_max - yllcorner
1362       
1363       
1364   
[1865]1365    if verbose: print 'Creating grid'
1366    ncols = int((xmax-xmin)/cellsize)+1
1367    nrows = int((ymax-ymin)/cellsize)+1
1368
[1875]1369
1370    #New absolute reference and coordinates
[1865]1371    newxllcorner = xmin+xllcorner
1372    newyllcorner = ymin+yllcorner
1373
1374    x = x+xllcorner-newxllcorner
1375    y = y+yllcorner-newyllcorner
1376
1377    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
1378    assert len(vertex_points.shape) == 2
1379
1380
1381    from Numeric import zeros, Float
1382    grid_points = zeros ( (ncols*nrows, 2), Float )
1383
1384
1385    for i in xrange(nrows):
1386        if format.lower() == 'asc':
1387            yg = i*cellsize
1388        else:   
1389            #this will flip the order of the y values for ers
1390            yg = (nrows-i)*cellsize
1391           
1392        for j in xrange(ncols):
1393            xg = j*cellsize
1394            k = i*ncols + j
1395
1396            grid_points[k,0] = xg
1397            grid_points[k,1] = yg
1398
1399    #Interpolate
1400    from least_squares import Interpolation
1401    from util import inside_polygon
1402
1403    #FIXME: This should be done with precrop = True (?), otherwise it'll
1404    #take forever. With expand_search  set to False, some grid points might
[1875]1405    #miss out.... This will be addressed though Duncan's refactoring of least_squares
[1865]1406
1407    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
1408                           precrop = False, expand_search = False,
1409                           verbose = verbose)
1410
1411    #Interpolate using quantity values
1412    if verbose: print 'Interpolating'
1413    grid_values = interp.interpolate(q).flat
1414
1415    if format.lower() == 'ers':
1416        # setup ERS header information
1417        grid_values = reshape(grid_values,(nrows, ncols))   
1418        NODATA_value = 0
1419        header = {}
1420        header['datum'] = '"' + datum + '"'
1421        # FIXME The use of hardwired UTM and zone number needs to be made optional
1422        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
1423        header['projection'] = '"UTM-' + str(zone) + '"'
1424        header['coordinatetype'] = 'EN'
1425        if header['coordinatetype'] == 'LL':
1426            header['longitude'] = str(newxllcorner)
1427            header['latitude'] = str(newyllcorner)
1428        elif header['coordinatetype'] == 'EN':
1429            header['eastings'] = str(newxllcorner)
1430            header['northings'] = str(newyllcorner)
1431        header['nullcellvalue'] = str(NODATA_value)
1432        header['xdimension'] = str(cellsize)
1433        header['ydimension'] = str(cellsize)
1434        header['value'] = '"' + quantity + '"'
1435
1436
1437        #Write
[1874]1438        if verbose: print 'Writing %s' %demfile
[1865]1439        import ermapper_grids
1440        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
1441
1442        fid.close()   
1443    else:
1444        #Write to Ascii format
1445       
1446        #Write prj file
1447        prjfile = basename_out + '.prj' 
1448       
1449        if verbose: print 'Writing %s' %prjfile
1450        prjid = open(prjfile, 'w')
1451        prjid.write('Projection    %s\n' %'UTM')
1452        prjid.write('Zone          %d\n' %zone)
1453        prjid.write('Datum         %s\n' %datum)
1454        prjid.write('Zunits        NO\n')
1455        prjid.write('Units         METERS\n')
1456        prjid.write('Spheroid      %s\n' %datum)
1457        prjid.write('Xshift        %d\n' %false_easting)
1458        prjid.write('Yshift        %d\n' %false_northing)
1459        prjid.write('Parameters\n')
1460        prjid.close()
1461
1462
1463   
1464        if verbose: print 'Writing %s' %ascfile
1465        NODATA_value = -9999
1466 
1467        ascid = open(demfile, 'w')
1468
1469        ascid.write('ncols         %d\n' %ncols)
1470        ascid.write('nrows         %d\n' %nrows)
1471        ascid.write('xllcorner     %d\n' %newxllcorner)
1472        ascid.write('yllcorner     %d\n' %newyllcorner)
1473        ascid.write('cellsize      %f\n' %cellsize)
1474        ascid.write('NODATA_value  %d\n' %NODATA_value)
1475
1476
1477        #Get bounding polygon from mesh
1478        P = interp.mesh.get_boundary_polygon()
1479        inside_indices = inside_polygon(grid_points, P)
1480
1481        for i in range(nrows):
1482            if verbose and i%((nrows+10)/10)==0:
1483                print 'Doing row %d of %d' %(i, nrows)
1484
1485            for j in range(ncols):
1486                index = (nrows-i-1)*ncols+j
1487
1488                if sometrue(inside_indices == index):
1489                    ascid.write('%f ' %grid_values[index])
1490                else:
1491                    ascid.write('%d ' %NODATA_value)
1492
1493            ascid.write('\n')
1494
1495        #Close
1496        ascid.close()
1497        fid.close()
1498
[1866]1499#Backwards compatibility       
1500def sww2asc(basename_in, basename_out = None,
1501            quantity = None,
1502            timestep = None,
1503            reduction = None,
1504            cellsize = 10,
1505            verbose = False,
1506            origin = None):
1507    print 'sww2asc will soon be obsoleted - please use sww2dem'
[1875]1508    sww2dem(basename_in,
1509            basename_out = basename_out,
1510            quantity = quantity,
1511            timestep = timestep,
1512            reduction = reduction,
1513            cellsize = cellsize,
1514            verbose = verbose,
1515            origin = origin,
[1866]1516            datum = 'WGS84',
1517            format = 'asc')
1518           
1519def sww2ers(basename_in, basename_out = None,
1520            quantity = None,
1521            timestep = None,
1522            reduction = None,
1523            cellsize = 10,
1524            verbose = False,
1525            origin = None,
1526            datum = 'WGS84'):
1527    print 'sww2ers will soon be obsoleted - please use sww2dem'
1528    sww2dem(basename_in, 
[1875]1529            basename_out = basename_out,
1530            quantity = quantity,
1531            timestep = timestep,
1532            reduction = reduction,
1533            cellsize = cellsize,
1534            verbose = verbose,
1535            origin = origin,
1536            datum = datum,
[1866]1537            format = 'ers')         
1538################################# END COMPATIBILITY ##############         
[1865]1539
1540
1541
[1866]1542
[1502]1543def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
1544                                  verbose=False):
1545    """Read Digitial Elevation model from the following ASCII format (.asc)
1546
1547    Example:
1548
1549    ncols         3121
1550    nrows         1800
1551    xllcorner     722000
1552    yllcorner     5893000
1553    cellsize      25
1554    NODATA_value  -9999
1555    138.3698 137.4194 136.5062 135.5558 ..........
1556
1557    Convert basename_in + '.asc' to NetCDF format (.dem)
1558    mimicking the ASCII format closely.
1559
1560
1561    An accompanying file with same basename_in but extension .prj must exist
1562    and is used to fix the UTM zone, datum, false northings and eastings.
1563
1564    The prj format is assumed to be as
1565
1566    Projection    UTM
1567    Zone          56
1568    Datum         WGS84
1569    Zunits        NO
1570    Units         METERS
1571    Spheroid      WGS84
1572    Xshift        0.0000000000
1573    Yshift        10000000.0000000000
1574    Parameters
1575    """
1576
1577    import os
1578    from Scientific.IO.NetCDF import NetCDFFile
1579    from Numeric import Float, array
1580
1581    #root, ext = os.path.splitext(basename_in)
1582    root = basename_in
1583
1584    ###########################################
1585    # Read Meta data
1586    if verbose: print 'Reading METADATA from %s' %root + '.prj'
1587    metadatafile = open(root + '.prj')
1588    metalines = metadatafile.readlines()
1589    metadatafile.close()
1590
1591    L = metalines[0].strip().split()
1592    assert L[0].strip().lower() == 'projection'
1593    projection = L[1].strip()                   #TEXT
1594
1595    L = metalines[1].strip().split()
1596    assert L[0].strip().lower() == 'zone'
1597    zone = int(L[1].strip())
1598
1599    L = metalines[2].strip().split()
1600    assert L[0].strip().lower() == 'datum'
1601    datum = L[1].strip()                        #TEXT
1602
1603    L = metalines[3].strip().split()
1604    assert L[0].strip().lower() == 'zunits'     #IGNORE
1605    zunits = L[1].strip()                       #TEXT
1606
1607    L = metalines[4].strip().split()
1608    assert L[0].strip().lower() == 'units'
1609    units = L[1].strip()                        #TEXT
1610
1611    L = metalines[5].strip().split()
1612    assert L[0].strip().lower() == 'spheroid'   #IGNORE
1613    spheroid = L[1].strip()                     #TEXT
1614
1615    L = metalines[6].strip().split()
1616    assert L[0].strip().lower() == 'xshift'
1617    false_easting = float(L[1].strip())
1618
1619    L = metalines[7].strip().split()
1620    assert L[0].strip().lower() == 'yshift'
1621    false_northing = float(L[1].strip())
1622
1623    #print false_easting, false_northing, zone, datum
1624
1625
1626    ###########################################
1627    #Read DEM data
1628
1629    datafile = open(basename_in + '.asc')
1630
1631    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
1632    lines = datafile.readlines()
1633    datafile.close()
1634
1635    if verbose: print 'Got', len(lines), ' lines'
1636
1637    ncols = int(lines[0].split()[1].strip())
1638    nrows = int(lines[1].split()[1].strip())
1639    xllcorner = float(lines[2].split()[1].strip())
1640    yllcorner = float(lines[3].split()[1].strip())
1641    cellsize = float(lines[4].split()[1].strip())
1642    NODATA_value = int(lines[5].split()[1].strip())
1643
1644    assert len(lines) == nrows + 6
1645
1646
1647    ##########################################
1648
1649
1650    if basename_out == None:
1651        netcdfname = root + '.dem'
1652    else:
1653        netcdfname = basename_out + '.dem'
1654
1655    if verbose: print 'Store to NetCDF file %s' %netcdfname
1656    # NetCDF file definition
1657    fid = NetCDFFile(netcdfname, 'w')
1658
1659    #Create new file
1660    fid.institution = 'Geoscience Australia'
1661    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
1662                      'of spatial point data'
1663
1664    fid.ncols = ncols
1665    fid.nrows = nrows
1666    fid.xllcorner = xllcorner
1667    fid.yllcorner = yllcorner
1668    fid.cellsize = cellsize
1669    fid.NODATA_value = NODATA_value
1670
1671    fid.zone = zone
1672    fid.false_easting = false_easting
1673    fid.false_northing = false_northing
1674    fid.projection = projection
1675    fid.datum = datum
1676    fid.units = units
1677
1678
1679    # dimension definitions
1680    fid.createDimension('number_of_rows', nrows)
1681    fid.createDimension('number_of_columns', ncols)
1682
1683    # variable definitions
1684    fid.createVariable('elevation', Float, ('number_of_rows',
1685                                            'number_of_columns'))
1686
1687    # Get handles to the variables
1688    elevation = fid.variables['elevation']
1689
1690    #Store data
[1787]1691    n = len(lines[6:])
[1502]1692    for i, line in enumerate(lines[6:]):
1693        fields = line.split()
[1787]1694        if verbose and i%((n+10)/10)==0:
1695            print 'Processing row %d of %d' %(i, nrows)           
[1502]1696
1697        elevation[i, :] = array([float(x) for x in fields])
1698
1699    fid.close()
1700
1701
1702
[1797]1703
1704
[1502]1705def ferret2sww(basename_in, basename_out = None,
1706               verbose = False,
1707               minlat = None, maxlat = None,
1708               minlon = None, maxlon = None,
1709               mint = None, maxt = None, mean_stage = 0,
1710               origin = None, zscale = 1,
1711               fail_on_NaN = True,
1712               NaN_filler = 0,
1713               elevation = None,
1714               inverted_bathymetry = False
1715               ): #FIXME: Bathymetry should be obtained
1716                                  #from MOST somehow.
1717                                  #Alternatively from elsewhere
1718                                  #or, as a last resort,
1719                                  #specified here.
1720                                  #The value of -100 will work
1721                                  #for the Wollongong tsunami
1722                                  #scenario but is very hacky
1723    """Convert 'Ferret' NetCDF format for wave propagation to
1724    sww format native to pyvolution.
1725
1726    Specify only basename_in and read files of the form
1727    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
1728    relative height, x-velocity and y-velocity, respectively.
1729
1730    Also convert latitude and longitude to UTM. All coordinates are
1731    assumed to be given in the GDA94 datum.
1732
1733    min's and max's: If omitted - full extend is used.
1734    To include a value min may equal it, while max must exceed it.
1735    Lat and lon are assuemd to be in decimal degrees
1736
1737    origin is a 3-tuple with geo referenced
1738    UTM coordinates (zone, easting, northing)
1739
1740    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
1741    which means that longitude is the fastest
1742    varying dimension (row major order, so to speak)
1743
1744    ferret2sww uses grid points as vertices in a triangular grid
1745    counting vertices from lower left corner upwards, then right
1746    """
1747
1748    import os
1749    from Scientific.IO.NetCDF import NetCDFFile
1750    from Numeric import Float, Int, Int32, searchsorted, zeros, array
[1854]1751    from Numeric import allclose, around
1752       
[1502]1753    precision = Float
1754
1755
1756    #Get NetCDF data
1757    if verbose: print 'Reading files %s_*.nc' %basename_in
1758    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
1759    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
1760    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
1761    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)
1762
1763    if basename_out is None:
1764        swwname = basename_in + '.sww'
1765    else:
1766        swwname = basename_out + '.sww'
1767
1768    times = file_h.variables['TIME']
1769    latitudes = file_h.variables['LAT']
1770    longitudes = file_h.variables['LON']
1771
[1797]1772
1773
[1854]1774    #Precision used by most for lat/lon is 4 or 5 decimals
1775    e_lat = around(file_e.variables['LAT'][:], 5)
1776    e_lon = around(file_e.variables['LON'][:], 5)
[1797]1777
1778    #Check that files are compatible
1779    assert allclose(latitudes, file_u.variables['LAT'])
[1854]1780    assert allclose(latitudes, file_v.variables['LAT'])
1781    assert allclose(latitudes, e_lat)
[1797]1782
1783    assert allclose(longitudes, file_u.variables['LON'])
1784    assert allclose(longitudes, file_v.variables['LON'])
[1854]1785    assert allclose(longitudes, e_lon)   
[1797]1786
1787
1788
[1502]1789    if mint == None:
1790        jmin = 0
1791    else:
1792        jmin = searchsorted(times, mint)
1793
1794    if maxt == None:
1795        jmax=len(times)
1796    else:
1797        jmax = searchsorted(times, maxt)
1798
1799    if minlat == None:
1800        kmin=0
1801    else:
1802        kmin = searchsorted(latitudes, minlat)
1803
1804    if maxlat == None:
1805        kmax = len(latitudes)
1806    else:
1807        kmax = searchsorted(latitudes, maxlat)
1808
1809    if minlon == None:
1810        lmin=0
1811    else:
1812        lmin = searchsorted(longitudes, minlon)
1813
1814    if maxlon == None:
1815        lmax = len(longitudes)
1816    else:
1817        lmax = searchsorted(longitudes, maxlon)
1818
1819
1820
1821    times = times[jmin:jmax]
1822    latitudes = latitudes[kmin:kmax]
1823    longitudes = longitudes[lmin:lmax]
1824
1825
1826    if verbose: print 'cropping'
[1797]1827    zname = 'ELEVATION'
1828
1829   
[1502]1830    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
1831    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
1832    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
[1797]1833    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
[1502]1834
1835#    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
1836#        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
1837#    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
1838#        from Numeric import asarray
1839#        elevations=elevations.tolist()
1840#        elevations.reverse()
1841#        elevations=asarray(elevations)
1842#    else:
1843#        from Numeric import asarray
1844#        elevations=elevations.tolist()
1845#        elevations.reverse()
1846#        elevations=asarray(elevations)
1847#        'print hmmm'
1848
1849
1850
1851    #Get missing values
1852    nan_ha = file_h.variables['HA'].missing_value[0]
1853    nan_ua = file_u.variables['UA'].missing_value[0]
1854    nan_va = file_v.variables['VA'].missing_value[0]
[1797]1855    if hasattr(file_e.variables[zname],'missing_value'):
1856        nan_e  = file_e.variables[zname].missing_value[0]
[1502]1857    else:
1858        nan_e = None
1859
1860    #Cleanup
1861    from Numeric import sometrue
1862
1863    missing = (amplitudes == nan_ha)
1864    if sometrue (missing):
1865        if fail_on_NaN:
1866            msg = 'NetCDFFile %s contains missing values'\
1867                  %(basename_in+'_ha.nc')
1868            raise msg
1869        else:
1870            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
1871
1872    missing = (uspeed == nan_ua)
1873    if sometrue (missing):
1874        if fail_on_NaN:
1875            msg = 'NetCDFFile %s contains missing values'\
1876                  %(basename_in+'_ua.nc')
1877            raise msg
1878        else:
1879            uspeed = uspeed*(missing==0) + missing*NaN_filler
1880
1881    missing = (vspeed == nan_va)
1882    if sometrue (missing):
1883        if fail_on_NaN:
1884            msg = 'NetCDFFile %s contains missing values'\
1885                  %(basename_in+'_va.nc')
1886            raise msg
1887        else:
1888            vspeed = vspeed*(missing==0) + missing*NaN_filler
1889
1890
1891    missing = (elevations == nan_e)
1892    if sometrue (missing):
1893        if fail_on_NaN:
1894            msg = 'NetCDFFile %s contains missing values'\
1895                  %(basename_in+'_e.nc')
1896            raise msg
1897        else:
1898            elevations = elevations*(missing==0) + missing*NaN_filler
1899
1900    #######
1901
1902
1903
1904    number_of_times = times.shape[0]
1905    number_of_latitudes = latitudes.shape[0]
1906    number_of_longitudes = longitudes.shape[0]
1907
1908    assert amplitudes.shape[0] == number_of_times
1909    assert amplitudes.shape[1] == number_of_latitudes
1910    assert amplitudes.shape[2] == number_of_longitudes
1911
1912    if verbose:
1913        print '------------------------------------------------'
1914        print 'Statistics:'
1915        print '  Extent (lat/lon):'
1916        print '    lat in [%f, %f], len(lat) == %d'\
1917              %(min(latitudes.flat), max(latitudes.flat),
1918                len(latitudes.flat))
1919        print '    lon in [%f, %f], len(lon) == %d'\
1920              %(min(longitudes.flat), max(longitudes.flat),
1921                len(longitudes.flat))
1922        print '    t in [%f, %f], len(t) == %d'\
1923              %(min(times.flat), max(times.flat), len(times.flat))
1924
1925        q = amplitudes.flat
1926        name = 'Amplitudes (ha) [cm]'
1927        print %s in [%f, %f]' %(name, min(q), max(q))
1928
1929        q = uspeed.flat
1930        name = 'Speeds (ua) [cm/s]'
1931        print %s in [%f, %f]' %(name, min(q), max(q))
1932
1933        q = vspeed.flat
1934        name = 'Speeds (va) [cm/s]'
1935        print %s in [%f, %f]' %(name, min(q), max(q))
1936
1937        q = elevations.flat
1938        name = 'Elevations (e) [m]'
1939        print %s in [%f, %f]' %(name, min(q), max(q))
1940
1941
1942    #print number_of_latitudes, number_of_longitudes
1943    number_of_points = number_of_latitudes*number_of_longitudes
1944    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
1945
1946
1947    file_h.close()
1948    file_u.close()
1949    file_v.close()
1950    file_e.close()
1951
1952
1953    # NetCDF file definition
1954    outfile = NetCDFFile(swwname, 'w')
1955
1956    #Create new file
1957    outfile.institution = 'Geoscience Australia'
1958    outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
1959                          %(basename_in + '_ha.nc',
1960                            basename_in + '_ua.nc',
1961                            basename_in + '_va.nc',
1962                            basename_in + '_e.nc')
1963
1964
1965    #For sww compatibility
1966    outfile.smoothing = 'Yes'
1967    outfile.order = 1
1968
1969    #Start time in seconds since the epoch (midnight 1/1/1970)
1970    outfile.starttime = starttime = times[0]
1971    times = times - starttime  #Store relative times
1972
1973    # dimension definitions
1974    outfile.createDimension('number_of_volumes', number_of_volumes)
1975
1976    outfile.createDimension('number_of_vertices', 3)
1977    outfile.createDimension('number_of_points', number_of_points)
1978
1979
1980    #outfile.createDimension('number_of_timesteps', len(times))
1981    outfile.createDimension('number_of_timesteps', len(times))
1982
1983    # variable definitions
1984    outfile.createVariable('x', precision, ('number_of_points',))
1985    outfile.createVariable('y', precision, ('number_of_points',))
1986    outfile.createVariable('elevation', precision, ('number_of_points',))
1987
1988    #FIXME: Backwards compatibility
1989    outfile.createVariable('z', precision, ('number_of_points',))
1990    #################################
1991
1992    outfile.createVariable('volumes', Int, ('number_of_volumes',
1993                                            'number_of_vertices'))
1994
1995    outfile.createVariable('time', precision,
1996                           ('number_of_timesteps',))
1997
1998    outfile.createVariable('stage', precision,
1999                           ('number_of_timesteps',
2000                            'number_of_points'))
2001
2002    outfile.createVariable('xmomentum', precision,
2003                           ('number_of_timesteps',
2004                            'number_of_points'))
2005
2006    outfile.createVariable('ymomentum', precision,
2007                           ('number_of_timesteps',
2008                            'number_of_points'))
2009
2010
2011    #Store
2012    from coordinate_transforms.redfearn import redfearn
2013    x = zeros(number_of_points, Float)  #Easting
2014    y = zeros(number_of_points, Float)  #Northing
2015
2016
2017    if verbose: print 'Making triangular grid'
2018    #Check zone boundaries
2019    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
2020
2021    vertices = {}
2022    i = 0
2023    for k, lat in enumerate(latitudes):       #Y direction
2024        for l, lon in enumerate(longitudes):  #X direction
2025
2026            vertices[l,k] = i
2027
2028            zone, easting, northing = redfearn(lat,lon)
2029
2030            msg = 'Zone boundary crossed at longitude =', lon
2031            #assert zone == refzone, msg
2032            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
2033            x[i] = easting
2034            y[i] = northing
2035            i += 1
2036
2037
2038    #Construct 2 triangles per 'rectangular' element
2039    volumes = []
2040    for l in range(number_of_longitudes-1):    #X direction
2041        for k in range(number_of_latitudes-1): #Y direction
2042            v1 = vertices[l,k+1]
2043            v2 = vertices[l,k]
2044            v3 = vertices[l+1,k+1]
2045            v4 = vertices[l+1,k]
2046
2047            volumes.append([v1,v2,v3]) #Upper element
2048            volumes.append([v4,v3,v2]) #Lower element
2049
2050    volumes = array(volumes)
2051
2052    if origin == None:
2053        zone = refzone
2054        xllcorner = min(x)
2055        yllcorner = min(y)
2056    else:
2057        zone = origin[0]
2058        xllcorner = origin[1]
2059        yllcorner = origin[2]
2060
2061
2062    outfile.xllcorner = xllcorner
2063    outfile.yllcorner = yllcorner
2064    outfile.zone = zone
2065
2066
2067    if elevation is not None:
2068        z = elevation
2069    else:
2070        if inverted_bathymetry:
2071            z = -1*elevations
2072        else:
2073            z = elevations
2074        #FIXME: z should be obtained from MOST and passed in here
2075
2076    from Numeric import resize
2077    z = resize(z,outfile.variables['z'][:].shape)
2078    outfile.variables['x'][:] = x - xllcorner
2079    outfile.variables['y'][:] = y - yllcorner
2080    outfile.variables['z'][:] = z
2081    outfile.variables['elevation'][:] = z  #FIXME HACK
2082    outfile.variables['time'][:] = times   #Store time relative
2083    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
2084
2085
2086
2087    #Time stepping
2088    stage = outfile.variables['stage']
2089    xmomentum = outfile.variables['xmomentum']
2090    ymomentum = outfile.variables['ymomentum']
2091
2092    if verbose: print 'Converting quantities'
2093    n = len(times)
2094    for j in range(n):
2095        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
2096        i = 0
2097        for k in range(number_of_latitudes):      #Y direction
2098            for l in range(number_of_longitudes): #X direction
2099                w = zscale*amplitudes[j,k,l]/100 + mean_stage
2100                stage[j,i] = w
2101                h = w - z[i]
2102                xmomentum[j,i] = uspeed[j,k,l]/100*h
2103                ymomentum[j,i] = vspeed[j,k,l]/100*h
2104                i += 1
2105
2106
2107    if verbose:
2108        x = outfile.variables['x'][:]
2109        y = outfile.variables['y'][:]
2110        print '------------------------------------------------'
2111        print 'Statistics of output file:'
2112        print '  Name: %s' %swwname
2113        print '  Reference:'
2114        print '    Lower left corner: [%f, %f]'\
2115              %(xllcorner, yllcorner)
2116        print '    Start time: %f' %starttime
2117        print '  Extent:'
2118        print '    x [m] in [%f, %f], len(x) == %d'\
2119              %(min(x.flat), max(x.flat), len(x.flat))
2120        print '    y [m] in [%f, %f], len(y) == %d'\
2121              %(min(y.flat), max(y.flat), len(y.flat))
2122        print '    t [s] in [%f, %f], len(t) == %d'\
2123              %(min(times), max(times), len(times))
2124        print '  Quantities [SI units]:'
[1797]2125        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2126            q = outfile.variables[name][:].flat
[1502]2127            print '    %s in [%f, %f]' %(name, min(q), max(q))
2128
2129
2130
2131
2132    outfile.close()
2133
2134
2135
[1671]2136
[1835]2137def timefile2netcdf(filename, quantity_names = None):
[1671]2138    """Template for converting typical text files with time series to
[1835]2139    NetCDF tms file.
[1671]2140
2141
2142    The file format is assumed to be either two fields separated by a comma:
2143
2144        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
2145
2146    E.g
2147
2148      31/08/04 00:00:00, 1.328223 0 0
2149      31/08/04 00:15:00, 1.292912 0 0
2150
2151    will provide a time dependent function f(t) with three attributes
2152
2153    filename is assumed to be the rootname with extenisons .txt and .sww
2154    """
2155
2156    import time, calendar
2157    from Numeric import array
2158    from config import time_format
2159    from util import ensure_numeric
2160
2161
2162    fid = open(filename + '.txt')
2163    line = fid.readline()
2164    fid.close()
2165
2166    fields = line.split(',')
2167    msg = 'File %s must have the format date, value0 value1 value2 ...'
2168    assert len(fields) == 2, msg
2169
2170    try:
2171        starttime = calendar.timegm(time.strptime(fields[0], time_format))
2172    except ValueError:
2173        msg = 'First field in file %s must be' %filename
2174        msg += ' date-time with format %s.\n' %time_format
2175        msg += 'I got %s instead.' %fields[0]
2176        raise msg
2177
2178
2179    #Split values
2180    values = []
2181    for value in fields[1].split():
2182        values.append(float(value))
2183
2184    q = ensure_numeric(values)
2185
2186    msg = 'ERROR: File must contain at least one independent value'
2187    assert len(q.shape) == 1, msg
2188
2189
2190
2191    #Read times proper
2192    from Numeric import zeros, Float, alltrue
2193    from config import time_format
2194    import time, calendar
2195
2196    fid = open(filename + '.txt')
2197    lines = fid.readlines()
2198    fid.close()
2199
2200    N = len(lines)
2201    d = len(q)
2202
2203    T = zeros(N, Float)       #Time
2204    Q = zeros((N, d), Float)  #Values
2205
2206    for i, line in enumerate(lines):
2207        fields = line.split(',')
2208        realtime = calendar.timegm(time.strptime(fields[0], time_format))
2209       
2210        T[i] = realtime - starttime
2211
2212        for j, value in enumerate(fields[1].split()):
2213            Q[i, j] = float(value)
2214
2215    msg = 'File %s must list time as a monotonuosly ' %filename
2216    msg += 'increasing sequence'
2217    assert alltrue( T[1:] - T[:-1] > 0 ), msg
2218
2219
[1835]2220    #Create NetCDF file
[1671]2221    from Scientific.IO.NetCDF import NetCDFFile
2222
[1835]2223    fid = NetCDFFile(filename + '.tms', 'w')
[1671]2224
2225
2226    fid.institution = 'Geoscience Australia'
2227    fid.description = 'Time series'
2228
2229
2230    #Reference point
2231    #Start time in seconds since the epoch (midnight 1/1/1970)
2232    #FIXME: Use Georef
2233    fid.starttime = starttime
2234   
2235   
2236    # dimension definitions
2237    #fid.createDimension('number_of_volumes', self.number_of_volumes)
2238    #fid.createDimension('number_of_vertices', 3)
2239
2240
2241    fid.createDimension('number_of_timesteps', len(T))
2242
2243    fid.createVariable('time', Float, ('number_of_timesteps',))
2244
2245    fid.variables['time'][:] = T
2246   
2247    for i in range(Q.shape[1]):
2248        try:
2249            name = quantity_names[i]
2250        except:   
2251            name = 'Attribute%d'%i
2252           
2253        fid.createVariable(name, Float, ('number_of_timesteps',))
2254        fid.variables[name][:] = Q[:,i]       
2255
2256    fid.close()
2257
2258
[1502]2259def extent_sww(file_name):
2260    """
2261    Read in an sww file.
2262
2263    Input;
2264    file_name - the sww file
2265
2266    Output;
2267    z - Vector of bed elevation
2268    volumes - Array.  Each row has 3 values, representing
2269    the vertices that define the volume
2270    time - Vector of the times where there is stage information
2271    stage - array with respect to time and vertices (x,y)
2272    """
2273
2274
2275    from Scientific.IO.NetCDF import NetCDFFile
2276
2277    #Check contents
2278    #Get NetCDF
2279    fid = NetCDFFile(file_name, 'r')
2280
2281    # Get the variables
2282    x = fid.variables['x'][:]
2283    y = fid.variables['y'][:]
2284    stage = fid.variables['stage'][:]
2285    #print "stage",stage
2286    #print "stage.shap",stage.shape
2287    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
2288    #print "min(stage)",min(stage)
2289
2290    fid.close()
2291
2292    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
2293
2294
2295def sww2domain(filename,boundary=None,t=None,\
2296               fail_if_NaN=True,NaN_filler=0\
2297               ,verbose = True,very_verbose = False):
2298    """
2299    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
2300
2301    Boundary is not recommended if domian.smooth is not selected, as it
2302    uses unique coordinates, but not unique boundaries. This means that
2303    the boundary file will not be compatable with the coordiantes, and will
2304    give a different final boundary, or crash.
2305    """
2306    NaN=9.969209968386869e+036
2307    #initialise NaN.
2308
2309    from Scientific.IO.NetCDF import NetCDFFile
2310    from shallow_water import Domain
2311    from Numeric import asarray, transpose, resize
2312
2313    if verbose: print 'Reading from ', filename
2314    fid = NetCDFFile(filename, 'r')    #Open existing file for read
2315    time = fid.variables['time']       #Timesteps
2316    if t is None:
2317        t = time[-1]
2318    time_interp = get_time_interp(time,t)
2319
2320    # Get the variables as Numeric arrays
2321    x = fid.variables['x'][:]             #x-coordinates of vertices
2322    y = fid.variables['y'][:]             #y-coordinates of vertices
2323    elevation = fid.variables['elevation']     #Elevation
2324    stage = fid.variables['stage']     #Water level
2325    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
2326    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
2327
2328    starttime = fid.starttime[0]
2329    volumes = fid.variables['volumes'][:] #Connectivity
2330    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
2331
2332    conserved_quantities = []
2333    interpolated_quantities = {}
2334    other_quantities = []
2335
2336    # get geo_reference
2337    #sww files don't have to have a geo_ref
2338    try:
2339        geo_reference = Geo_reference(NetCDFObject=fid)
2340    except: #AttributeError, e:
2341        geo_reference = None
2342
2343    if verbose: print '    getting quantities'
2344    for quantity in fid.variables.keys():
2345        dimensions = fid.variables[quantity].dimensions
2346        if 'number_of_timesteps' in dimensions:
2347            conserved_quantities.append(quantity)
2348            interpolated_quantities[quantity]=\
2349                  interpolated_quantity(fid.variables[quantity][:],time_interp)
2350        else: other_quantities.append(quantity)
2351
2352    other_quantities.remove('x')
2353    other_quantities.remove('y')
2354    other_quantities.remove('z')
2355    other_quantities.remove('volumes')
2356
2357    conserved_quantities.remove('time')
2358
2359    if verbose: print '    building domain'
2360#    From domain.Domain:
2361#    domain = Domain(coordinates, volumes,\
2362#                    conserved_quantities = conserved_quantities,\
2363#                    other_quantities = other_quantities,zone=zone,\
2364#                    xllcorner=xllcorner, yllcorner=yllcorner)
2365
2366#   From shallow_water.Domain:
2367    coordinates=coordinates.tolist()
2368    volumes=volumes.tolist()
2369    #FIXME:should this be in mesh?(peter row)
2370    if fid.smoothing == 'Yes': unique = False
2371    else: unique = True
2372    if unique:
2373        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
2374
2375
2376    domain = Domain(coordinates, volumes, boundary)
2377
2378    if not boundary is None:
2379        domain.boundary = boundary
2380
2381    domain.geo_reference = geo_reference
2382
2383    domain.starttime=float(starttime)+float(t)
2384    domain.time=0.0
2385
2386    for quantity in other_quantities:
2387        try:
2388            NaN = fid.variables[quantity].missing_value
2389        except:
2390            pass #quantity has no missing_value number
2391        X = fid.variables[quantity][:]
2392        if very_verbose:
2393            print '       ',quantity
2394            print '        NaN =',NaN
2395            print '        max(X)'
2396            print '       ',max(X)
2397            print '        max(X)==NaN'
2398            print '       ',max(X)==NaN
2399            print ''
2400        if (max(X)==NaN) or (min(X)==NaN):
2401            if fail_if_NaN:
2402                msg = 'quantity "%s" contains no_data entry'%quantity
2403                raise msg
2404            else:
2405                data = (X<>NaN)
2406                X = (X*data)+(data==0)*NaN_filler
2407        if unique:
2408            X = resize(X,(len(X)/3,3))
2409        domain.set_quantity(quantity,X)
2410#
2411    for quantity in conserved_quantities:
2412        try:
2413            NaN = fid.variables[quantity].missing_value
2414        except:
2415            pass #quantity has no missing_value number
2416        X = interpolated_quantities[quantity]
2417        if very_verbose:
2418            print '       ',quantity
2419            print '        NaN =',NaN
2420            print '        max(X)'
2421            print '       ',max(X)
2422            print '        max(X)==NaN'
2423            print '       ',max(X)==NaN
2424            print ''
2425        if (max(X)==NaN) or (min(X)==NaN):
2426            if fail_if_NaN:
2427                msg = 'quantity "%s" contains no_data entry'%quantity
2428                raise msg
2429            else:
2430                data = (X<>NaN)
2431                X = (X*data)+(data==0)*NaN_filler
2432        if unique:
2433            X = resize(X,(X.shape[0]/3,3))
2434        domain.set_quantity(quantity,X)
2435    fid.close()
2436    return domain
2437
2438def interpolated_quantity(saved_quantity,time_interp):
2439
2440    #given an index and ratio, interpolate quantity with respect to time.
2441    index,ratio = time_interp
2442    Q = saved_quantity
2443    if ratio > 0:
2444        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
2445    else:
2446        q = Q[index]
2447    #Return vector of interpolated values
2448    return q
2449
2450def get_time_interp(time,t=None):
2451    #Finds the ratio and index for time interpolation.
2452    #It is borrowed from previous pyvolution code.
2453    if t is None:
2454        t=time[-1]
2455        index = -1
2456        ratio = 0.
2457    else:
2458        T = time
2459        tau = t
2460        index=0
2461        msg = 'Time interval derived from file %s [%s:%s]'\
2462            %('FIXMEfilename', T[0], T[-1])
2463        msg += ' does not match model time: %s' %tau
2464        if tau < time[0]: raise msg
2465        if tau > time[-1]: raise msg
2466        while tau > time[index]: index += 1
2467        while tau < time[index]: index -= 1
2468        if tau == time[index]:
2469            #Protect against case where tau == time[-1] (last time)
2470            # - also works in general when tau == time[i]
2471            ratio = 0
2472        else:
2473            #t is now between index and index+1
2474            ratio = (tau - time[index])/(time[index+1] - time[index])
2475    return (index,ratio)
2476
2477
2478def weed(coordinates,volumes,boundary = None):
2479    if type(coordinates)=='array':
2480        coordinates = coordinates.tolist()
2481    if type(volumes)=='array':
2482        volumes = volumes.tolist()
2483
2484    unique = False
2485    point_dict = {}
2486    same_point = {}
2487    for i in range(len(coordinates)):
2488        point = tuple(coordinates[i])
2489        if point_dict.has_key(point):
2490            unique = True
2491            same_point[i]=point
2492            #to change all point i references to point j
2493        else:
2494            point_dict[point]=i
2495            same_point[i]=point
2496
2497    coordinates = []
2498    i = 0
2499    for point in point_dict.keys():
2500        point = tuple(point)
2501        coordinates.append(list(point))
2502        point_dict[point]=i
2503        i+=1
2504
2505
2506    for volume in volumes:
2507        for i in range(len(volume)):
2508            index = volume[i]
2509            if index>-1:
2510                volume[i]=point_dict[same_point[index]]
2511
2512    new_boundary = {}
2513    if not boundary is None:
2514        for segment in boundary.keys():
2515            point0 = point_dict[same_point[segment[0]]]
2516            point1 = point_dict[same_point[segment[1]]]
2517            label = boundary[segment]
2518            #FIXME should the bounday attributes be concaterated
2519            #('exterior, pond') or replaced ('pond')(peter row)
2520
2521            if new_boundary.has_key((point0,point1)):
2522                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
2523                                              #+','+label
2524
2525            elif new_boundary.has_key((point1,point0)):
2526                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
2527                                              #+','+label
2528            else: new_boundary[(point0,point1)]=label
2529
2530        boundary = new_boundary
2531
2532    return coordinates,volumes,boundary
2533
2534
2535
[1835]2536
2537
2538
2539
2540
2541
2542
2543###############################################
[1502]2544#OBSOLETE STUFF
2545#Native checkpoint format.
2546#Information needed to recreate a state is preserved
2547#FIXME: Rethink and maybe use netcdf format
2548def cpt_variable_writer(filename, t, v0, v1, v2):
2549    """Store all conserved quantities to file
2550    """
2551
2552    M, N = v0.shape
2553
2554    FN = create_filename(filename, 'cpt', M, t)
2555    #print 'Writing to %s' %FN
2556
2557    fid = open(FN, 'w')
2558    for i in range(M):
2559        for j in range(N):
2560            fid.write('%.16e ' %v0[i,j])
2561        for j in range(N):
2562            fid.write('%.16e ' %v1[i,j])
2563        for j in range(N):
2564            fid.write('%.16e ' %v2[i,j])
2565
2566        fid.write('\n')
2567    fid.close()
2568
2569
2570def cpt_variable_reader(filename, t, v0, v1, v2):
2571    """Store all conserved quantities to file
2572    """
2573
2574    M, N = v0.shape
2575
2576    FN = create_filename(filename, 'cpt', M, t)
2577    #print 'Reading from %s' %FN
2578
2579    fid = open(FN)
2580
2581
2582    for i in range(M):
2583        values = fid.readline().split() #Get one line
2584
2585        for j in range(N):
2586            v0[i,j] = float(values[j])
2587            v1[i,j] = float(values[3+j])
2588            v2[i,j] = float(values[6+j])
2589
2590    fid.close()
2591
2592def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
2593    """Writes x,y,z,z,z coordinates of triangles constituting the bed
2594    elevation.
[1835]2595    FIXME: Not in use pt
[1502]2596    """
2597
2598    M, N = v0.shape
2599
[1835]2600
[1502]2601    print X0
2602    import sys; sys.exit()
2603    FN = create_filename(filename, 'cpt', M)
2604    print 'Writing to %s' %FN
2605
2606    fid = open(FN, 'w')
2607    for i in range(M):
2608        for j in range(2):
2609            fid.write('%.16e ' %X0[i,j])   #x, y
2610        for j in range(N):
2611            fid.write('%.16e ' %v0[i,j])       #z,z,z,
2612
2613        for j in range(2):
2614            fid.write('%.16e ' %X1[i,j])   #x, y
2615        for j in range(N):
2616            fid.write('%.16e ' %v1[i,j])
2617
2618        for j in range(2):
2619            fid.write('%.16e ' %X2[i,j])   #x, y
2620        for j in range(N):
2621            fid.write('%.16e ' %v2[i,j])
2622
2623        fid.write('\n')
2624    fid.close()
2625
2626
2627
2628#Function for storing out to e.g. visualisation
2629#FIXME: Do we want this?
2630#FIXME: Not done yet for this version
2631def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
2632    """Writes x,y,z coordinates of triangles constituting the bed elevation.
2633    """
2634
2635    M, N = v0.shape
2636
2637    FN = create_filename(filename, 'dat', M)
2638    #print 'Writing to %s' %FN
2639
2640    fid = open(FN, 'w')
2641    for i in range(M):
2642        for j in range(2):
2643            fid.write('%f ' %X0[i,j])   #x, y
2644        fid.write('%f ' %v0[i,0])       #z
2645
2646        for j in range(2):
2647            fid.write('%f ' %X1[i,j])   #x, y
2648        fid.write('%f ' %v1[i,0])       #z
2649
2650        for j in range(2):
2651            fid.write('%f ' %X2[i,j])   #x, y
2652        fid.write('%f ' %v2[i,0])       #z
2653
2654        fid.write('\n')
2655    fid.close()
2656
2657
2658
2659def dat_variable_writer(filename, t, v0, v1, v2):
2660    """Store water height to file
2661    """
2662
2663    M, N = v0.shape
2664
2665    FN = create_filename(filename, 'dat', M, t)
2666    #print 'Writing to %s' %FN
2667
2668    fid = open(FN, 'w')
2669    for i in range(M):
2670        fid.write('%.4f ' %v0[i,0])
2671        fid.write('%.4f ' %v1[i,0])
2672        fid.write('%.4f ' %v2[i,0])
2673
2674        fid.write('\n')
2675    fid.close()
2676
2677
2678def read_sww(filename):
2679    """Read sww Net CDF file containing Shallow Water Wave simulation
2680
2681    The integer array volumes is of shape Nx3 where N is the number of
2682    triangles in the mesh.
2683
2684    Each entry in volumes is an index into the x,y arrays (the location).
2685
2686    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
2687    number_of_timesteps, number_of_points.
2688
2689    The momentum is not always stored.
2690
2691    """
2692    from Scientific.IO.NetCDF import NetCDFFile
2693    print 'Reading from ', filename
2694    fid = NetCDFFile(filename, 'r')    #Open existing file for read
2695#latitude, longitude
2696    # Get the variables as Numeric arrays
2697    x = fid.variables['x']             #x-coordinates of vertices
2698    y = fid.variables['y']             #y-coordinates of vertices
2699    z = fid.variables['elevation']     #Elevation
2700    time = fid.variables['time']       #Timesteps
2701    stage = fid.variables['stage']     #Water level
2702    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
2703    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
2704
2705    volumes = fid.variables['volumes'] #Connectivity
[1660]2706
2707
[1666]2708def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
2709                 verbose=False):
2710    """Read Digitial Elevation model from the following NetCDF format (.dem)
2711
2712    Example:
2713
2714    ncols         3121
2715    nrows         1800
2716    xllcorner     722000
2717    yllcorner     5893000
2718    cellsize      25
2719    NODATA_value  -9999
2720    138.3698 137.4194 136.5062 135.5558 ..........
2721
2722    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
2723    """
2724
2725    import os
2726    from Scientific.IO.NetCDF import NetCDFFile
[1669]2727    from Numeric import Float, zeros, sum, reshape, equal
[1666]2728
2729    root = basename_in
2730    inname = root + '.dem'
2731
2732    #Open existing netcdf file to read
2733    infile = NetCDFFile(inname, 'r')
2734    if verbose: print 'Reading DEM from %s' %inname
2735
2736    #Read metadata
2737    ncols = infile.ncols[0]
2738    nrows = infile.nrows[0]
2739    xllcorner = infile.xllcorner[0]
2740    yllcorner = infile.yllcorner[0]
2741    cellsize = infile.cellsize[0]
2742    NODATA_value = infile.NODATA_value[0]
2743    zone = infile.zone[0]
2744    false_easting = infile.false_easting[0]
2745    false_northing = infile.false_northing[0]
2746    projection = infile.projection
2747    datum = infile.datum
2748    units = infile.units
2749
2750    dem_elevation = infile.variables['elevation']
2751
2752    #Get output file name
2753    if basename_out == None:
2754        outname = root + '_' + repr(cellsize_new) + '.dem'
2755    else:
2756        outname = basename_out + '.dem'
2757
2758    if verbose: print 'Write decimated NetCDF file to %s' %outname
2759
2760    #Determine some dimensions for decimated grid
2761    (nrows_stencil, ncols_stencil) = stencil.shape
2762    x_offset = ncols_stencil / 2
2763    y_offset = nrows_stencil / 2
[1667]2764    cellsize_ratio = int(cellsize_new / cellsize)
[1666]2765    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
2766    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
2767
2768    #Open netcdf file for output
2769    outfile = NetCDFFile(outname, 'w')
2770
2771    #Create new file
2772    outfile.institution = 'Geoscience Australia'
2773    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
2774                           'of spatial point data'
2775    #Georeferencing
2776    outfile.zone = zone
2777    outfile.projection = projection
2778    outfile.datum = datum
2779    outfile.units = units
2780
2781    outfile.cellsize = cellsize_new
2782    outfile.NODATA_value = NODATA_value
2783    outfile.false_easting = false_easting
2784    outfile.false_northing = false_northing
2785
2786    outfile.xllcorner = xllcorner + (x_offset * cellsize)
2787    outfile.yllcorner = yllcorner + (y_offset * cellsize)
2788    outfile.ncols = ncols_new
2789    outfile.nrows = nrows_new
2790
2791    # dimension definition
2792    outfile.createDimension('number_of_points', nrows_new*ncols_new)
2793
2794    # variable definition
2795    outfile.createVariable('elevation', Float, ('number_of_points',))
2796
2797    # Get handle to the variable
2798    elevation = outfile.variables['elevation']
2799
2800    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
2801
2802    #Store data
2803    global_index = 0
2804    for i in range(nrows_new):
[1667]2805        if verbose: print 'Processing row %d of %d' %(i, nrows_new)
[1666]2806        lower_index = global_index
2807        telev =  zeros(ncols_new, Float)
2808        local_index = 0
2809        trow = i * cellsize_ratio
2810
2811        for j in range(ncols_new):
2812            tcol = j * cellsize_ratio
[1669]2813            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
2814
2815            #if dem contains 1 or more NODATA_values set value in
2816            #decimated dem to NODATA_value, else compute decimated
2817            #value using stencil
2818            if sum(sum(equal(tmp, NODATA_value))) > 0:
2819                telev[local_index] = NODATA_value
2820            else:
2821                telev[local_index] = sum(sum(tmp * stencil))
2822
[1666]2823            global_index += 1
2824            local_index += 1
2825
2826        upper_index = global_index
2827
2828        elevation[lower_index:upper_index] = telev
2829
2830    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
2831
2832    infile.close()
2833    outfile.close()
2834
[1866]2835
2836   
2837def sww2asc_obsolete(basename_in, basename_out = None,
2838            quantity = None,
2839            timestep = None,
2840            reduction = None,
2841            cellsize = 10,
2842            verbose = False,
2843            origin = None):
2844    """Read SWW file and convert to Digitial Elevation model format (.asc)
2845
2846    Example:
2847
2848    ncols         3121
2849    nrows         1800
2850    xllcorner     722000
2851    yllcorner     5893000
2852    cellsize      25
2853    NODATA_value  -9999
2854    138.3698 137.4194 136.5062 135.5558 ..........
2855
2856    Also write accompanying file with same basename_in but extension .prj
2857    used to fix the UTM zone, datum, false northings and eastings.
2858
2859    The prj format is assumed to be as
2860
2861    Projection    UTM
2862    Zone          56
2863    Datum         WGS84
2864    Zunits        NO
2865    Units         METERS
2866    Spheroid      WGS84
2867    Xshift        0.0000000000
2868    Yshift        10000000.0000000000
2869    Parameters
2870
2871
2872    if quantity is given, out values from quantity otherwise default to
2873    elevation
2874
2875    if timestep (an index) is given, output quantity at that timestep
2876
2877    if reduction is given use that to reduce quantity over all timesteps.
2878
2879    """
2880    from Numeric import array, Float, concatenate, NewAxis, zeros,\
2881         sometrue
2882
2883
2884    #FIXME: Should be variable
2885    datum = 'WGS84'
2886    false_easting = 500000
2887    false_northing = 10000000
2888
2889    if quantity is None:
2890        quantity = 'elevation'
2891
2892    if reduction is None:
2893        reduction = max
2894
2895    if basename_out is None:
2896        basename_out = basename_in + '_%s' %quantity
2897
2898    swwfile = basename_in + '.sww'
2899    ascfile = basename_out + '.asc'
2900    prjfile = basename_out + '.prj'
2901
2902
2903    if verbose: print 'Reading from %s' %swwfile
2904    #Read sww file
2905    from Scientific.IO.NetCDF import NetCDFFile
2906    fid = NetCDFFile(swwfile)
2907
2908    #Get extent and reference
2909    x = fid.variables['x'][:]
2910    y = fid.variables['y'][:]
2911    volumes = fid.variables['volumes'][:]
2912
2913    ymin = min(y); ymax = max(y)
2914    xmin = min(x); xmax = max(x)
2915
2916    number_of_timesteps = fid.dimensions['number_of_timesteps']
2917    number_of_points = fid.dimensions['number_of_points']
2918    if origin is None:
2919
2920        #Get geo_reference
2921        #sww files don't have to have a geo_ref
2922        try:
2923            geo_reference = Geo_reference(NetCDFObject=fid)
2924        except AttributeError, e:
2925            geo_reference = Geo_reference() #Default georef object
2926
2927        xllcorner = geo_reference.get_xllcorner()
2928        yllcorner = geo_reference.get_yllcorner()
2929        zone = geo_reference.get_zone()
2930    else:
2931        zone = origin[0]
2932        xllcorner = origin[1]
2933        yllcorner = origin[2]
2934
2935
2936    #Get quantity and reduce if applicable
2937    if verbose: print 'Reading quantity %s' %quantity
2938
2939    if quantity.lower() == 'depth':
2940        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
2941    else:
2942        q = fid.variables[quantity][:]
2943
2944
2945    if len(q.shape) == 2:
2946        if verbose: print 'Reducing quantity %s' %quantity
2947        q_reduced = zeros( number_of_points, Float )
2948
2949        for k in range(number_of_points):
2950            q_reduced[k] = reduction( q[:,k] )
2951
2952        q = q_reduced
2953
2954    #Now q has dimension: number_of_points
2955
2956    #Create grid and update xll/yll corner and x,y
2957    if verbose: print 'Creating grid'
2958    ncols = int((xmax-xmin)/cellsize)+1
2959    nrows = int((ymax-ymin)/cellsize)+1
2960
2961    newxllcorner = xmin+xllcorner
2962    newyllcorner = ymin+yllcorner
2963
2964    x = x+xllcorner-newxllcorner
2965    y = y+yllcorner-newyllcorner
2966
2967    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2968    assert len(vertex_points.shape) == 2
2969
2970
2971    from Numeric import zeros, Float
2972    grid_points = zeros ( (ncols*nrows, 2), Float )
2973
2974
2975    for i in xrange(nrows):
2976        yg = i*cellsize
2977        for j in xrange(ncols):
2978            xg = j*cellsize
2979            k = i*ncols + j
2980
2981            grid_points[k,0] = xg
2982            grid_points[k,1] = yg
2983
2984    #Interpolate
2985    from least_squares import Interpolation
2986    from util import inside_polygon
2987
2988    #FIXME: This should be done with precrop = True, otherwise it'll
2989    #take forever. With expand_search  set to False, some grid points might
2990    #miss out....
2991
2992    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
2993                           precrop = False, expand_search = False,
2994                           verbose = verbose)
2995
2996    #Interpolate using quantity values
2997    if verbose: print 'Interpolating'
2998    grid_values = interp.interpolate(q).flat
2999
3000    #Write
3001    #Write prj file
3002    if verbose: print 'Writing %s' %prjfile
3003    prjid = open(prjfile, 'w')
3004    prjid.write('Projection    %s\n' %'UTM')
3005    prjid.write('Zone          %d\n' %zone)
3006    prjid.write('Datum         %s\n' %datum)
3007    prjid.write('Zunits        NO\n')
3008    prjid.write('Units         METERS\n')
3009    prjid.write('Spheroid      %s\n' %datum)
3010    prjid.write('Xshift        %d\n' %false_easting)
3011    prjid.write('Yshift        %d\n' %false_northing)
3012    prjid.write('Parameters\n')
3013    prjid.close()
3014
3015
3016   
3017    if verbose: print 'Writing %s' %ascfile
3018    NODATA_value = -9999
3019 
3020    ascid = open(ascfile, 'w')
3021
3022    ascid.write('ncols         %d\n' %ncols)
3023    ascid.write('nrows         %d\n' %nrows)
3024    ascid.write('xllcorner     %d\n' %newxllcorner)
3025    ascid.write('yllcorner     %d\n' %newyllcorner)
3026    ascid.write('cellsize      %f\n' %cellsize)
3027    ascid.write('NODATA_value  %d\n' %NODATA_value)
3028
3029
3030    #Get bounding polygon from mesh
3031    P = interp.mesh.get_boundary_polygon()
3032    inside_indices = inside_polygon(grid_points, P)
3033
3034    for i in range(nrows):
3035        if verbose and i%((nrows+10)/10)==0:
3036            print 'Doing row %d of %d' %(i, nrows)
3037
3038        for j in range(ncols):
3039            index = (nrows-i-1)*ncols+j
3040
3041            if sometrue(inside_indices == index):
3042                ascid.write('%f ' %grid_values[index])
3043            else:
3044                ascid.write('%d ' %NODATA_value)
3045
3046        ascid.write('\n')
3047
3048    #Close
3049    ascid.close()
3050    fid.close()
3051
3052def sww2ers_obsolete(basename_in, basename_out = None,
3053            quantity = None,
3054            timestep = None,
3055            reduction = None,
3056            cellsize = 10,
3057            verbose = False,
3058            origin = None,
3059            datum = 'WGS84'):
3060   
3061    """Read SWW file and convert to Digitial Elevation model format (.asc)
3062
3063    Example:
3064
3065    ncols         3121
3066    nrows         1800
3067    xllcorner     722000
3068    yllcorner     5893000
3069    cellsize      25
3070    NODATA_value  -9999
3071    138.3698 137.4194 136.5062 135.5558 ..........
3072
3073    Also write accompanying file with same basename_in but extension .prj
3074    used to fix the UTM zone, datum, false northings and eastings.
3075
3076    The prj format is assumed to be as
3077
3078    Projection    UTM
3079    Zone          56
3080    Datum         WGS84
3081    Zunits        NO
3082    Units         METERS
3083    Spheroid      WGS84
3084    Xshift        0.0000000000
3085    Yshift        10000000.0000000000
3086    Parameters
3087
3088
3089    if quantity is given, out values from quantity otherwise default to
3090    elevation
3091
3092    if timestep (an index) is given, output quantity at that timestep
3093
3094    if reduction is given use that to reduce quantity over all timesteps.
3095
3096    """
3097    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
3098    import ermapper_grids
3099
3100    header = {}
3101    false_easting = 500000
3102    false_northing = 10000000
3103    NODATA_value = 0
3104
3105    if quantity is None:
3106        quantity = 'elevation'
3107
3108    if reduction is None:
3109        reduction = max
3110
3111    if basename_out is None:
3112        basename_out = basename_in + '_%s' %quantity
3113 
3114    swwfile = basename_in + '.sww'
3115    # Note the use of a .ers extension is optional (write_ermapper_grid will
3116    # deal with either option
3117    ersfile = basename_out
3118##    prjfile = basename_out + '.prj'
3119
3120
3121    if verbose: print 'Reading from %s' %swwfile
3122    #Read sww file
3123    from Scientific.IO.NetCDF import NetCDFFile
3124    fid = NetCDFFile(swwfile)
3125
3126    #Get extent and reference
3127    x = fid.variables['x'][:]
3128    y = fid.variables['y'][:]
3129    volumes = fid.variables['volumes'][:]
3130
3131    ymin = min(y); ymax = max(y)
3132    xmin = min(x); xmax = max(x)
3133
3134    number_of_timesteps = fid.dimensions['number_of_timesteps']
3135    number_of_points = fid.dimensions['number_of_points']
3136    if origin is None:
3137
3138        #Get geo_reference
3139        #sww files don't have to have a geo_ref
3140        try:
3141            geo_reference = Geo_reference(NetCDFObject=fid)
3142        except AttributeError, e:
3143            geo_reference = Geo_reference() #Default georef object
3144
3145        xllcorner = geo_reference.get_xllcorner()
3146        yllcorner = geo_reference.get_yllcorner()
3147        zone = geo_reference.get_zone()
3148    else:
3149        zone = origin[0]
3150        xllcorner = origin[1]
3151        yllcorner = origin[2]
3152
3153
3154    #Get quantity and reduce if applicable
3155    if verbose: print 'Reading quantity %s' %quantity
3156
3157    if quantity.lower() == 'depth':
3158        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
3159    else:
3160        q = fid.variables[quantity][:]
3161
3162
3163    if len(q.shape) == 2:
3164        if verbose: print 'Reducing quantity %s' %quantity
3165        q_reduced = zeros( number_of_points, Float )
3166
3167        for k in range(number_of_points):
3168            q_reduced[k] = reduction( q[:,k] )
3169
3170        q = q_reduced
3171
3172    #Now q has dimension: number_of_points
3173
3174    #Create grid and update xll/yll corner and x,y
3175    if verbose: print 'Creating grid'
3176    ncols = int((xmax-xmin)/cellsize)+1
3177    nrows = int((ymax-ymin)/cellsize)+1
3178
3179    newxllcorner = xmin+xllcorner
3180    newyllcorner = ymin+yllcorner
3181
3182    x = x+xllcorner-newxllcorner
3183    y = y+yllcorner-newyllcorner
3184
3185    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
3186    assert len(vertex_points.shape) == 2
3187
3188
3189    from Numeric import zeros, Float
3190    grid_points = zeros ( (ncols*nrows, 2), Float )
3191
3192
3193    for i in xrange(nrows):
3194        yg = (nrows-i)*cellsize # this will flip the order of the y values
3195        for j in xrange(ncols):
3196            xg = j*cellsize
3197            k = i*ncols + j
3198
3199            grid_points[k,0] = xg
3200            grid_points[k,1] = yg
3201
3202    #Interpolate
3203    from least_squares import Interpolation
3204    from util import inside_polygon
3205
3206    #FIXME: This should be done with precrop = True (?), otherwise it'll
3207    #take forever. With expand_search  set to False, some grid points might
3208    #miss out....
3209
3210    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
3211                           precrop = False, expand_search = False,
3212                           verbose = verbose)
3213
3214    #Interpolate using quantity values
3215    if verbose: print 'Interpolating'
3216    grid_values = interp.interpolate(q).flat
3217    grid_values = reshape(grid_values,(nrows, ncols))
3218
3219
3220    # setup header information
3221    header['datum'] = '"' + datum + '"'
3222    # FIXME The use of hardwired UTM and zone number needs to be made optional
3223    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
3224    header['projection'] = '"UTM-' + str(zone) + '"'
3225    header['coordinatetype'] = 'EN'
3226    if header['coordinatetype'] == 'LL':
3227        header['longitude'] = str(newxllcorner)
3228        header['latitude'] = str(newyllcorner)
3229    elif header['coordinatetype'] == 'EN':
3230        header['eastings'] = str(newxllcorner)
3231        header['northings'] = str(newyllcorner)
3232    header['nullcellvalue'] = str(NODATA_value)
3233    header['xdimension'] = str(cellsize)
3234    header['ydimension'] = str(cellsize)
3235    header['value'] = '"' + quantity + '"'
3236
3237
3238    #Write
3239    if verbose: print 'Writing %s' %ersfile
3240    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
3241
3242    fid.close()   
3243
3244##    ascid = open(ascfile, 'w')
3245##
3246##    ascid.write('ncols         %d\n' %ncols)
3247##    ascid.write('nrows         %d\n' %nrows)
3248##    ascid.write('xllcorner     %d\n' %newxllcorner)
3249##    ascid.write('yllcorner     %d\n' %newyllcorner)
3250##    ascid.write('cellsize      %f\n' %cellsize)
3251##    ascid.write('NODATA_value  %d\n' %NODATA_value)
3252##
3253##
3254##    #Get bounding polygon from mesh
3255##    P = interp.mesh.get_boundary_polygon()
3256##    inside_indices = inside_polygon(grid_points, P)
3257##
3258##    for i in range(nrows):
3259##        if verbose and i%((nrows+10)/10)==0:
3260##            print 'Doing row %d of %d' %(i, nrows)
3261##
3262##        for j in range(ncols):
3263##            index = (nrows-i-1)*ncols+j
3264##
3265##            if sometrue(inside_indices == index):
3266##                ascid.write('%f ' %grid_values[index])
3267##            else:
3268##                ascid.write('%d ' %NODATA_value)
3269##
3270##        ascid.write('\n')
3271##
3272##    #Close
3273##    ascid.close()
3274##    fid.close()
3275
Note: See TracBrowser for help on using the repository browser.