source: development/pyvolution-1d/data_manager.py @ 3322

Last change on this file since 3322 was 2792, checked in by jakeman, 19 years ago

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