source: anuga_core/source/anuga/shallow_water/data_manager.py @ 4776

Last change on this file since 4776 was 4776, checked in by ole, 16 years ago

Cosmetics according to style guide.

File size: 188.4 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.csv: 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 anuga.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 abstract_2d_finite_volumes
48
49"""
50
51import exceptions
52class TitleValueError(exceptions.Exception): pass
53class DataMissingValuesError(exceptions.Exception): pass
54class DataFileNotOpenError(exceptions.Exception): pass
55class DataTimeError(exceptions.Exception): pass
56class DataDomainError(exceptions.Exception): pass
57class NewQuantity(exceptions.Exception): pass
58
59
60
61import csv
62import os, sys
63import shutil
64from struct import unpack
65import array as p_array
66#import time, os
67from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
68
69
70from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
71     searchsorted, zeros, allclose, around, reshape, transpose, sort, \
72     NewAxis, ArrayType, compress, take, arange, argmax, alltrue,shape,Float
73
74import string
75
76from Scientific.IO.NetCDF import NetCDFFile
77#from shutil import copy
78from os.path import exists, basename, join
79from os import getcwd
80
81
82from anuga.coordinate_transforms.redfearn import redfearn, \
83     convert_from_latlon_to_utm
84from anuga.coordinate_transforms.geo_reference import Geo_reference, \
85     write_NetCDF_georeference, ensure_geo_reference
86from anuga.geospatial_data.geospatial_data import Geospatial_data,\
87     ensure_absolute
88from anuga.config import minimum_storable_height as default_minimum_storable_height
89from anuga.config import max_float
90from anuga.utilities.numerical_tools import ensure_numeric,  mean
91from anuga.caching.caching import myhash
92from anuga.utilities.anuga_exceptions import ANUGAError
93from anuga.shallow_water import Domain
94from anuga.abstract_2d_finite_volumes.pmesh2domain import \
95     pmesh_to_domain_instance
96from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
97     remove_lone_verts, sww2timeseries, get_centroid_values
98from anuga.load_mesh.loadASCII import export_mesh_file
99# formula mappings
100
101quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
102                    'depth':'stage-elevation',
103                    'speed': \
104 '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
105
106
107   
108def make_filename(s):
109    """Transform argument string into a Sexsuitable filename
110    """
111
112    s = s.strip()
113    s = s.replace(' ', '_')
114    s = s.replace('(', '')
115    s = s.replace(')', '')
116    s = s.replace('__', '_')
117
118    return s
119
120
121def check_dir(path, verbose=None):
122    """Check that specified path exists.
123    If path does not exist it will be created if possible
124
125    USAGE:
126       checkdir(path, verbose):
127
128    ARGUMENTS:
129        path -- Directory
130        verbose -- Flag verbose output (default: None)
131
132    RETURN VALUE:
133        Verified path including trailing separator
134
135    """
136
137    import os.path
138
139    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
140        unix = 0
141    else:
142        unix = 1
143
144
145    if path[-1] != os.sep:
146        path = path + os.sep  # Add separator for directories
147
148    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
149    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
150        try:
151            exitcode=os.mkdir(path)
152
153            # Change access rights if possible
154            #
155            if unix:
156                exitcode=os.system('chmod 775 '+path)
157            else:
158                pass  # FIXME: What about acces rights under Windows?
159
160            if verbose: print 'MESSAGE: Directory', path, 'created.'
161
162        except:
163            print 'WARNING: Directory', path, 'could not be created.'
164            if unix:
165                path = '/tmp/'
166            else:
167                path = 'C:'
168
169            print 'Using directory %s instead' %path
170
171    return(path)
172
173
174
175def del_dir(path):
176    """Recursively delete directory path and all its contents
177    """
178
179    import os
180
181    if os.path.isdir(path):
182        for file in os.listdir(path):
183            X = os.path.join(path, file)
184
185
186            if os.path.isdir(X) and not os.path.islink(X):
187                del_dir(X)
188            else:
189                try:
190                    os.remove(X)
191                except:
192                    print "Could not remove file %s" %X
193
194        os.rmdir(path)
195       
196       
197# ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007   
198def rmgeneric(path, __func__,verbose=False):
199    ERROR_STR= """Error removing %(path)s, %(error)s """
200
201    try:
202        __func__(path)
203        if verbose: print 'Removed ', path
204    except OSError, (errno, strerror):
205        print ERROR_STR % {'path' : path, 'error': strerror }
206           
207def removeall(path,verbose=False):
208
209    if not os.path.isdir(path):
210        return
211   
212    files=os.listdir(path)
213
214    for x in files:
215        fullpath=os.path.join(path, x)
216        if os.path.isfile(fullpath):
217            f=os.remove
218            rmgeneric(fullpath, f)
219        elif os.path.isdir(fullpath):
220            removeall(fullpath)
221            f=os.rmdir
222            rmgeneric(fullpath, f,verbose)
223
224
225
226def create_filename(datadir, filename, format, size=None, time=None):
227
228    import os
229    #from anuga.config import data_dir
230
231    FN = check_dir(datadir) + filename
232
233    if size is not None:
234        FN += '_size%d' %size
235
236    if time is not None:
237        FN += '_time%.2f' %time
238
239    FN += '.' + format
240    return FN
241
242
243def get_files(datadir, filename, format, size):
244    """Get all file (names) with given name, size and format
245    """
246
247    import glob
248
249    import os
250    #from anuga.config import data_dir
251
252    dir = check_dir(datadir)
253
254    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
255    return glob.glob(pattern)
256
257
258
259#Generic class for storing output to e.g. visualisation or checkpointing
260class Data_format:
261    """Generic interface to data formats
262    """
263
264
265    def __init__(self, domain, extension, mode = 'w'):
266        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
267                                        '''   'w' (write)'''+\
268                                        '''   'r' (read)''' +\
269                                        '''   'a' (append)'''
270
271        #Create filename
272        self.filename = create_filename(domain.get_datadir(),
273                                        domain.get_name(), extension)
274
275        #print 'F', self.filename
276        self.timestep = 0
277        self.domain = domain
278       
279
280
281        # Exclude ghosts in case this is a parallel domain
282        self.number_of_nodes = domain.number_of_full_nodes       
283        self.number_of_volumes = domain.number_of_full_triangles
284        #self.number_of_volumes = len(domain)       
285
286
287
288
289        #FIXME: Should we have a general set_precision function?
290
291
292
293#Class for storing output to e.g. visualisation
294class Data_format_sww(Data_format):
295    """Interface to native NetCDF format (.sww) for storing model output
296
297    There are two kinds of data
298
299    1: Constant data: Vertex coordinates and field values. Stored once
300    2: Variable data: Conserved quantities. Stored once per timestep.
301
302    All data is assumed to reside at vertex locations.
303    """
304
305
306    def __init__(self, domain, mode = 'w',\
307                 max_size = 2000000000,
308                 recursion = False):
309        from Scientific.IO.NetCDF import NetCDFFile
310        from Numeric import Int, Float, Float32
311
312        self.precision = Float32 #Use single precision for quantities
313        if hasattr(domain, 'max_size'):
314            self.max_size = domain.max_size #file size max is 2Gig
315        else:
316            self.max_size = max_size
317        self.recursion = recursion
318        self.mode = mode
319
320        Data_format.__init__(self, domain, 'sww', mode)
321
322        if hasattr(domain, 'minimum_storable_height'):
323            self.minimum_storable_height = domain.minimum_storable_height
324        else:
325            self.minimum_storable_height = default_minimum_storable_height
326
327        # NetCDF file definition
328        fid = NetCDFFile(self.filename, mode)
329
330        if mode == 'w':
331            description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
332            self.writer = Write_sww()
333            self.writer.store_header(fid,
334                                     domain.starttime,
335                                     self.number_of_volumes,
336                                     self.domain.number_of_full_nodes,
337                                     description=description,
338                                     smoothing=domain.smooth,
339                                     order=domain.default_order)
340
341            # Extra optional information
342            if hasattr(domain, 'texture'):
343                fid.texture = domain.texture
344
345            if domain.quantities_to_be_monitored is not None:
346                fid.createDimension('singleton', 1)
347                fid.createDimension('two', 2)               
348
349                poly = domain.monitor_polygon
350                if poly is not None:
351                    N = len(poly)
352                    fid.createDimension('polygon_length', N)
353                    fid.createVariable('extrema.polygon',
354                                       self.precision,
355                                       ('polygon_length',
356                                        'two'))
357                    fid.variables['extrema.polygon'][:] = poly                                   
358
359                   
360                interval = domain.monitor_time_interval
361                if interval is not None:
362                    fid.createVariable('extrema.time_interval',
363                                       self.precision,
364                                       ('two',))
365                    fid.variables['extrema.time_interval'][:] = interval
366
367               
368                for q in domain.quantities_to_be_monitored:
369                    #print 'doing', q
370                    fid.createVariable(q+'.extrema', self.precision,
371                                       ('numbers_in_range',))
372                    fid.createVariable(q+'.min_location', self.precision,
373                                       ('numbers_in_range',))
374                    fid.createVariable(q+'.max_location', self.precision,
375                                       ('numbers_in_range',))
376                    fid.createVariable(q+'.min_time', self.precision,
377                                       ('singleton',))
378                    fid.createVariable(q+'.max_time', self.precision,
379                                       ('singleton',))
380
381                   
382        fid.close()
383
384
385    def store_connectivity(self):
386        """Specialisation of store_connectivity for net CDF format
387
388        Writes x,y,z coordinates of triangles constituting
389        the bed elevation.
390        """
391
392        from Scientific.IO.NetCDF import NetCDFFile
393
394        from Numeric import concatenate, Int
395
396        domain = self.domain
397
398        #Get NetCDF
399        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
400
401        # Get the variables
402        x = fid.variables['x']
403        y = fid.variables['y']
404        z = fid.variables['elevation']
405
406        volumes = fid.variables['volumes']
407
408        # Get X, Y and bed elevation Z
409        Q = domain.quantities['elevation']
410        X,Y,Z,V = Q.get_vertex_values(xy=True,
411                                      precision=self.precision)
412
413        #
414        points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
415        self.writer.store_triangulation(fid,
416                                        points,
417                                        V.astype(volumes.typecode()),
418                                        Z,
419                                        points_georeference= \
420                                        domain.geo_reference)
421
422        # Close
423        fid.close()
424
425
426    def store_timestep(self, names):
427        """Store time and named quantities to file
428        """
429        from Scientific.IO.NetCDF import NetCDFFile
430        import types
431        from time import sleep
432        from os import stat
433
434        from Numeric import choose
435       
436        # Get NetCDF       
437        retries = 0
438        file_open = False
439        while not file_open and retries < 10:
440            try:
441                fid = NetCDFFile(self.filename, 'a') # Open existing file
442            except IOError:
443                # This could happen if someone was reading the file.
444                # In that case, wait a while and try again
445                msg = 'Warning (store_timestep): File %s could not be opened'\
446                      %self.filename
447                msg += ' - trying step %s again' %self.domain.time
448                print msg
449                retries += 1
450                sleep(1)
451            else:
452                file_open = True
453
454        if not file_open:
455            msg = 'File %s could not be opened for append' %self.filename
456            raise DataFileNotOpenError, msg
457
458
459
460        # Check to see if the file is already too big:
461        time = fid.variables['time']
462        i = len(time)+1
463        file_size = stat(self.filename)[6]
464        file_size_increase =  file_size/i
465        if file_size + file_size_increase > self.max_size*(2**self.recursion):
466            # In order to get the file name and start time correct,
467            # I change the domain.filename and domain.starttime.
468            # This is the only way to do this without changing
469            # other modules (I think).
470
471            # Write a filename addon that won't break swollens reader
472            # (10.sww is bad)
473            filename_ext = '_time_%s'%self.domain.time
474            filename_ext = filename_ext.replace('.', '_')
475           
476            # Remember the old filename, then give domain a
477            # name with the extension
478            old_domain_filename = self.domain.get_name()
479            if not self.recursion:
480                self.domain.set_name(old_domain_filename+filename_ext)
481
482
483            # Change the domain starttime to the current time
484            old_domain_starttime = self.domain.starttime
485            self.domain.starttime = self.domain.time
486
487            # Build a new data_structure.
488            next_data_structure=\
489                Data_format_sww(self.domain, mode=self.mode,\
490                                max_size = self.max_size,\
491                                recursion = self.recursion+1)
492            if not self.recursion:
493                print '    file_size = %s'%file_size
494                print '    saving file to %s'%next_data_structure.filename
495            #set up the new data_structure
496            self.domain.writer = next_data_structure
497
498            #FIXME - could be cleaner to use domain.store_timestep etc.
499            next_data_structure.store_connectivity()
500            next_data_structure.store_timestep(names)
501            fid.sync()
502            fid.close()
503
504            #restore the old starttime and filename
505            self.domain.starttime = old_domain_starttime
506            self.domain.set_name(old_domain_filename)           
507        else:
508            self.recursion = False
509            domain = self.domain
510
511            # Get the variables
512            time = fid.variables['time']
513            stage = fid.variables['stage']
514            xmomentum = fid.variables['xmomentum']
515            ymomentum = fid.variables['ymomentum']
516            i = len(time)
517            if type(names) not in [types.ListType, types.TupleType]:
518                names = [names]
519
520            if 'stage' in names and 'xmomentum' in names and \
521               'ymomentum' in names:
522
523                # Get stage
524                Q = domain.quantities['stage']
525                A,_ = Q.get_vertex_values(xy = False,
526                                          precision = self.precision)               
527                z = fid.variables['elevation']
528                stage = choose(A-z[:] >= self.minimum_storable_height,
529                           (z[:], A))
530               
531                # Get xmomentum
532                Q = domain.quantities['xmomentum']
533                xmomentum, _ = Q.get_vertex_values(xy = False,
534                                          precision = self.precision)
535               
536                # Get ymomentum
537                Q = domain.quantities['ymomentum']
538                ymomentum, _ = Q.get_vertex_values(xy = False,
539                                          precision = self.precision)
540               
541                # Write quantities to NetCDF
542                self.writer.store_quantities(fid, 
543                                             time=self.domain.time,
544                                             precision=self.precision,
545                                             stage=stage,
546                                             xmomentum=xmomentum,
547                                             ymomentum=ymomentum)
548            else:
549                # This is producing a sww that is not standard.
550                # Store time
551                time[i] = self.domain.time
552               
553                for name in names:
554                    # Get quantity
555                    Q = domain.quantities[name]
556                    A,V = Q.get_vertex_values(xy = False,
557                                              precision = self.precision)
558
559                    # FIXME: Make this general (see below)
560                    if name == 'stage':
561                        z = fid.variables['elevation']
562                        A = choose(A-z[:] >= self.minimum_storable_height,
563                                   (z[:], A))
564                        stage[i,:] = A.astype(self.precision)
565                    elif name == 'xmomentum':
566                        xmomentum[i,:] = A.astype(self.precision)
567                    elif name == 'ymomentum':
568                        ymomentum[i,:] = A.astype(self.precision)
569
570                   #As in....
571                   #eval( name + '[i,:] = A.astype(self.precision)' )
572                   #FIXME (Ole): But we need a UNIT test for that before
573                   # refactoring
574
575
576
577            # Update extrema if requested
578            domain = self.domain
579            if domain.quantities_to_be_monitored is not None:
580                for q, info in domain.quantities_to_be_monitored.items():
581
582                    if info['min'] is not None:
583                        fid.variables[q + '.extrema'][0] = info['min']
584                        fid.variables[q + '.min_location'][:] =\
585                                        info['min_location']
586                        fid.variables[q + '.min_time'][0] = info['min_time']
587                       
588                    if info['max'] is not None:
589                        fid.variables[q + '.extrema'][1] = info['max']
590                        fid.variables[q + '.max_location'][:] =\
591                                        info['max_location']
592                        fid.variables[q + '.max_time'][0] = info['max_time']
593
594           
595
596            #Flush and close
597            fid.sync()
598            fid.close()
599
600
601
602# Class for handling checkpoints data
603class Data_format_cpt(Data_format):
604    """Interface to native NetCDF format (.cpt)
605    """
606
607
608    def __init__(self, domain, mode = 'w'):
609        from Scientific.IO.NetCDF import NetCDFFile
610        from Numeric import Int, Float, Float
611
612        self.precision = Float #Use full precision
613
614        Data_format.__init__(self, domain, 'sww', mode)
615
616
617        # NetCDF file definition
618        fid = NetCDFFile(self.filename, mode)
619
620        if mode == 'w':
621            #Create new file
622            fid.institution = 'Geoscience Australia'
623            fid.description = 'Checkpoint data'
624            #fid.smooth = domain.smooth
625            fid.order = domain.default_order
626
627            # dimension definitions
628            fid.createDimension('number_of_volumes', self.number_of_volumes)
629            fid.createDimension('number_of_vertices', 3)
630
631            #Store info at all vertices (no smoothing)
632            fid.createDimension('number_of_points', 3*self.number_of_volumes)
633            fid.createDimension('number_of_timesteps', None) #extensible
634
635            # variable definitions
636
637            #Mesh
638            fid.createVariable('x', self.precision, ('number_of_points',))
639            fid.createVariable('y', self.precision, ('number_of_points',))
640
641
642            fid.createVariable('volumes', Int, ('number_of_volumes',
643                                                'number_of_vertices'))
644
645            fid.createVariable('time', self.precision,
646                               ('number_of_timesteps',))
647
648            #Allocate space for all quantities
649            for name in domain.quantities.keys():
650                fid.createVariable(name, self.precision,
651                                   ('number_of_timesteps',
652                                    'number_of_points'))
653
654        #Close
655        fid.close()
656
657
658    def store_checkpoint(self):
659        """
660        Write x,y coordinates of triangles.
661        Write connectivity (
662        constituting
663        the bed elevation.
664        """
665
666        from Scientific.IO.NetCDF import NetCDFFile
667
668        from Numeric import concatenate
669
670        domain = self.domain
671
672        #Get NetCDF
673        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
674
675        # Get the variables
676        x = fid.variables['x']
677        y = fid.variables['y']
678
679        volumes = fid.variables['volumes']
680
681        # Get X, Y and bed elevation Z
682        Q = domain.quantities['elevation']
683        X,Y,Z,V = Q.get_vertex_values(xy=True,
684                      precision = self.precision)
685
686
687
688        x[:] = X.astype(self.precision)
689        y[:] = Y.astype(self.precision)
690        z[:] = Z.astype(self.precision)
691
692        volumes[:] = V
693
694        #Close
695        fid.close()
696
697
698    def store_timestep(self, name):
699        """Store time and named quantity to file
700        """
701        from Scientific.IO.NetCDF import NetCDFFile
702        from time import sleep
703
704        #Get NetCDF
705        retries = 0
706        file_open = False
707        while not file_open and retries < 10:
708            try:
709                fid = NetCDFFile(self.filename, 'a') #Open existing file
710            except IOError:
711                #This could happen if someone was reading the file.
712                #In that case, wait a while and try again
713                msg = 'Warning (store_timestep): File %s could not be opened'\
714                  %self.filename
715                msg += ' - trying again'
716                print msg
717                retries += 1
718                sleep(1)
719            else:
720                file_open = True
721
722        if not file_open:
723            msg = 'File %s could not be opened for append' %self.filename
724            raise DataFileNotOPenError, msg
725
726
727        domain = self.domain
728
729        # Get the variables
730        time = fid.variables['time']
731        stage = fid.variables['stage']
732        i = len(time)
733
734        #Store stage
735        time[i] = self.domain.time
736
737        # Get quantity
738        Q = domain.quantities[name]
739        A,V = Q.get_vertex_values(xy=False,
740                                  precision = self.precision)
741
742        stage[i,:] = A.astype(self.precision)
743
744        #Flush and close
745        fid.sync()
746        fid.close()
747
748
749#### NED is national exposure database (name changed to NEXIS)
750   
751LAT_TITLE = 'LATITUDE'
752LONG_TITLE = 'LONGITUDE'
753X_TITLE = 'x'
754Y_TITLE = 'y'
755class Exposure_csv:
756    def __init__(self,file_name, latitude_title=LAT_TITLE,
757                 longitude_title=LONG_TITLE, is_x_y_locations=None,
758                 x_title=X_TITLE, y_title=Y_TITLE,
759                 refine_polygon=None, title_check_list=None):
760        """
761        This class is for handling the exposure csv file.
762        It reads the file in and converts the lats and longs to a geospatial
763        data object.
764        Use the methods to read and write columns.
765
766        The format of the csv files it reads is;
767           The first row is a title row.
768           comma's are the delimiters
769           each column is a 'set' of data
770
771        Feel free to use/expand it to read other csv files.
772           
773           
774        It is not for adding and deleting rows
775       
776        Can geospatial handle string attributes? It's not made for them.
777        Currently it can't load and save string att's.
778
779        So just use geospatial to hold the x, y and georef? Bad, since
780        different att's are in diferent structures.  Not so bad, the info
781        to write if the .csv file is saved is in attribute_dic
782
783        The location info is in the geospatial attribute.
784       
785       
786        """
787        self._file_name = file_name
788        self._geospatial = None #
789
790        # self._attribute_dic is a dictionary.
791        #The keys are the column titles.
792        #The values are lists of column data
793       
794        # self._title_index_dic is a dictionary.
795        #The keys are the column titles.
796        #The values are the index positions of file columns.
797        self._attribute_dic, self._title_index_dic = \
798            csv2dict(self._file_name, title_check_list=title_check_list)
799        try:
800            #Have code here that handles caps or lower
801            lats = self._attribute_dic[latitude_title]
802            longs = self._attribute_dic[longitude_title]
803           
804        except KeyError:
805            # maybe a warning..
806            #Let's see if this works..
807            if False != is_x_y_locations:
808                is_x_y_locations = True
809            pass
810        else:
811            self._geospatial = Geospatial_data(latitudes = lats,
812                 longitudes = longs)
813
814        if is_x_y_locations is True:
815            if self._geospatial is not None:
816                pass #fixme throw an error
817            try:
818                xs = self._attribute_dic[x_title]
819                ys = self._attribute_dic[y_title]
820                points = [[float(i),float(j)] for i,j in map(None,xs,ys)]
821            except KeyError:
822                # maybe a warning..
823                msg = "Could not find location information."
824                raise TitleValueError, msg
825            else:
826                self._geospatial = Geospatial_data(data_points=points)
827           
828        # create a list of points that are in the refining_polygon
829        # described by a list of indexes representing the points
830
831    def __cmp__(self, other):
832        #print "self._attribute_dic",self._attribute_dic
833        #print "other._attribute_dic",other._attribute_dic
834        #print "self._title_index_dic", self._title_index_dic
835        #print "other._title_index_dic", other._title_index_dic
836       
837        #check that a is an instance of this class
838        if isinstance(self, type(other)):
839            result = cmp(self._attribute_dic, other._attribute_dic)
840            if result <>0:
841                return result
842            # The order of the columns is important. Therefore..
843            result = cmp(self._title_index_dic, other._title_index_dic)
844            if result <>0:
845                return result
846            for self_ls, other_ls in map(None,self._attribute_dic, \
847                    other._attribute_dic):
848                result = cmp(self._attribute_dic[self_ls],
849                             other._attribute_dic[other_ls])
850                if result <>0:
851                    return result
852            return 0
853        else:
854            return 1
855   
856
857    def get_column(self, column_name, use_refind_polygon=False):
858        """
859        Given a column name return a list of the column values
860
861        Note, the type of the values will be String!
862        do this to change a list of strings to a list of floats
863        time = [float(x) for x in time]
864       
865        Not implemented:
866        if use_refind_polygon is True, only return values in the
867        refined polygon
868        """
869        if not self._attribute_dic.has_key(column_name):
870            msg = 'Therer is  no column called %s!' %column_name
871            raise TitleValueError, msg
872        return self._attribute_dic[column_name]
873
874
875    def get_value(self, value_column_name,
876                  known_column_name,
877                  known_values,
878                  use_refind_polygon=False):
879        """
880        Do linear interpolation on the known_colum, using the known_value,
881        to return a value of the column_value_name.
882        """
883        pass
884
885
886    def get_location(self, use_refind_polygon=False):
887        """
888        Return a geospatial object which describes the
889        locations of the location file.
890
891        Note, if there is not location info, this returns None.
892       
893        Not implemented:
894        if use_refind_polygon is True, only return values in the
895        refined polygon
896        """
897        return self._geospatial
898
899    def set_column(self, column_name, column_values, overwrite=False):
900        """
901        Add a column to the 'end' (with the right most column being the end)
902        of the csv file.
903
904        Set overwrite to True if you want to overwrite a column.
905       
906        Note, in column_name white space is removed and case is not checked.
907        Precondition
908        The column_name and column_values cannot have comma's in it.
909        """
910       
911        value_row_count = \
912            len(self._attribute_dic[self._title_index_dic.keys()[0]])
913        if len(column_values) <> value_row_count: 
914            msg = 'The number of column values must equal the number of rows.'
915            raise DataMissingValuesError, msg
916       
917        if self._attribute_dic.has_key(column_name):
918            if not overwrite:
919                msg = 'Column name %s already in use!' %column_name
920                raise TitleValueError, msg
921        else:
922            # New title.  Add it to the title index.
923            self._title_index_dic[column_name] = len(self._title_index_dic)
924        self._attribute_dic[column_name] = column_values
925        #print "self._title_index_dic[column_name]",self._title_index_dic
926
927    def save(self, file_name=None):
928        """
929        Save the exposure csv file
930        """
931        if file_name is None:
932            file_name = self._file_name
933       
934        fd = open(file_name,'wb')
935        writer = csv.writer(fd)
936       
937        #Write the title to a cvs file
938        line = [None]* len(self._title_index_dic)
939        for title in self._title_index_dic.iterkeys():
940            line[self._title_index_dic[title]]= title
941        writer.writerow(line)
942       
943        # Write the values to a cvs file
944        value_row_count = \
945            len(self._attribute_dic[self._title_index_dic.keys()[0]])
946        for row_i in range(value_row_count):
947            line = [None]* len(self._title_index_dic)
948            for title in self._title_index_dic.iterkeys():
949                line[self._title_index_dic[title]]= \
950                     self._attribute_dic[title][row_i]
951            writer.writerow(line)
952
953
954def csv2dict(file_name, title_check_list=None):
955    """
956    Load in the csv as a dic, title as key and column info as value, .
957    Also, create a dic, title as key and column index as value,
958    to keep track of the column order.
959
960    Two dictionaries are returned.
961   
962    WARNING: Vaules are returned as strings.
963    do this to change a list of strings to a list of floats
964        time = [float(x) for x in time]
965
966       
967    """
968   
969    #
970    attribute_dic = {}
971    title_index_dic = {}
972    titles_stripped = [] # list of titles
973    reader = csv.reader(file(file_name))
974
975    # Read in and manipulate the title info
976    titles = reader.next()
977    for i,title in enumerate(titles):
978        titles_stripped.append(title.strip())
979        title_index_dic[title.strip()] = i
980    title_count = len(titles_stripped)       
981    #print "title_index_dic",title_index_dic
982    if title_check_list is not None:
983        for title_check in title_check_list:
984            #msg = "Reading error.  This row is not present ", title_check
985            #assert title_index_dic.has_key(title_check), msg
986            if not title_index_dic.has_key(title_check):
987                #reader.close()
988                msg = "Reading error.  This row is not present ", \
989                      title_check                     
990                raise IOError, msg
991               
992   
993   
994    #create a dic of colum values, indexed by column title
995    for line in reader:
996        if len(line) <> title_count:
997            raise IOError #FIXME make this nicer
998        for i, value in enumerate(line):
999            attribute_dic.setdefault(titles_stripped[i],[]).append(value)
1000       
1001    return attribute_dic, title_index_dic
1002
1003
1004#Auxiliary
1005def write_obj(filename,x,y,z):
1006    """Store x,y,z vectors into filename (obj format)
1007       Vectors are assumed to have dimension (M,3) where
1008       M corresponds to the number elements.
1009       triangles are assumed to be disconnected
1010
1011       The three numbers in each vector correspond to three vertices,
1012
1013       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
1014
1015    """
1016    #print 'Writing obj to %s' % filename
1017
1018    import os.path
1019
1020    root, ext = os.path.splitext(filename)
1021    if ext == '.obj':
1022        FN = filename
1023    else:
1024        FN = filename + '.obj'
1025
1026
1027    outfile = open(FN, 'wb')
1028    outfile.write("# Triangulation as an obj file\n")
1029
1030    M, N = x.shape
1031    assert N==3  #Assuming three vertices per element
1032
1033    for i in range(M):
1034        for j in range(N):
1035            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
1036
1037    for i in range(M):
1038        base = i*N
1039        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
1040
1041    outfile.close()
1042
1043
1044#########################################################
1045#Conversion routines
1046########################################################
1047
1048def sww2obj(basefilename, size):
1049    """Convert netcdf based data output to obj
1050    """
1051    from Scientific.IO.NetCDF import NetCDFFile
1052
1053    from Numeric import Float, zeros
1054
1055    #Get NetCDF
1056    FN = create_filename('.', basefilename, 'sww', size)
1057    print 'Reading from ', FN
1058    fid = NetCDFFile(FN, 'r')  #Open existing file for read
1059
1060
1061    # Get the variables
1062    x = fid.variables['x']
1063    y = fid.variables['y']
1064    z = fid.variables['elevation']
1065    time = fid.variables['time']
1066    stage = fid.variables['stage']
1067
1068    M = size  #Number of lines
1069    xx = zeros((M,3), Float)
1070    yy = zeros((M,3), Float)
1071    zz = zeros((M,3), Float)
1072
1073    for i in range(M):
1074        for j in range(3):
1075            xx[i,j] = x[i+j*M]
1076            yy[i,j] = y[i+j*M]
1077            zz[i,j] = z[i+j*M]
1078
1079    #Write obj for bathymetry
1080    FN = create_filename('.', basefilename, 'obj', size)
1081    write_obj(FN,xx,yy,zz)
1082
1083
1084    #Now read all the data with variable information, combine with
1085    #x,y info and store as obj
1086
1087    for k in range(len(time)):
1088        t = time[k]
1089        print 'Processing timestep %f' %t
1090
1091        for i in range(M):
1092            for j in range(3):
1093                zz[i,j] = stage[k,i+j*M]
1094
1095
1096        #Write obj for variable data
1097        #FN = create_filename(basefilename, 'obj', size, time=t)
1098        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
1099        write_obj(FN,xx,yy,zz)
1100
1101
1102def dat2obj(basefilename):
1103    """Convert line based data output to obj
1104    FIXME: Obsolete?
1105    """
1106
1107    import glob, os
1108    from anuga.config import data_dir
1109
1110
1111    #Get bathymetry and x,y's
1112    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
1113
1114    from Numeric import zeros, Float
1115
1116    M = len(lines)  #Number of lines
1117    x = zeros((M,3), Float)
1118    y = zeros((M,3), Float)
1119    z = zeros((M,3), Float)
1120
1121    ##i = 0
1122    for i, line in enumerate(lines):
1123        tokens = line.split()
1124        values = map(float,tokens)
1125
1126        for j in range(3):
1127            x[i,j] = values[j*3]
1128            y[i,j] = values[j*3+1]
1129            z[i,j] = values[j*3+2]
1130
1131        ##i += 1
1132
1133
1134    #Write obj for bathymetry
1135    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
1136
1137
1138    #Now read all the data files with variable information, combine with
1139    #x,y info
1140    #and store as obj
1141
1142    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
1143
1144    for filename in files:
1145        print 'Processing %s' % filename
1146
1147        lines = open(data_dir+os.sep+filename,'r').readlines()
1148        assert len(lines) == M
1149        root, ext = os.path.splitext(filename)
1150
1151        #Get time from filename
1152        i0 = filename.find('_time=')
1153        if i0 == -1:
1154            #Skip bathymetry file
1155            continue
1156
1157        i0 += 6  #Position where time starts
1158        i1 = filename.find('.dat')
1159
1160        if i1 > i0:
1161            t = float(filename[i0:i1])
1162        else:
1163            raise DataTimeError, 'Hmmmm'
1164
1165
1166
1167        ##i = 0
1168        for i, line in enumerate(lines):
1169            tokens = line.split()
1170            values = map(float,tokens)
1171
1172            for j in range(3):
1173                z[i,j] = values[j]
1174
1175            ##i += 1
1176
1177        #Write obj for variable data
1178        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
1179
1180
1181def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
1182    """Read netcdf filename1, pick timesteps first:step:last and save to
1183    nettcdf file filename2
1184    """
1185    from Scientific.IO.NetCDF import NetCDFFile
1186
1187    #Get NetCDF
1188    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
1189    outfile = NetCDFFile(filename2, 'w')  #Open new file
1190
1191
1192    #Copy dimensions
1193    for d in infile.dimensions:
1194        outfile.createDimension(d, infile.dimensions[d])
1195
1196    for name in infile.variables:
1197        var = infile.variables[name]
1198        outfile.createVariable(name, var.typecode(), var.dimensions)
1199
1200
1201    #Copy the static variables
1202    for name in infile.variables:
1203        if name == 'time' or name == 'stage':
1204            pass
1205        else:
1206            #Copy
1207            outfile.variables[name][:] = infile.variables[name][:]
1208
1209    #Copy selected timesteps
1210    time = infile.variables['time']
1211    stage = infile.variables['stage']
1212
1213    newtime = outfile.variables['time']
1214    newstage = outfile.variables['stage']
1215
1216    if last is None:
1217        last = len(time)
1218
1219    selection = range(first, last, step)
1220    for i, j in enumerate(selection):
1221        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
1222        newtime[i] = time[j]
1223        newstage[i,:] = stage[j,:]
1224
1225    #Close
1226    infile.close()
1227    outfile.close()
1228
1229
1230#Get data objects
1231def get_dataobject(domain, mode='w'):
1232    """Return instance of class of given format using filename
1233    """
1234
1235    cls = eval('Data_format_%s' %domain.format)
1236    return cls(domain, mode)
1237
1238
1239
1240
1241def dem2pts(basename_in, basename_out=None,
1242            easting_min=None, easting_max=None,
1243            northing_min=None, northing_max=None,
1244            use_cache=False, verbose=False,):
1245    """Read Digitial Elevation model from the following NetCDF format (.dem)
1246
1247    Example:
1248
1249    ncols         3121
1250    nrows         1800
1251    xllcorner     722000
1252    yllcorner     5893000
1253    cellsize      25
1254    NODATA_value  -9999
1255    138.3698 137.4194 136.5062 135.5558 ..........
1256
1257    Convert to NetCDF pts format which is
1258
1259    points:  (Nx2) Float array
1260    elevation: N Float array
1261    """
1262
1263
1264
1265    kwargs = {'basename_out': basename_out,
1266              'easting_min': easting_min,
1267              'easting_max': easting_max,
1268              'northing_min': northing_min,
1269              'northing_max': northing_max,
1270              'verbose': verbose}
1271
1272    if use_cache is True:
1273        from caching import cache
1274        result = cache(_dem2pts, basename_in, kwargs,
1275                       dependencies = [basename_in + '.dem'],
1276                       verbose = verbose)
1277
1278    else:
1279        result = apply(_dem2pts, [basename_in], kwargs)
1280
1281    return result
1282
1283
1284def _dem2pts(basename_in, basename_out=None, verbose=False,
1285            easting_min=None, easting_max=None,
1286            northing_min=None, northing_max=None):
1287    """Read Digitial Elevation model from the following NetCDF format (.dem)
1288
1289    Internal function. See public function dem2pts for details.
1290    """
1291
1292    # FIXME: Can this be written feasibly using write_pts?
1293
1294    import os
1295    from Scientific.IO.NetCDF import NetCDFFile
1296    from Numeric import Float, zeros, reshape, sum
1297
1298    root = basename_in
1299
1300    # Get NetCDF
1301    infile = NetCDFFile(root + '.dem', 'r')  # Open existing netcdf file for read
1302
1303    if verbose: print 'Reading DEM from %s' %(root + '.dem')
1304
1305    ncols = infile.ncols[0]
1306    nrows = infile.nrows[0]
1307    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
1308    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
1309    cellsize = infile.cellsize[0]
1310    NODATA_value = infile.NODATA_value[0]
1311    dem_elevation = infile.variables['elevation']
1312
1313    zone = infile.zone[0]
1314    false_easting = infile.false_easting[0]
1315    false_northing = infile.false_northing[0]
1316
1317    # Text strings
1318    projection = infile.projection
1319    datum = infile.datum
1320    units = infile.units
1321
1322
1323    # Get output file
1324    if basename_out == None:
1325        ptsname = root + '.pts'
1326    else:
1327        ptsname = basename_out + '.pts'
1328
1329    if verbose: print 'Store to NetCDF file %s' %ptsname
1330    # NetCDF file definition
1331    outfile = NetCDFFile(ptsname, 'w')
1332
1333    # Create new file
1334    outfile.institution = 'Geoscience Australia'
1335    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
1336                      'of spatial point data'
1337    # Assign default values
1338    if easting_min is None: easting_min = xllcorner
1339    if easting_max is None: easting_max = xllcorner + ncols*cellsize
1340    if northing_min is None: northing_min = yllcorner
1341    if northing_max is None: northing_max = yllcorner + nrows*cellsize
1342
1343    # Compute offsets to update georeferencing
1344    easting_offset = xllcorner - easting_min
1345    northing_offset = yllcorner - northing_min
1346
1347    # Georeferencing
1348    outfile.zone = zone
1349    outfile.xllcorner = easting_min # Easting of lower left corner
1350    outfile.yllcorner = northing_min # Northing of lower left corner
1351    outfile.false_easting = false_easting
1352    outfile.false_northing = false_northing
1353
1354    outfile.projection = projection
1355    outfile.datum = datum
1356    outfile.units = units
1357
1358
1359    # Grid info (FIXME: probably not going to be used, but heck)
1360    outfile.ncols = ncols
1361    outfile.nrows = nrows
1362
1363    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
1364    totalnopoints = nrows*ncols
1365
1366    # Calculating number of NODATA_values for each row in clipped region
1367    # FIXME: use array operations to do faster
1368    nn = 0
1369    k = 0
1370    i1_0 = 0
1371    j1_0 = 0
1372    thisj = 0
1373    thisi = 0
1374    for i in range(nrows):
1375        y = (nrows-i-1)*cellsize + yllcorner
1376        for j in range(ncols):
1377            x = j*cellsize + xllcorner
1378            if easting_min <= x <= easting_max and \
1379               northing_min <= y <= northing_max:
1380                thisj = j
1381                thisi = i
1382                if dem_elevation_r[i,j] == NODATA_value: nn += 1
1383
1384                if k == 0:
1385                    i1_0 = i
1386                    j1_0 = j
1387                k += 1
1388
1389    index1 = j1_0
1390    index2 = thisj
1391
1392    # Dimension definitions
1393    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
1394    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
1395
1396    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
1397    nopoints = clippednopoints-nn
1398
1399    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
1400
1401    if verbose:
1402        print 'There are %d values in the elevation' %totalnopoints
1403        print 'There are %d values in the clipped elevation' %clippednopoints
1404        print 'There are %d NODATA_values in the clipped elevation' %nn
1405
1406    outfile.createDimension('number_of_points', nopoints)
1407    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1408
1409    # Variable definitions
1410    outfile.createVariable('points', Float, ('number_of_points',
1411                                             'number_of_dimensions'))
1412    outfile.createVariable('elevation', Float, ('number_of_points',))
1413
1414    # Get handles to the variables
1415    points = outfile.variables['points']
1416    elevation = outfile.variables['elevation']
1417
1418    lenv = index2-index1+1
1419    # Store data
1420    global_index = 0
1421    # for i in range(nrows):
1422    for i in range(i1_0,thisi+1,1):
1423        if verbose and i%((nrows+10)/10)==0:
1424            print 'Processing row %d of %d' %(i, nrows)
1425
1426        lower_index = global_index
1427
1428        v = dem_elevation_r[i,index1:index2+1]
1429        no_NODATA = sum(v == NODATA_value)
1430        if no_NODATA > 0:
1431            newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA
1432        else:
1433            newcols = lenv # ncols_in_bounding_box
1434
1435        telev = zeros(newcols, Float)
1436        tpoints = zeros((newcols, 2), Float)
1437
1438        local_index = 0
1439
1440        y = (nrows-i-1)*cellsize + yllcorner
1441        #for j in range(ncols):
1442        for j in range(j1_0,index2+1,1):
1443
1444            x = j*cellsize + xllcorner
1445            if easting_min <= x <= easting_max and \
1446               northing_min <= y <= northing_max and \
1447               dem_elevation_r[i,j] <> NODATA_value:
1448                tpoints[local_index, :] = [x-easting_min,y-northing_min]
1449                telev[local_index] = dem_elevation_r[i, j]
1450                global_index += 1
1451                local_index += 1
1452
1453        upper_index = global_index
1454
1455        if upper_index == lower_index + newcols:
1456            points[lower_index:upper_index, :] = tpoints
1457            elevation[lower_index:upper_index] = telev
1458
1459    assert global_index == nopoints, 'index not equal to number of points'
1460
1461    infile.close()
1462    outfile.close()
1463
1464
1465
1466def _read_hecras_cross_sections(lines):
1467    """Return block of surface lines for each cross section
1468    Starts with SURFACE LINE,
1469    Ends with END CROSS-SECTION
1470    """
1471
1472    points = []
1473
1474    reading_surface = False
1475    for i, line in enumerate(lines):
1476
1477        if len(line.strip()) == 0:    #Ignore blanks
1478            continue
1479
1480        if lines[i].strip().startswith('SURFACE LINE'):
1481            reading_surface = True
1482            continue
1483
1484        if lines[i].strip().startswith('END') and reading_surface:
1485            yield points
1486            reading_surface = False
1487            points = []
1488
1489        if reading_surface:
1490            fields = line.strip().split(',')
1491            easting = float(fields[0])
1492            northing = float(fields[1])
1493            elevation = float(fields[2])
1494            points.append([easting, northing, elevation])
1495
1496
1497
1498
1499def hecras_cross_sections2pts(basename_in,
1500                              basename_out=None,
1501                              verbose=False):
1502    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
1503
1504    Example:
1505
1506
1507# RAS export file created on Mon 15Aug2005 11:42
1508# by HEC-RAS Version 3.1.1
1509
1510BEGIN HEADER:
1511  UNITS: METRIC
1512  DTM TYPE: TIN
1513  DTM: v:\1\cit\perth_topo\river_tin
1514  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
1515  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
1516  MAP PROJECTION: UTM
1517  PROJECTION ZONE: 50
1518  DATUM: AGD66
1519  VERTICAL DATUM:
1520  NUMBER OF REACHES:  19
1521  NUMBER OF CROSS-SECTIONS:  14206
1522END HEADER:
1523
1524
1525Only the SURFACE LINE data of the following form will be utilised
1526
1527  CROSS-SECTION:
1528    STREAM ID:Southern-Wungong
1529    REACH ID:Southern-Wungong
1530    STATION:19040.*
1531    CUT LINE:
1532      405548.671603161 , 6438142.7594925
1533      405734.536092045 , 6438326.10404912
1534      405745.130459356 , 6438331.48627354
1535      405813.89633823 , 6438368.6272789
1536    SURFACE LINE:
1537     405548.67,   6438142.76,   35.37
1538     405552.24,   6438146.28,   35.41
1539     405554.78,   6438148.78,   35.44
1540     405555.80,   6438149.79,   35.44
1541     405559.37,   6438153.31,   35.45
1542     405560.88,   6438154.81,   35.44
1543     405562.93,   6438156.83,   35.42
1544     405566.50,   6438160.35,   35.38
1545     405566.99,   6438160.83,   35.37
1546     ...
1547   END CROSS-SECTION
1548
1549    Convert to NetCDF pts format which is
1550
1551    points:  (Nx2) Float array
1552    elevation: N Float array
1553    """
1554
1555    import os
1556    from Scientific.IO.NetCDF import NetCDFFile
1557    from Numeric import Float, zeros, reshape
1558
1559    root = basename_in
1560
1561    #Get ASCII file
1562    infile = open(root + '.sdf', 'r')  #Open SDF file for read
1563
1564    if verbose: print 'Reading DEM from %s' %(root + '.sdf')
1565
1566    lines = infile.readlines()
1567    infile.close()
1568
1569    if verbose: print 'Converting to pts format'
1570
1571    i = 0
1572    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
1573        i += 1
1574
1575    assert lines[i].strip().upper() == 'BEGIN HEADER:'
1576    i += 1
1577
1578    assert lines[i].strip().upper().startswith('UNITS:')
1579    units = lines[i].strip().split()[1]
1580    i += 1
1581
1582    assert lines[i].strip().upper().startswith('DTM TYPE:')
1583    i += 1
1584
1585    assert lines[i].strip().upper().startswith('DTM:')
1586    i += 1
1587
1588    assert lines[i].strip().upper().startswith('STREAM')
1589    i += 1
1590
1591    assert lines[i].strip().upper().startswith('CROSS')
1592    i += 1
1593
1594    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
1595    projection = lines[i].strip().split(':')[1]
1596    i += 1
1597
1598    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
1599    zone = int(lines[i].strip().split(':')[1])
1600    i += 1
1601
1602    assert lines[i].strip().upper().startswith('DATUM:')
1603    datum = lines[i].strip().split(':')[1]
1604    i += 1
1605
1606    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
1607    i += 1
1608
1609    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
1610    i += 1
1611
1612    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
1613    number_of_cross_sections = int(lines[i].strip().split(':')[1])
1614    i += 1
1615
1616
1617    #Now read all points
1618    points = []
1619    elevation = []
1620    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
1621        for k, entry in enumerate(entries):
1622            points.append(entry[:2])
1623            elevation.append(entry[2])
1624
1625
1626    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
1627          %(j+1, number_of_cross_sections)
1628    assert j+1 == number_of_cross_sections, msg
1629
1630    #Get output file
1631    if basename_out == None:
1632        ptsname = root + '.pts'
1633    else:
1634        ptsname = basename_out + '.pts'
1635
1636    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
1637    geo = Geospatial_data(points, {"elevation":elevation},
1638                          verbose=verbose, geo_reference=geo_ref)
1639    geo.export_points_file(ptsname)
1640
1641def export_grid(basename_in, extra_name_out = None,
1642                quantities = None, # defaults to elevation
1643                timestep = None,
1644                reduction = None,
1645                cellsize = 10,
1646                NODATA_value = -9999,
1647                easting_min = None,
1648                easting_max = None,
1649                northing_min = None,
1650                northing_max = None,
1651                verbose = False,
1652                origin = None,
1653                datum = 'WGS84',
1654                format = 'ers'):
1655    """
1656   
1657    Wrapper for sww2dem. - see sww2dem to find out what most of the
1658    parameters do.
1659
1660    Quantities is a list of quantities.  Each quantity will be
1661    calculated for each sww file.
1662
1663    This returns the basenames of the files returned, which is made up
1664    of the dir and all of the file name, except the extension.
1665
1666    This function returns the names of the files produced.
1667
1668    It will also produce as many output files as there are input sww files.
1669    """
1670   
1671    if quantities is None:
1672        quantities = ['elevation']
1673       
1674    if type(quantities) is str:
1675            quantities = [quantities]
1676
1677    # How many sww files are there?
1678    dir, base = os.path.split(basename_in)
1679    #print "basename_in",basename_in
1680    #print "base",base
1681
1682    iterate_over = get_all_swwfiles(dir,base,verbose)
1683   
1684    if dir == "":
1685        dir = "." # Unix compatibility
1686   
1687    files_out = []
1688    #print 'sww_file',iterate_over
1689    for sww_file in iterate_over:
1690        for quantity in quantities:
1691            if extra_name_out is None:
1692                basename_out = sww_file + '_' + quantity
1693            else:
1694                basename_out = sww_file + '_' + quantity + '_' \
1695                               + extra_name_out
1696#            print "basename_out", basename_out
1697       
1698            file_out = sww2dem(dir+sep+sww_file, dir+sep+basename_out,
1699                               quantity, 
1700                               timestep,
1701                               reduction,
1702                               cellsize,
1703                               NODATA_value,
1704                               easting_min,
1705                               easting_max,
1706                               northing_min,
1707                               northing_max,
1708                               verbose,
1709                               origin,
1710                               datum,
1711                               format)
1712            files_out.append(file_out)
1713    #print "basenames_out after",basenames_out
1714    return files_out
1715
1716
1717def get_timeseries(production_dirs, output_dir, scenario_name, gauges_dir_name,
1718                   plot_quantity, generate_fig = False,
1719                   reportname = None, surface = False, time_min = None,
1720                   time_max = None, title_on = False, verbose = True,
1721                   nodes=None):
1722    """
1723    nodes - number of processes used.
1724
1725    warning - this function has no tests
1726    """
1727    if reportname == None:
1728        report = False
1729    else:
1730        report = True
1731       
1732    if nodes is None:
1733        is_parallel = False
1734    else:
1735        is_parallel = True
1736       
1737    # Generate figures
1738    swwfiles = {}
1739   
1740    if is_parallel is True:   
1741        for i in range(nodes):
1742            print 'Sending node %d of %d' %(i,nodes)
1743            swwfiles = {}
1744            if not reportname == None:
1745                reportname = report_name + '_%s' %(i)
1746            for label_id in production_dirs.keys():
1747                if label_id == 'boundaries':
1748                    swwfile = best_boundary_sww
1749                else:
1750                    file_loc = output_dir + label_id + sep
1751                    sww_extra = '_P%s_%s' %(i,nodes)
1752                    swwfile = file_loc + scenario_name + sww_extra + '.sww'
1753                    print 'swwfile',swwfile
1754                    swwfiles[swwfile] = label_id
1755
1756                texname, elev_output = sww2timeseries(swwfiles,
1757                                              gauges_dir_name,
1758                                              production_dirs,
1759                                              report = report,
1760                                              reportname = reportname,
1761                                              plot_quantity = plot_quantity,
1762                                              generate_fig = generate_fig,
1763                                              surface = surface,
1764                                              time_min = time_min,
1765                                              time_max = time_max,
1766                                              title_on = title_on,
1767                                              verbose = verbose)
1768    else:   
1769        for label_id in production_dirs.keys():       
1770            if label_id == 'boundaries':
1771                print 'boundaries'
1772                file_loc = project.boundaries_in_dir
1773                swwfile = project.boundaries_dir_name3 + '.sww'
1774                #  swwfile = boundary_dir_filename
1775            else:
1776                file_loc = output_dir + label_id + sep
1777                swwfile = file_loc + scenario_name + '.sww'
1778            swwfiles[swwfile] = label_id
1779       
1780        texname, elev_output = sww2timeseries(swwfiles,
1781                                              gauges_dir_name,
1782                                              production_dirs,
1783                                              report = report,
1784                                              reportname = reportname,
1785                                              plot_quantity = plot_quantity,
1786                                              generate_fig = generate_fig,
1787                                              surface = surface,
1788                                              time_min = time_min,
1789                                              time_max = time_max,
1790                                              title_on = title_on,
1791                                              verbose = verbose)
1792                                         
1793
1794   
1795def sww2dem(basename_in, basename_out = None,
1796            quantity = None, # defaults to elevation
1797            timestep = None,
1798            reduction = None,
1799            cellsize = 10,
1800            NODATA_value = -9999,
1801            easting_min = None,
1802            easting_max = None,
1803            northing_min = None,
1804            northing_max = None,
1805            verbose = False,
1806            origin = None,
1807            datum = 'WGS84',
1808            format = 'ers'):
1809
1810    """Read SWW file and convert to Digitial Elevation model format
1811    (.asc or .ers)
1812
1813    Example (ASC):
1814
1815    ncols         3121
1816    nrows         1800
1817    xllcorner     722000
1818    yllcorner     5893000
1819    cellsize      25
1820    NODATA_value  -9999
1821    138.3698 137.4194 136.5062 135.5558 ..........
1822
1823    Also write accompanying file with same basename_in but extension .prj
1824    used to fix the UTM zone, datum, false northings and eastings.
1825
1826    The prj format is assumed to be as
1827
1828    Projection    UTM
1829    Zone          56
1830    Datum         WGS84
1831    Zunits        NO
1832    Units         METERS
1833    Spheroid      WGS84
1834    Xshift        0.0000000000
1835    Yshift        10000000.0000000000
1836    Parameters
1837
1838
1839    The parameter quantity must be the name of an existing quantity or
1840    an expression involving existing quantities. The default is
1841    'elevation'. Quantity is not a list of quantities.
1842
1843    if timestep (an index) is given, output quantity at that timestep
1844
1845    if reduction is given use that to reduce quantity over all timesteps.
1846
1847    datum
1848
1849    format can be either 'asc' or 'ers'
1850    """
1851
1852    import sys
1853    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, \
1854         sometrue
1855    from Numeric import array2string
1856
1857    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
1858         separate_points_by_polygon
1859    from anuga.abstract_2d_finite_volumes.util import \
1860         apply_expression_to_dictionary
1861
1862    msg = 'Format must be either asc or ers'
1863    assert format.lower() in ['asc', 'ers'], msg
1864
1865
1866    false_easting = 500000
1867    false_northing = 10000000
1868
1869    if quantity is None:
1870        quantity = 'elevation'
1871       
1872    if reduction is None:
1873        reduction = max
1874
1875    if basename_out is None:
1876        basename_out = basename_in + '_%s' %quantity
1877
1878    if quantity_formula.has_key(quantity):
1879        quantity = quantity_formula[quantity]
1880       
1881    swwfile = basename_in + '.sww'
1882    demfile = basename_out + '.' + format
1883    # Note the use of a .ers extension is optional (write_ermapper_grid will
1884    # deal with either option
1885   
1886    #if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors
1887   
1888    # Read sww file
1889    if verbose: 
1890        print 'Reading from %s' %swwfile
1891        print 'Output directory is %s' %basename_out
1892   
1893    from Scientific.IO.NetCDF import NetCDFFile
1894    fid = NetCDFFile(swwfile)
1895
1896    #Get extent and reference
1897    x = fid.variables['x'][:]
1898    y = fid.variables['y'][:]
1899    volumes = fid.variables['volumes'][:]
1900    if timestep is not None:
1901        times = fid.variables['time'][timestep]
1902    else:
1903        times = fid.variables['time'][:]
1904
1905    number_of_timesteps = fid.dimensions['number_of_timesteps']
1906    number_of_points = fid.dimensions['number_of_points']
1907   
1908    if origin is None:
1909
1910        # Get geo_reference
1911        # sww files don't have to have a geo_ref
1912        try:
1913            geo_reference = Geo_reference(NetCDFObject=fid)
1914        except AttributeError, e:
1915            geo_reference = Geo_reference() # Default georef object
1916
1917        xllcorner = geo_reference.get_xllcorner()
1918        yllcorner = geo_reference.get_yllcorner()
1919        zone = geo_reference.get_zone()
1920    else:
1921        zone = origin[0]
1922        xllcorner = origin[1]
1923        yllcorner = origin[2]
1924
1925
1926
1927    # FIXME: Refactor using code from Interpolation_function.statistics
1928    # (in interpolate.py)
1929    # Something like print swwstats(swwname)
1930    if verbose:
1931        print '------------------------------------------------'
1932        print 'Statistics of SWW file:'
1933        print '  Name: %s' %swwfile
1934        print '  Reference:'
1935        print '    Lower left corner: [%f, %f]'\
1936              %(xllcorner, yllcorner)
1937        if timestep is not None:
1938            print '    Time: %f' %(times)
1939        else:
1940            print '    Start time: %f' %fid.starttime[0]
1941        print '  Extent:'
1942        print '    x [m] in [%f, %f], len(x) == %d'\
1943              %(min(x.flat), max(x.flat), len(x.flat))
1944        print '    y [m] in [%f, %f], len(y) == %d'\
1945              %(min(y.flat), max(y.flat), len(y.flat))
1946        if timestep is not None:
1947            print '    t [s] = %f, len(t) == %d' %(times, 1)
1948        else:
1949            print '    t [s] in [%f, %f], len(t) == %d'\
1950              %(min(times), max(times), len(times))
1951        print '  Quantities [SI units]:'
1952        for name in ['stage', 'xmomentum', 'ymomentum']:
1953            q = fid.variables[name][:].flat
1954            if timestep is not None:
1955                q = q[timestep*len(x):(timestep+1)*len(x)]
1956            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1957        for name in ['elevation']:
1958            q = fid.variables[name][:].flat
1959            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1960           
1961    # Get quantity and reduce if applicable
1962    if verbose: print 'Processing quantity %s' %quantity
1963
1964    # Turn NetCDF objects into Numeric arrays
1965    quantity_dict = {}
1966    for name in fid.variables.keys():
1967        quantity_dict[name] = fid.variables[name][:] 
1968
1969
1970    # Convert quantity expression to quantities found in sww file   
1971    q = apply_expression_to_dictionary(quantity, quantity_dict)
1972
1973    if len(q.shape) == 2:
1974        #q has a time component and needs to be reduced along
1975        #the temporal dimension
1976        if verbose: print 'Reducing quantity %s' %quantity
1977        q_reduced = zeros( number_of_points, Float )
1978       
1979        if timestep is not None:
1980            for k in range(number_of_points):
1981                q_reduced[k] = q[timestep,k] 
1982        else:
1983            for k in range(number_of_points):
1984                q_reduced[k] = reduction( q[:,k] )
1985
1986        q = q_reduced
1987
1988    #Post condition: Now q has dimension: number_of_points
1989    assert len(q.shape) == 1
1990    assert q.shape[0] == number_of_points
1991
1992
1993    if verbose:
1994        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
1995
1996
1997    #Create grid and update xll/yll corner and x,y
1998
1999    #Relative extent
2000    if easting_min is None:
2001        xmin = min(x)
2002    else:
2003        xmin = easting_min - xllcorner
2004
2005    if easting_max is None:
2006        xmax = max(x)
2007    else:
2008        xmax = easting_max - xllcorner
2009
2010    if northing_min is None:
2011        ymin = min(y)
2012    else:
2013        ymin = northing_min - yllcorner
2014
2015    if northing_max is None:
2016        ymax = max(y)
2017    else:
2018        ymax = northing_max - yllcorner
2019
2020
2021
2022    if verbose: print 'Creating grid'
2023    ncols = int((xmax-xmin)/cellsize)+1
2024    nrows = int((ymax-ymin)/cellsize)+1
2025
2026
2027    #New absolute reference and coordinates
2028    newxllcorner = xmin+xllcorner
2029    newyllcorner = ymin+yllcorner
2030
2031    x = x+xllcorner-newxllcorner
2032    y = y+yllcorner-newyllcorner
2033   
2034    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2035    assert len(vertex_points.shape) == 2
2036
2037    grid_points = zeros ( (ncols*nrows, 2), Float )
2038
2039
2040    for i in xrange(nrows):
2041        if format.lower() == 'asc':
2042            yg = i*cellsize
2043        else:
2044        #this will flip the order of the y values for ers
2045            yg = (nrows-i)*cellsize
2046
2047        for j in xrange(ncols):
2048            xg = j*cellsize
2049            k = i*ncols + j
2050
2051            grid_points[k,0] = xg
2052            grid_points[k,1] = yg
2053
2054    #Interpolate
2055    from anuga.fit_interpolate.interpolate import Interpolate
2056
2057    # Remove loners from vertex_points, volumes here
2058    vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
2059    #export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
2060    #import sys; sys.exit()
2061    interp = Interpolate(vertex_points, volumes, verbose = verbose)
2062
2063    #Interpolate using quantity values
2064    if verbose: print 'Interpolating'
2065    grid_values = interp.interpolate(q, grid_points).flat
2066
2067
2068    if verbose:
2069        print 'Interpolated values are in [%f, %f]' %(min(grid_values),
2070                                                      max(grid_values))
2071
2072    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2073    P = interp.mesh.get_boundary_polygon()
2074    outside_indices = outside_polygon(grid_points, P, closed=True)
2075
2076    for i in outside_indices:
2077        grid_values[i] = NODATA_value
2078
2079
2080
2081
2082    if format.lower() == 'ers':
2083        # setup ERS header information
2084        grid_values = reshape(grid_values,(nrows, ncols))
2085        header = {}
2086        header['datum'] = '"' + datum + '"'
2087        # FIXME The use of hardwired UTM and zone number needs to be made optional
2088        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
2089        header['projection'] = '"UTM-' + str(zone) + '"'
2090        header['coordinatetype'] = 'EN'
2091        if header['coordinatetype'] == 'LL':
2092            header['longitude'] = str(newxllcorner)
2093            header['latitude'] = str(newyllcorner)
2094        elif header['coordinatetype'] == 'EN':
2095            header['eastings'] = str(newxllcorner)
2096            header['northings'] = str(newyllcorner)
2097        header['nullcellvalue'] = str(NODATA_value)
2098        header['xdimension'] = str(cellsize)
2099        header['ydimension'] = str(cellsize)
2100        header['value'] = '"' + quantity + '"'
2101        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
2102
2103
2104        #Write
2105        if verbose: print 'Writing %s' %demfile
2106        import ermapper_grids
2107        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
2108
2109        fid.close()
2110    else:
2111        #Write to Ascii format
2112
2113        #Write prj file
2114        prjfile = basename_out + '.prj'
2115
2116        if verbose: print 'Writing %s' %prjfile
2117        prjid = open(prjfile, 'w')
2118        prjid.write('Projection    %s\n' %'UTM')
2119        prjid.write('Zone          %d\n' %zone)
2120        prjid.write('Datum         %s\n' %datum)
2121        prjid.write('Zunits        NO\n')
2122        prjid.write('Units         METERS\n')
2123        prjid.write('Spheroid      %s\n' %datum)
2124        prjid.write('Xshift        %d\n' %false_easting)
2125        prjid.write('Yshift        %d\n' %false_northing)
2126        prjid.write('Parameters\n')
2127        prjid.close()
2128
2129
2130
2131        if verbose: print 'Writing %s' %demfile
2132
2133        ascid = open(demfile, 'w')
2134
2135        ascid.write('ncols         %d\n' %ncols)
2136        ascid.write('nrows         %d\n' %nrows)
2137        ascid.write('xllcorner     %d\n' %newxllcorner)
2138        ascid.write('yllcorner     %d\n' %newyllcorner)
2139        ascid.write('cellsize      %f\n' %cellsize)
2140        ascid.write('NODATA_value  %d\n' %NODATA_value)
2141
2142
2143        #Get bounding polygon from mesh
2144        #P = interp.mesh.get_boundary_polygon()
2145        #inside_indices = inside_polygon(grid_points, P)
2146
2147        for i in range(nrows):
2148            if verbose and i%((nrows+10)/10)==0:
2149                print 'Doing row %d of %d' %(i, nrows)
2150
2151            base_index = (nrows-i-1)*ncols
2152
2153            slice = grid_values[base_index:base_index+ncols]
2154            s = array2string(slice, max_line_width=sys.maxint)
2155            ascid.write(s[1:-1] + '\n')
2156
2157
2158            #print
2159            #for j in range(ncols):
2160            #    index = base_index+j#
2161            #    print grid_values[index],
2162            #    ascid.write('%f ' %grid_values[index])
2163            #ascid.write('\n')
2164
2165        #Close
2166        ascid.close()
2167        fid.close()
2168        return basename_out
2169
2170#Backwards compatibility
2171def sww2asc(basename_in, basename_out = None,
2172            quantity = None,
2173            timestep = None,
2174            reduction = None,
2175            cellsize = 10,
2176            verbose = False,
2177            origin = None):
2178    print 'sww2asc will soon be obsoleted - please use sww2dem'
2179    sww2dem(basename_in,
2180            basename_out = basename_out,
2181            quantity = quantity,
2182            timestep = timestep,
2183            reduction = reduction,
2184            cellsize = cellsize,
2185            verbose = verbose,
2186            origin = origin,
2187        datum = 'WGS84',
2188        format = 'asc')
2189
2190def sww2ers(basename_in, basename_out = None,
2191            quantity = None,
2192            timestep = None,
2193            reduction = None,
2194            cellsize = 10,
2195            verbose = False,
2196            origin = None,
2197            datum = 'WGS84'):
2198    print 'sww2ers will soon be obsoleted - please use sww2dem'
2199    sww2dem(basename_in,
2200            basename_out = basename_out,
2201            quantity = quantity,
2202            timestep = timestep,
2203            reduction = reduction,
2204            cellsize = cellsize,
2205            verbose = verbose,
2206            origin = origin,
2207            datum = datum,
2208            format = 'ers')
2209################################# END COMPATIBILITY ##############
2210
2211
2212
2213def sww2pts(basename_in, basename_out=None,
2214            data_points=None,
2215            quantity=None,
2216            timestep=None,
2217            reduction=None,
2218            NODATA_value=-9999,
2219            verbose=False,
2220            origin=None):
2221            #datum = 'WGS84')
2222
2223
2224    """Read SWW file and convert to interpolated values at selected points
2225
2226    The parameter quantity' must be the name of an existing quantity or
2227    an expression involving existing quantities. The default is
2228    'elevation'.
2229
2230    if timestep (an index) is given, output quantity at that timestep
2231
2232    if reduction is given use that to reduce quantity over all timesteps.
2233
2234    data_points (Nx2 array) give locations of points where quantity is to be computed.
2235   
2236    """
2237
2238    import sys
2239    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
2240    from Numeric import array2string
2241
2242    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
2243    from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
2244
2245    from anuga.geospatial_data.geospatial_data import Geospatial_data
2246
2247    if quantity is None:
2248        quantity = 'elevation'
2249
2250    if reduction is None:
2251        reduction = max
2252
2253    if basename_out is None:
2254        basename_out = basename_in + '_%s' %quantity
2255
2256    swwfile = basename_in + '.sww'
2257    ptsfile = basename_out + '.pts'
2258
2259    # Read sww file
2260    if verbose: print 'Reading from %s' %swwfile
2261    from Scientific.IO.NetCDF import NetCDFFile
2262    fid = NetCDFFile(swwfile)
2263
2264    # Get extent and reference
2265    x = fid.variables['x'][:]
2266    y = fid.variables['y'][:]
2267    volumes = fid.variables['volumes'][:]
2268
2269    number_of_timesteps = fid.dimensions['number_of_timesteps']
2270    number_of_points = fid.dimensions['number_of_points']
2271    if origin is None:
2272
2273        # Get geo_reference
2274        # sww files don't have to have a geo_ref
2275        try:
2276            geo_reference = Geo_reference(NetCDFObject=fid)
2277        except AttributeError, e:
2278            geo_reference = Geo_reference() #Default georef object
2279
2280        xllcorner = geo_reference.get_xllcorner()
2281        yllcorner = geo_reference.get_yllcorner()
2282        zone = geo_reference.get_zone()
2283    else:
2284        zone = origin[0]
2285        xllcorner = origin[1]
2286        yllcorner = origin[2]
2287
2288
2289
2290    # FIXME: Refactor using code from file_function.statistics
2291    # Something like print swwstats(swwname)
2292    if verbose:
2293        x = fid.variables['x'][:]
2294        y = fid.variables['y'][:]
2295        times = fid.variables['time'][:]
2296        print '------------------------------------------------'
2297        print 'Statistics of SWW file:'
2298        print '  Name: %s' %swwfile
2299        print '  Reference:'
2300        print '    Lower left corner: [%f, %f]'\
2301              %(xllcorner, yllcorner)
2302        print '    Start time: %f' %fid.starttime[0]
2303        print '  Extent:'
2304        print '    x [m] in [%f, %f], len(x) == %d'\
2305              %(min(x.flat), max(x.flat), len(x.flat))
2306        print '    y [m] in [%f, %f], len(y) == %d'\
2307              %(min(y.flat), max(y.flat), len(y.flat))
2308        print '    t [s] in [%f, %f], len(t) == %d'\
2309              %(min(times), max(times), len(times))
2310        print '  Quantities [SI units]:'
2311        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2312            q = fid.variables[name][:].flat
2313            print '    %s in [%f, %f]' %(name, min(q), max(q))
2314
2315
2316
2317    # Get quantity and reduce if applicable
2318    if verbose: print 'Processing quantity %s' %quantity
2319
2320    # Turn NetCDF objects into Numeric arrays
2321    quantity_dict = {}
2322    for name in fid.variables.keys():
2323        quantity_dict[name] = fid.variables[name][:]
2324
2325
2326
2327    # Convert quantity expression to quantities found in sww file   
2328    q = apply_expression_to_dictionary(quantity, quantity_dict)
2329
2330
2331
2332    if len(q.shape) == 2:
2333        # q has a time component and needs to be reduced along
2334        # the temporal dimension
2335        if verbose: print 'Reducing quantity %s' %quantity
2336        q_reduced = zeros( number_of_points, Float )
2337
2338        for k in range(number_of_points):
2339            q_reduced[k] = reduction( q[:,k] )
2340
2341        q = q_reduced
2342
2343    # Post condition: Now q has dimension: number_of_points
2344    assert len(q.shape) == 1
2345    assert q.shape[0] == number_of_points
2346
2347
2348    if verbose:
2349        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
2350
2351
2352    # Create grid and update xll/yll corner and x,y
2353    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2354    assert len(vertex_points.shape) == 2
2355
2356    # Interpolate
2357    from anuga.fit_interpolate.interpolate import Interpolate
2358    interp = Interpolate(vertex_points, volumes, verbose = verbose)
2359
2360    # Interpolate using quantity values
2361    if verbose: print 'Interpolating'
2362    interpolated_values = interp.interpolate(q, data_points).flat
2363
2364
2365    if verbose:
2366        print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
2367                                                      max(interpolated_values))
2368
2369    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2370    P = interp.mesh.get_boundary_polygon()
2371    outside_indices = outside_polygon(data_points, P, closed=True)
2372
2373    for i in outside_indices:
2374        interpolated_values[i] = NODATA_value
2375
2376
2377    # Store results   
2378    G = Geospatial_data(data_points=data_points,
2379                        attributes=interpolated_values)
2380
2381    G.export_points_file(ptsfile, absolute = True)
2382
2383    fid.close()
2384
2385
2386def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2387                                  use_cache = False,
2388                                  verbose = False):
2389    """Read Digitial Elevation model from the following ASCII format (.asc)
2390
2391    Example:
2392
2393    ncols         3121
2394    nrows         1800
2395    xllcorner     722000
2396    yllcorner     5893000
2397    cellsize      25
2398    NODATA_value  -9999
2399    138.3698 137.4194 136.5062 135.5558 ..........
2400
2401    Convert basename_in + '.asc' to NetCDF format (.dem)
2402    mimicking the ASCII format closely.
2403
2404
2405    An accompanying file with same basename_in but extension .prj must exist
2406    and is used to fix the UTM zone, datum, false northings and eastings.
2407
2408    The prj format is assumed to be as
2409
2410    Projection    UTM
2411    Zone          56
2412    Datum         WGS84
2413    Zunits        NO
2414    Units         METERS
2415    Spheroid      WGS84
2416    Xshift        0.0000000000
2417    Yshift        10000000.0000000000
2418    Parameters
2419    """
2420
2421
2422
2423    kwargs = {'basename_out': basename_out, 'verbose': verbose}
2424
2425    if use_cache is True:
2426        from caching import cache
2427        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
2428                       dependencies = [basename_in + '.asc',
2429                                       basename_in + '.prj'],
2430                       verbose = verbose)
2431
2432    else:
2433        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
2434
2435    return result
2436
2437
2438
2439
2440
2441
2442def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2443                                  verbose = False):
2444    """Read Digitial Elevation model from the following ASCII format (.asc)
2445
2446    Internal function. See public function convert_dem_from_ascii2netcdf for details.
2447    """
2448
2449    import os
2450    from Scientific.IO.NetCDF import NetCDFFile
2451    from Numeric import Float, array
2452
2453    #root, ext = os.path.splitext(basename_in)
2454    root = basename_in
2455
2456    ###########################################
2457    # Read Meta data
2458    if verbose: print 'Reading METADATA from %s' %root + '.prj'
2459    metadatafile = open(root + '.prj')
2460    metalines = metadatafile.readlines()
2461    metadatafile.close()
2462
2463    L = metalines[0].strip().split()
2464    assert L[0].strip().lower() == 'projection'
2465    projection = L[1].strip()                   #TEXT
2466
2467    L = metalines[1].strip().split()
2468    assert L[0].strip().lower() == 'zone'
2469    zone = int(L[1].strip())
2470
2471    L = metalines[2].strip().split()
2472    assert L[0].strip().lower() == 'datum'
2473    datum = L[1].strip()                        #TEXT
2474
2475    L = metalines[3].strip().split()
2476    assert L[0].strip().lower() == 'zunits'     #IGNORE
2477    zunits = L[1].strip()                       #TEXT
2478
2479    L = metalines[4].strip().split()
2480    assert L[0].strip().lower() == 'units'
2481    units = L[1].strip()                        #TEXT
2482
2483    L = metalines[5].strip().split()
2484    assert L[0].strip().lower() == 'spheroid'   #IGNORE
2485    spheroid = L[1].strip()                     #TEXT
2486
2487    L = metalines[6].strip().split()
2488    assert L[0].strip().lower() == 'xshift'
2489    false_easting = float(L[1].strip())
2490
2491    L = metalines[7].strip().split()
2492    assert L[0].strip().lower() == 'yshift'
2493    false_northing = float(L[1].strip())
2494
2495    #print false_easting, false_northing, zone, datum
2496
2497
2498    ###########################################
2499    #Read DEM data
2500
2501    datafile = open(basename_in + '.asc')
2502
2503    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
2504    lines = datafile.readlines()
2505    datafile.close()
2506
2507    if verbose: print 'Got', len(lines), ' lines'
2508
2509    ncols = int(lines[0].split()[1].strip())
2510    nrows = int(lines[1].split()[1].strip())
2511    xllcorner = float(lines[2].split()[1].strip())
2512    yllcorner = float(lines[3].split()[1].strip())
2513    cellsize = float(lines[4].split()[1].strip())
2514    NODATA_value = int(lines[5].split()[1].strip())
2515
2516    assert len(lines) == nrows + 6
2517
2518
2519    ##########################################
2520
2521
2522    if basename_out == None:
2523        netcdfname = root + '.dem'
2524    else:
2525        netcdfname = basename_out + '.dem'
2526
2527    if verbose: print 'Store to NetCDF file %s' %netcdfname
2528    # NetCDF file definition
2529    fid = NetCDFFile(netcdfname, 'w')
2530
2531    #Create new file
2532    fid.institution = 'Geoscience Australia'
2533    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
2534                      'of spatial point data'
2535
2536    fid.ncols = ncols
2537    fid.nrows = nrows
2538    fid.xllcorner = xllcorner
2539    fid.yllcorner = yllcorner
2540    fid.cellsize = cellsize
2541    fid.NODATA_value = NODATA_value
2542
2543    fid.zone = zone
2544    fid.false_easting = false_easting
2545    fid.false_northing = false_northing
2546    fid.projection = projection
2547    fid.datum = datum
2548    fid.units = units
2549
2550
2551    # dimension definitions
2552    fid.createDimension('number_of_rows', nrows)
2553    fid.createDimension('number_of_columns', ncols)
2554
2555    # variable definitions
2556    fid.createVariable('elevation', Float, ('number_of_rows',
2557                                            'number_of_columns'))
2558
2559    # Get handles to the variables
2560    elevation = fid.variables['elevation']
2561
2562    #Store data
2563    n = len(lines[6:])
2564    for i, line in enumerate(lines[6:]):
2565        fields = line.split()
2566        if verbose and i%((n+10)/10)==0:
2567            print 'Processing row %d of %d' %(i, nrows)
2568
2569        elevation[i, :] = array([float(x) for x in fields])
2570
2571    fid.close()
2572
2573
2574
2575
2576
2577def ferret2sww(basename_in, basename_out = None,
2578               verbose = False,
2579               minlat = None, maxlat = None,
2580               minlon = None, maxlon = None,
2581               mint = None, maxt = None, mean_stage = 0,
2582               origin = None, zscale = 1,
2583               fail_on_NaN = True,
2584               NaN_filler = 0,
2585               elevation = None,
2586               inverted_bathymetry = True
2587               ): #FIXME: Bathymetry should be obtained
2588                                  #from MOST somehow.
2589                                  #Alternatively from elsewhere
2590                                  #or, as a last resort,
2591                                  #specified here.
2592                                  #The value of -100 will work
2593                                  #for the Wollongong tsunami
2594                                  #scenario but is very hacky
2595    """Convert MOST and 'Ferret' NetCDF format for wave propagation to
2596    sww format native to abstract_2d_finite_volumes.
2597
2598    Specify only basename_in and read files of the form
2599    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
2600    relative height, x-velocity and y-velocity, respectively.
2601
2602    Also convert latitude and longitude to UTM. All coordinates are
2603    assumed to be given in the GDA94 datum.
2604
2605    min's and max's: If omitted - full extend is used.
2606    To include a value min may equal it, while max must exceed it.
2607    Lat and lon are assuemd to be in decimal degrees
2608
2609    origin is a 3-tuple with geo referenced
2610    UTM coordinates (zone, easting, northing)
2611
2612    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
2613    which means that longitude is the fastest
2614    varying dimension (row major order, so to speak)
2615
2616    ferret2sww uses grid points as vertices in a triangular grid
2617    counting vertices from lower left corner upwards, then right
2618    """
2619
2620    import os
2621    from Scientific.IO.NetCDF import NetCDFFile
2622    from Numeric import Float, Int, Int32, searchsorted, zeros, array
2623    from Numeric import allclose, around
2624
2625    precision = Float
2626
2627    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
2628
2629    if minlat != None:
2630        assert -90 < minlat < 90 , msg
2631    if maxlat != None:
2632        assert -90 < maxlat < 90 , msg
2633        if minlat != None:
2634            assert maxlat > minlat
2635    if minlon != None:
2636        assert -180 < minlon < 180 , msg
2637    if maxlon != None:
2638        assert -180 < maxlon < 180 , msg
2639        if minlon != None:
2640            assert maxlon > minlon
2641       
2642
2643
2644    #Get NetCDF data
2645    if verbose: print 'Reading files %s_*.nc' %basename_in
2646    #print "basename_in + '_ha.nc'",basename_in + '_ha.nc'
2647    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
2648    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
2649    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
2650    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
2651
2652    if basename_out is None:
2653        swwname = basename_in + '.sww'
2654    else:
2655        swwname = basename_out + '.sww'
2656
2657    # Get dimensions of file_h
2658    for dimension in file_h.dimensions.keys():
2659        if dimension[:3] == 'LON':
2660            dim_h_longitude = dimension
2661        if dimension[:3] == 'LAT':
2662            dim_h_latitude = dimension
2663        if dimension[:4] == 'TIME':
2664            dim_h_time = dimension
2665
2666#    print 'long:', dim_h_longitude
2667#    print 'lats:', dim_h_latitude
2668#    print 'times:', dim_h_time
2669
2670    times = file_h.variables[dim_h_time]
2671    latitudes = file_h.variables[dim_h_latitude]
2672    longitudes = file_h.variables[dim_h_longitude]
2673   
2674    # get dimensions for file_e
2675    for dimension in file_e.dimensions.keys():
2676        if dimension[:3] == 'LON':
2677            dim_e_longitude = dimension
2678        if dimension[:3] == 'LAT':
2679            dim_e_latitude = dimension
2680
2681    # get dimensions for file_u
2682    for dimension in file_u.dimensions.keys():
2683        if dimension[:3] == 'LON':
2684            dim_u_longitude = dimension
2685        if dimension[:3] == 'LAT':
2686            dim_u_latitude = dimension
2687        if dimension[:4] == 'TIME':
2688            dim_u_time = dimension
2689           
2690    # get dimensions for file_v
2691    for dimension in file_v.dimensions.keys():
2692        if dimension[:3] == 'LON':
2693            dim_v_longitude = dimension
2694        if dimension[:3] == 'LAT':
2695            dim_v_latitude = dimension
2696        if dimension[:4] == 'TIME':
2697            dim_v_time = dimension
2698
2699
2700    # Precision used by most for lat/lon is 4 or 5 decimals
2701    e_lat = around(file_e.variables[dim_e_latitude][:], 5)
2702    e_lon = around(file_e.variables[dim_e_longitude][:], 5)
2703
2704    # Check that files are compatible
2705    assert allclose(latitudes, file_u.variables[dim_u_latitude])
2706    assert allclose(latitudes, file_v.variables[dim_v_latitude])
2707    assert allclose(latitudes, e_lat)
2708
2709    assert allclose(longitudes, file_u.variables[dim_u_longitude])
2710    assert allclose(longitudes, file_v.variables[dim_v_longitude])
2711    assert allclose(longitudes, e_lon)
2712
2713    if mint is None:
2714        jmin = 0
2715        mint = times[0]       
2716    else:
2717        jmin = searchsorted(times, mint)
2718       
2719    if maxt is None:
2720        jmax = len(times)
2721        maxt = times[-1]
2722    else:
2723        jmax = searchsorted(times, maxt)
2724
2725    #print "latitudes[:]",latitudes[:]
2726    #print "longitudes[:]",longitudes [:]
2727    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
2728                                                  longitudes[:],
2729                                                  minlat, maxlat,
2730                                                  minlon, maxlon)
2731
2732
2733    times = times[jmin:jmax]
2734    latitudes = latitudes[kmin:kmax]
2735    longitudes = longitudes[lmin:lmax]
2736
2737    #print "latitudes[:]",latitudes[:]
2738    #print "longitudes[:]",longitudes [:]
2739
2740    if verbose: print 'cropping'
2741    zname = 'ELEVATION'
2742
2743    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
2744    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
2745    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
2746    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
2747
2748    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
2749    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
2750    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
2751    #        from Numeric import asarray
2752    #        elevations=elevations.tolist()
2753    #        elevations.reverse()
2754    #        elevations=asarray(elevations)
2755    #    else:
2756    #        from Numeric import asarray
2757    #        elevations=elevations.tolist()
2758    #        elevations.reverse()
2759    #        elevations=asarray(elevations)
2760    #        'print hmmm'
2761
2762
2763
2764    #Get missing values
2765    nan_ha = file_h.variables['HA'].missing_value[0]
2766    nan_ua = file_u.variables['UA'].missing_value[0]
2767    nan_va = file_v.variables['VA'].missing_value[0]
2768    if hasattr(file_e.variables[zname],'missing_value'):
2769        nan_e  = file_e.variables[zname].missing_value[0]
2770    else:
2771        nan_e = None
2772
2773    #Cleanup
2774    from Numeric import sometrue
2775
2776    missing = (amplitudes == nan_ha)
2777    if sometrue (missing):
2778        if fail_on_NaN:
2779            msg = 'NetCDFFile %s contains missing values'\
2780                  %(basename_in+'_ha.nc')
2781            raise DataMissingValuesError, msg
2782        else:
2783            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
2784
2785    missing = (uspeed == nan_ua)
2786    if sometrue (missing):
2787        if fail_on_NaN:
2788            msg = 'NetCDFFile %s contains missing values'\
2789                  %(basename_in+'_ua.nc')
2790            raise DataMissingValuesError, msg
2791        else:
2792            uspeed = uspeed*(missing==0) + missing*NaN_filler
2793
2794    missing = (vspeed == nan_va)
2795    if sometrue (missing):
2796        if fail_on_NaN:
2797            msg = 'NetCDFFile %s contains missing values'\
2798                  %(basename_in+'_va.nc')
2799            raise DataMissingValuesError, msg
2800        else:
2801            vspeed = vspeed*(missing==0) + missing*NaN_filler
2802
2803
2804    missing = (elevations == nan_e)
2805    if sometrue (missing):
2806        if fail_on_NaN:
2807            msg = 'NetCDFFile %s contains missing values'\
2808                  %(basename_in+'_e.nc')
2809            raise DataMissingValuesError, msg
2810        else:
2811            elevations = elevations*(missing==0) + missing*NaN_filler
2812
2813    #######
2814
2815
2816
2817    number_of_times = times.shape[0]
2818    number_of_latitudes = latitudes.shape[0]
2819    number_of_longitudes = longitudes.shape[0]
2820
2821    assert amplitudes.shape[0] == number_of_times
2822    assert amplitudes.shape[1] == number_of_latitudes
2823    assert amplitudes.shape[2] == number_of_longitudes
2824
2825    if verbose:
2826        print '------------------------------------------------'
2827        print 'Statistics:'
2828        print '  Extent (lat/lon):'
2829        print '    lat in [%f, %f], len(lat) == %d'\
2830              %(min(latitudes.flat), max(latitudes.flat),
2831                len(latitudes.flat))
2832        print '    lon in [%f, %f], len(lon) == %d'\
2833              %(min(longitudes.flat), max(longitudes.flat),
2834                len(longitudes.flat))
2835        print '    t in [%f, %f], len(t) == %d'\
2836              %(min(times.flat), max(times.flat), len(times.flat))
2837
2838        q = amplitudes.flat
2839        name = 'Amplitudes (ha) [cm]'
2840        print %s in [%f, %f]' %(name, min(q), max(q))
2841
2842        q = uspeed.flat
2843        name = 'Speeds (ua) [cm/s]'
2844        print %s in [%f, %f]' %(name, min(q), max(q))
2845
2846        q = vspeed.flat
2847        name = 'Speeds (va) [cm/s]'
2848        print %s in [%f, %f]' %(name, min(q), max(q))
2849
2850        q = elevations.flat
2851        name = 'Elevations (e) [m]'
2852        print %s in [%f, %f]' %(name, min(q), max(q))
2853
2854
2855    # print number_of_latitudes, number_of_longitudes
2856    number_of_points = number_of_latitudes*number_of_longitudes
2857    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
2858
2859
2860    file_h.close()
2861    file_u.close()
2862    file_v.close()
2863    file_e.close()
2864
2865
2866    # NetCDF file definition
2867    outfile = NetCDFFile(swwname, 'w')
2868
2869    description = 'Converted from Ferret files: %s, %s, %s, %s'\
2870                  %(basename_in + '_ha.nc',
2871                    basename_in + '_ua.nc',
2872                    basename_in + '_va.nc',
2873                    basename_in + '_e.nc')
2874   
2875    # Create new file
2876    starttime = times[0]
2877    sww = Write_sww()
2878    sww.store_header(outfile, times, number_of_volumes,
2879                     number_of_points, description=description,
2880                     verbose=verbose)
2881
2882    # Store
2883    from anuga.coordinate_transforms.redfearn import redfearn
2884    x = zeros(number_of_points, Float)  #Easting
2885    y = zeros(number_of_points, Float)  #Northing
2886
2887
2888    if verbose: print 'Making triangular grid'
2889
2890    # Check zone boundaries
2891    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
2892
2893    vertices = {}
2894    i = 0
2895    for k, lat in enumerate(latitudes):       #Y direction
2896        for l, lon in enumerate(longitudes):  #X direction
2897
2898            vertices[l,k] = i
2899
2900            zone, easting, northing = redfearn(lat,lon)
2901
2902            msg = 'Zone boundary crossed at longitude =', lon
2903            #assert zone == refzone, msg
2904            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
2905            x[i] = easting
2906            y[i] = northing
2907            i += 1
2908
2909    #Construct 2 triangles per 'rectangular' element
2910    volumes = []
2911    for l in range(number_of_longitudes-1):    #X direction
2912        for k in range(number_of_latitudes-1): #Y direction
2913            v1 = vertices[l,k+1]
2914            v2 = vertices[l,k]
2915            v3 = vertices[l+1,k+1]
2916            v4 = vertices[l+1,k]
2917
2918            volumes.append([v1,v2,v3]) #Upper element
2919            volumes.append([v4,v3,v2]) #Lower element
2920
2921    volumes = array(volumes)
2922
2923    if origin is None:
2924        origin = Geo_reference(refzone,min(x),min(y))
2925    geo_ref = write_NetCDF_georeference(origin, outfile)
2926   
2927    if elevation is not None:
2928        z = elevation
2929    else:
2930        if inverted_bathymetry:
2931            z = -1*elevations
2932        else:
2933            z = elevations
2934    #FIXME: z should be obtained from MOST and passed in here
2935
2936    from Numeric import resize
2937    z = resize(z,outfile.variables['z'][:].shape)
2938    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
2939    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
2940    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
2941    outfile.variables['elevation'][:] = z
2942    outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
2943
2944
2945
2946    #Time stepping
2947    stage = outfile.variables['stage']
2948    xmomentum = outfile.variables['xmomentum']
2949    ymomentum = outfile.variables['ymomentum']
2950
2951    if verbose: print 'Converting quantities'
2952    n = len(times)
2953    for j in range(n):
2954        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
2955        i = 0
2956        for k in range(number_of_latitudes):      #Y direction
2957            for l in range(number_of_longitudes): #X direction
2958                w = zscale*amplitudes[j,k,l]/100 + mean_stage
2959                stage[j,i] = w
2960                h = w - z[i]
2961                xmomentum[j,i] = uspeed[j,k,l]/100*h
2962                ymomentum[j,i] = vspeed[j,k,l]/100*h
2963                i += 1
2964
2965    #outfile.close()
2966
2967    #FIXME: Refactor using code from file_function.statistics
2968    #Something like print swwstats(swwname)
2969    if verbose:
2970        x = outfile.variables['x'][:]
2971        y = outfile.variables['y'][:]
2972        print '------------------------------------------------'
2973        print 'Statistics of output file:'
2974        print '  Name: %s' %swwname
2975        print '  Reference:'
2976        print '    Lower left corner: [%f, %f]'\
2977              %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
2978        print ' Start time: %f' %starttime
2979        print '    Min time: %f' %mint
2980        print '    Max time: %f' %maxt
2981        print '  Extent:'
2982        print '    x [m] in [%f, %f], len(x) == %d'\
2983              %(min(x.flat), max(x.flat), len(x.flat))
2984        print '    y [m] in [%f, %f], len(y) == %d'\
2985              %(min(y.flat), max(y.flat), len(y.flat))
2986        print '    t [s] in [%f, %f], len(t) == %d'\
2987              %(min(times), max(times), len(times))
2988        print '  Quantities [SI units]:'
2989        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2990            q = outfile.variables[name][:].flat
2991            print '    %s in [%f, %f]' %(name, min(q), max(q))
2992
2993
2994
2995    outfile.close()
2996
2997
2998
2999
3000
3001def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
3002    """Template for converting typical text files with time series to
3003    NetCDF tms file.
3004
3005
3006    The file format is assumed to be either two fields separated by a comma:
3007
3008        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
3009
3010    E.g
3011
3012      31/08/04 00:00:00, 1.328223 0 0
3013      31/08/04 00:15:00, 1.292912 0 0
3014
3015    or time (seconds), value0 value1 value2 ...
3016   
3017      0.0, 1.328223 0 0
3018      0.1, 1.292912 0 0
3019
3020    will provide a time dependent function f(t) with three attributes
3021
3022    filename is assumed to be the rootname with extenisons .txt and .sww
3023    """
3024
3025    import time, calendar
3026    from Numeric import array
3027    from anuga.config import time_format
3028    from anuga.utilities.numerical_tools import ensure_numeric
3029
3030
3031
3032    fid = open(filename + '.txt')
3033    line = fid.readline()
3034    fid.close()
3035
3036    fields = line.split(',')
3037    msg = 'File %s must have the format date, value0 value1 value2 ...'
3038    assert len(fields) == 2, msg
3039
3040    if not time_as_seconds:
3041        try:
3042            starttime = calendar.timegm(time.strptime(fields[0], time_format))
3043        except ValueError:
3044            msg = 'First field in file %s must be' %filename
3045            msg += ' date-time with format %s.\n' %time_format
3046            msg += 'I got %s instead.' %fields[0]
3047            raise DataTimeError, msg
3048    else:
3049        try:
3050            starttime = float(fields[0])
3051        except Error:
3052            msg = "Bad time format"
3053            raise DataTimeError, msg
3054
3055
3056    #Split values
3057    values = []
3058    for value in fields[1].split():
3059        values.append(float(value))
3060
3061    q = ensure_numeric(values)
3062
3063    msg = 'ERROR: File must contain at least one independent value'
3064    assert len(q.shape) == 1, msg
3065
3066
3067
3068    #Read times proper
3069    from Numeric import zeros, Float, alltrue
3070    from anuga.config import time_format
3071    import time, calendar
3072
3073    fid = open(filename + '.txt')
3074    lines = fid.readlines()
3075    fid.close()
3076
3077    N = len(lines)
3078    d = len(q)
3079
3080    T = zeros(N, Float)       #Time
3081    Q = zeros((N, d), Float)  #Values
3082
3083    for i, line in enumerate(lines):
3084        fields = line.split(',')
3085        if not time_as_seconds:
3086            realtime = calendar.timegm(time.strptime(fields[0], time_format))
3087        else:
3088             realtime = float(fields[0])
3089        T[i] = realtime - starttime
3090
3091        for j, value in enumerate(fields[1].split()):
3092            Q[i, j] = float(value)
3093
3094    msg = 'File %s must list time as a monotonuosly ' %filename
3095    msg += 'increasing sequence'
3096    assert alltrue( T[1:] - T[:-1] > 0 ), msg
3097
3098    #Create NetCDF file
3099    from Scientific.IO.NetCDF import NetCDFFile
3100
3101    fid = NetCDFFile(filename + '.tms', 'w')
3102
3103
3104    fid.institution = 'Geoscience Australia'
3105    fid.description = 'Time series'
3106
3107
3108    #Reference point
3109    #Start time in seconds since the epoch (midnight 1/1/1970)
3110    #FIXME: Use Georef
3111    fid.starttime = starttime
3112
3113    # dimension definitions
3114    #fid.createDimension('number_of_volumes', self.number_of_volumes)
3115    #fid.createDimension('number_of_vertices', 3)
3116
3117
3118    fid.createDimension('number_of_timesteps', len(T))
3119
3120    fid.createVariable('time', Float, ('number_of_timesteps',))
3121
3122    fid.variables['time'][:] = T
3123
3124    for i in range(Q.shape[1]):
3125        try:
3126            name = quantity_names[i]
3127        except:
3128            name = 'Attribute%d'%i
3129
3130        fid.createVariable(name, Float, ('number_of_timesteps',))
3131        fid.variables[name][:] = Q[:,i]
3132
3133    fid.close()
3134
3135
3136def extent_sww(file_name):
3137    """
3138    Read in an sww file.
3139
3140    Input;
3141    file_name - the sww file
3142
3143    Output;
3144    z - Vector of bed elevation
3145    volumes - Array.  Each row has 3 values, representing
3146    the vertices that define the volume
3147    time - Vector of the times where there is stage information
3148    stage - array with respect to time and vertices (x,y)
3149    """
3150
3151
3152    from Scientific.IO.NetCDF import NetCDFFile
3153
3154    #Check contents
3155    #Get NetCDF
3156    fid = NetCDFFile(file_name, 'r')
3157
3158    # Get the variables
3159    x = fid.variables['x'][:]
3160    y = fid.variables['y'][:]
3161    stage = fid.variables['stage'][:]
3162    #print "stage",stage
3163    #print "stage.shap",stage.shape
3164    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
3165    #print "min(stage)",min(stage)
3166
3167    fid.close()
3168
3169    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
3170
3171
3172def sww2domain(filename,boundary=None,t=None,\
3173               fail_if_NaN=True,NaN_filler=0\
3174               ,verbose = False,very_verbose = False):
3175    """
3176    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
3177
3178    Boundary is not recommended if domain.smooth is not selected, as it
3179    uses unique coordinates, but not unique boundaries. This means that
3180    the boundary file will not be compatable with the coordinates, and will
3181    give a different final boundary, or crash.
3182    """
3183    NaN=9.969209968386869e+036
3184    #initialise NaN.
3185
3186    from Scientific.IO.NetCDF import NetCDFFile
3187    from shallow_water import Domain
3188    from Numeric import asarray, transpose, resize
3189
3190    if verbose: print 'Reading from ', filename
3191    fid = NetCDFFile(filename, 'r')    #Open existing file for read
3192    time = fid.variables['time']       #Timesteps
3193    if t is None:
3194        t = time[-1]
3195    time_interp = get_time_interp(time,t)
3196
3197    # Get the variables as Numeric arrays
3198    x = fid.variables['x'][:]             #x-coordinates of vertices
3199    y = fid.variables['y'][:]             #y-coordinates of vertices
3200    elevation = fid.variables['elevation']     #Elevation
3201    stage = fid.variables['stage']     #Water level
3202    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
3203    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
3204
3205    starttime = fid.starttime[0]
3206    volumes = fid.variables['volumes'][:] #Connectivity
3207    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
3208
3209    conserved_quantities = []
3210    interpolated_quantities = {}
3211    other_quantities = []
3212
3213    # get geo_reference
3214    #sww files don't have to have a geo_ref
3215    try:
3216        geo_reference = Geo_reference(NetCDFObject=fid)
3217    except: #AttributeError, e:
3218        geo_reference = None
3219
3220    if verbose: print '    getting quantities'
3221    for quantity in fid.variables.keys():
3222        dimensions = fid.variables[quantity].dimensions
3223        if 'number_of_timesteps' in dimensions:
3224            conserved_quantities.append(quantity)
3225            interpolated_quantities[quantity]=\
3226                  interpolated_quantity(fid.variables[quantity][:],time_interp)
3227        else: other_quantities.append(quantity)
3228
3229    other_quantities.remove('x')
3230    other_quantities.remove('y')
3231    other_quantities.remove('z')
3232    other_quantities.remove('volumes')
3233    try:
3234        other_quantities.remove('stage_range')
3235        other_quantities.remove('xmomentum_range')
3236        other_quantities.remove('ymomentum_range')
3237        other_quantities.remove('elevation_range')
3238    except:
3239        pass
3240       
3241
3242    conserved_quantities.remove('time')
3243
3244    if verbose: print '    building domain'
3245    #    From domain.Domain:
3246    #    domain = Domain(coordinates, volumes,\
3247    #                    conserved_quantities = conserved_quantities,\
3248    #                    other_quantities = other_quantities,zone=zone,\
3249    #                    xllcorner=xllcorner, yllcorner=yllcorner)
3250
3251    #   From shallow_water.Domain:
3252    coordinates=coordinates.tolist()
3253    volumes=volumes.tolist()
3254    #FIXME:should this be in mesh?(peter row)
3255    if fid.smoothing == 'Yes': unique = False
3256    else: unique = True
3257    if unique:
3258        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
3259
3260
3261    try:
3262        domain = Domain(coordinates, volumes, boundary)
3263    except AssertionError, e:
3264        fid.close()
3265        msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
3266        raise DataDomainError, msg
3267
3268    if not boundary is None:
3269        domain.boundary = boundary
3270
3271    domain.geo_reference = geo_reference
3272
3273    domain.starttime=float(starttime)+float(t)
3274    domain.time=0.0
3275
3276    for quantity in other_quantities:
3277        try:
3278            NaN = fid.variables[quantity].missing_value
3279        except:
3280            pass #quantity has no missing_value number
3281        X = fid.variables[quantity][:]
3282        if very_verbose:
3283            print '       ',quantity
3284            print '        NaN =',NaN
3285            print '        max(X)'
3286            print '       ',max(X)
3287            print '        max(X)==NaN'
3288            print '       ',max(X)==NaN
3289            print ''
3290        if (max(X)==NaN) or (min(X)==NaN):
3291            if fail_if_NaN:
3292                msg = 'quantity "%s" contains no_data entry'%quantity
3293                raise DataMissingValuesError, msg
3294            else:
3295                data = (X<>NaN)
3296                X = (X*data)+(data==0)*NaN_filler
3297        if unique:
3298            X = resize(X,(len(X)/3,3))
3299        domain.set_quantity(quantity,X)
3300    #
3301    for quantity in conserved_quantities:
3302        try:
3303            NaN = fid.variables[quantity].missing_value
3304        except:
3305            pass #quantity has no missing_value number
3306        X = interpolated_quantities[quantity]
3307        if very_verbose:
3308            print '       ',quantity
3309            print '        NaN =',NaN
3310            print '        max(X)'
3311            print '       ',max(X)
3312            print '        max(X)==NaN'
3313            print '       ',max(X)==NaN
3314            print ''
3315        if (max(X)==NaN) or (min(X)==NaN):
3316            if fail_if_NaN:
3317                msg = 'quantity "%s" contains no_data entry'%quantity
3318                raise DataMissingValuesError, msg
3319            else:
3320                data = (X<>NaN)
3321                X = (X*data)+(data==0)*NaN_filler
3322        if unique:
3323            X = resize(X,(X.shape[0]/3,3))
3324        domain.set_quantity(quantity,X)
3325
3326    fid.close()
3327    return domain
3328
3329def interpolated_quantity(saved_quantity,time_interp):
3330
3331    #given an index and ratio, interpolate quantity with respect to time.
3332    index,ratio = time_interp
3333    Q = saved_quantity
3334    if ratio > 0:
3335        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
3336    else:
3337        q = Q[index]
3338    #Return vector of interpolated values
3339    return q
3340
3341def get_time_interp(time,t=None):
3342    #Finds the ratio and index for time interpolation.
3343    #It is borrowed from previous abstract_2d_finite_volumes code.
3344    if t is None:
3345        t=time[-1]
3346        index = -1
3347        ratio = 0.
3348    else:
3349        T = time
3350        tau = t
3351        index=0
3352        msg = 'Time interval derived from file %s [%s:%s]'\
3353            %('FIXMEfilename', T[0], T[-1])
3354        msg += ' does not match model time: %s' %tau
3355        if tau < time[0]: raise DataTimeError, msg
3356        if tau > time[-1]: raise DataTimeError, msg
3357        while tau > time[index]: index += 1
3358        while tau < time[index]: index -= 1
3359        if tau == time[index]:
3360            #Protect against case where tau == time[-1] (last time)
3361            # - also works in general when tau == time[i]
3362            ratio = 0
3363        else:
3364            #t is now between index and index+1
3365            ratio = (tau - time[index])/(time[index+1] - time[index])
3366    return (index,ratio)
3367
3368
3369def weed(coordinates,volumes,boundary = None):
3370    if type(coordinates)==ArrayType:
3371        coordinates = coordinates.tolist()
3372    if type(volumes)==ArrayType:
3373        volumes = volumes.tolist()
3374
3375    unique = False
3376    point_dict = {}
3377    same_point = {}
3378    for i in range(len(coordinates)):
3379        point = tuple(coordinates[i])
3380        if point_dict.has_key(point):
3381            unique = True
3382            same_point[i]=point
3383            #to change all point i references to point j
3384        else:
3385            point_dict[point]=i
3386            same_point[i]=point
3387
3388    coordinates = []
3389    i = 0
3390    for point in point_dict.keys():
3391        point = tuple(point)
3392        coordinates.append(list(point))
3393        point_dict[point]=i
3394        i+=1
3395
3396
3397    for volume in volumes:
3398        for i in range(len(volume)):
3399            index = volume[i]
3400            if index>-1:
3401                volume[i]=point_dict[same_point[index]]
3402
3403    new_boundary = {}
3404    if not boundary is None:
3405        for segment in boundary.keys():
3406            point0 = point_dict[same_point[segment[0]]]
3407            point1 = point_dict[same_point[segment[1]]]
3408            label = boundary[segment]
3409            #FIXME should the bounday attributes be concaterated
3410            #('exterior, pond') or replaced ('pond')(peter row)
3411
3412            if new_boundary.has_key((point0,point1)):
3413                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
3414                                              #+','+label
3415
3416            elif new_boundary.has_key((point1,point0)):
3417                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
3418                                              #+','+label
3419            else: new_boundary[(point0,point1)]=label
3420
3421        boundary = new_boundary
3422
3423    return coordinates,volumes,boundary
3424
3425
3426def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
3427                 verbose=False):
3428    """Read Digitial Elevation model from the following NetCDF format (.dem)
3429
3430    Example:
3431
3432    ncols         3121
3433    nrows         1800
3434    xllcorner     722000
3435    yllcorner     5893000
3436    cellsize      25
3437    NODATA_value  -9999
3438    138.3698 137.4194 136.5062 135.5558 ..........
3439
3440    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
3441    """
3442
3443    import os
3444    from Scientific.IO.NetCDF import NetCDFFile
3445    from Numeric import Float, zeros, sum, reshape, equal
3446
3447    root = basename_in
3448    inname = root + '.dem'
3449
3450    #Open existing netcdf file to read
3451    infile = NetCDFFile(inname, 'r')
3452    if verbose: print 'Reading DEM from %s' %inname
3453
3454    #Read metadata
3455    ncols = infile.ncols[0]
3456    nrows = infile.nrows[0]
3457    xllcorner = infile.xllcorner[0]
3458    yllcorner = infile.yllcorner[0]
3459    cellsize = infile.cellsize[0]
3460    NODATA_value = infile.NODATA_value[0]
3461    zone = infile.zone[0]
3462    false_easting = infile.false_easting[0]
3463    false_northing = infile.false_northing[0]
3464    projection = infile.projection
3465    datum = infile.datum
3466    units = infile.units
3467
3468    dem_elevation = infile.variables['elevation']
3469
3470    #Get output file name
3471    if basename_out == None:
3472        outname = root + '_' + repr(cellsize_new) + '.dem'
3473    else:
3474        outname = basename_out + '.dem'
3475
3476    if verbose: print 'Write decimated NetCDF file to %s' %outname
3477
3478    #Determine some dimensions for decimated grid
3479    (nrows_stencil, ncols_stencil) = stencil.shape
3480    x_offset = ncols_stencil / 2
3481    y_offset = nrows_stencil / 2
3482    cellsize_ratio = int(cellsize_new / cellsize)
3483    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
3484    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
3485
3486    #Open netcdf file for output
3487    outfile = NetCDFFile(outname, 'w')
3488
3489    #Create new file
3490    outfile.institution = 'Geoscience Australia'
3491    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
3492                           'of spatial point data'
3493    #Georeferencing
3494    outfile.zone = zone
3495    outfile.projection = projection
3496    outfile.datum = datum
3497    outfile.units = units
3498
3499    outfile.cellsize = cellsize_new
3500    outfile.NODATA_value = NODATA_value
3501    outfile.false_easting = false_easting
3502    outfile.false_northing = false_northing
3503
3504    outfile.xllcorner = xllcorner + (x_offset * cellsize)
3505    outfile.yllcorner = yllcorner + (y_offset * cellsize)
3506    outfile.ncols = ncols_new
3507    outfile.nrows = nrows_new
3508
3509    # dimension definition
3510    outfile.createDimension('number_of_points', nrows_new*ncols_new)
3511
3512    # variable definition
3513    outfile.createVariable('elevation', Float, ('number_of_points',))
3514
3515    # Get handle to the variable
3516    elevation = outfile.variables['elevation']
3517
3518    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
3519
3520    #Store data
3521    global_index = 0
3522    for i in range(nrows_new):
3523        if verbose: print 'Processing row %d of %d' %(i, nrows_new)
3524        lower_index = global_index
3525        telev =  zeros(ncols_new, Float)
3526        local_index = 0
3527        trow = i * cellsize_ratio
3528
3529        for j in range(ncols_new):
3530            tcol = j * cellsize_ratio
3531            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
3532
3533            #if dem contains 1 or more NODATA_values set value in
3534            #decimated dem to NODATA_value, else compute decimated
3535            #value using stencil
3536            if sum(sum(equal(tmp, NODATA_value))) > 0:
3537                telev[local_index] = NODATA_value
3538            else:
3539                telev[local_index] = sum(sum(tmp * stencil))
3540
3541            global_index += 1
3542            local_index += 1
3543
3544        upper_index = global_index
3545
3546        elevation[lower_index:upper_index] = telev
3547
3548    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
3549
3550    infile.close()
3551    outfile.close()
3552
3553
3554
3555
3556def tsh2sww(filename, verbose=False): #test_tsh2sww
3557    """
3558    to check if a tsh/msh file 'looks' good.
3559    """
3560
3561
3562    if verbose == True:print 'Creating domain from', filename
3563    domain = pmesh_to_domain_instance(filename, Domain)
3564    if verbose == True:print "Number of triangles = ", len(domain)
3565
3566    domain.smooth = True
3567    domain.format = 'sww'   #Native netcdf visualisation format
3568    file_path, filename = path.split(filename)
3569    filename, ext = path.splitext(filename)
3570    domain.set_name(filename)   
3571    domain.reduction = mean
3572    if verbose == True:print "file_path",file_path
3573    if file_path == "":file_path = "."
3574    domain.set_datadir(file_path)
3575
3576    if verbose == True:
3577        print "Output written to " + domain.get_datadir() + sep + \
3578              domain.get_name() + "." + domain.format
3579    sww = get_dataobject(domain)
3580    sww.store_connectivity()
3581    sww.store_timestep('stage')
3582
3583
3584def asc_csiro2sww(bath_dir,
3585                  elevation_dir,
3586                  ucur_dir,
3587                  vcur_dir,
3588                  sww_file,
3589                  minlat = None, maxlat = None,
3590                  minlon = None, maxlon = None,
3591                  zscale=1,
3592                  mean_stage = 0,
3593                  fail_on_NaN = True,
3594                  elevation_NaN_filler = 0,
3595                  bath_prefix='ba',
3596                  elevation_prefix='el',
3597                  verbose=False):
3598    """
3599    Produce an sww boundary file, from esri ascii data from CSIRO.
3600
3601    Also convert latitude and longitude to UTM. All coordinates are
3602    assumed to be given in the GDA94 datum.
3603
3604    assume:
3605    All files are in esri ascii format
3606
3607    4 types of information
3608    bathymetry
3609    elevation
3610    u velocity
3611    v velocity
3612
3613    Assumptions
3614    The metadata of all the files is the same
3615    Each type is in a seperate directory
3616    One bath file with extention .000
3617    The time period is less than 24hrs and uniform.
3618    """
3619    from Scientific.IO.NetCDF import NetCDFFile
3620
3621    from anuga.coordinate_transforms.redfearn import redfearn
3622
3623    precision = Float # So if we want to change the precision its done here
3624
3625    # go in to the bath dir and load the only file,
3626    bath_files = os.listdir(bath_dir)
3627
3628    bath_file = bath_files[0]
3629    bath_dir_file =  bath_dir + os.sep + bath_file
3630    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
3631
3632    #Use the date.time of the bath file as a basis for
3633    #the start time for other files
3634    base_start = bath_file[-12:]
3635
3636    #go into the elevation dir and load the 000 file
3637    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
3638                         + base_start
3639
3640    elevation_files = os.listdir(elevation_dir)
3641    ucur_files = os.listdir(ucur_dir)
3642    vcur_files = os.listdir(vcur_dir)
3643    elevation_files.sort()
3644    # the first elevation file should be the
3645    # file with the same base name as the bath data
3646    assert elevation_files[0] == 'el' + base_start
3647
3648    number_of_latitudes = bath_grid.shape[0]
3649    number_of_longitudes = bath_grid.shape[1]
3650    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3651
3652    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
3653                  for x in range(number_of_longitudes)]
3654    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
3655                 for y in range(number_of_latitudes)]
3656
3657     # reverse order of lat, so the fist lat represents the first grid row
3658    latitudes.reverse()
3659
3660    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:],
3661                                                 minlat=minlat, maxlat=maxlat,
3662                                                 minlon=minlon, maxlon=maxlon)
3663
3664
3665    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
3666    latitudes = latitudes[kmin:kmax]
3667    longitudes = longitudes[lmin:lmax]
3668    number_of_latitudes = len(latitudes)
3669    number_of_longitudes = len(longitudes)
3670    number_of_times = len(os.listdir(elevation_dir))
3671    number_of_points = number_of_latitudes*number_of_longitudes
3672    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3673
3674    #Work out the times
3675    if len(elevation_files) > 1:
3676        # Assume: The time period is less than 24hrs.
3677        time_period = (int(elevation_files[1][-3:]) - \
3678                      int(elevation_files[0][-3:]))*60*60
3679        times = [x*time_period for x in range(len(elevation_files))]
3680    else:
3681        times = [0.0]
3682
3683
3684    if verbose:
3685        print '------------------------------------------------'
3686        print 'Statistics:'
3687        print '  Extent (lat/lon):'
3688        print '    lat in [%f, %f], len(lat) == %d'\
3689              %(min(latitudes), max(latitudes),
3690                len(latitudes))
3691        print '    lon in [%f, %f], len(lon) == %d'\
3692              %(min(longitudes), max(longitudes),
3693                len(longitudes))
3694        print '    t in [%f, %f], len(t) == %d'\
3695              %(min(times), max(times), len(times))
3696
3697    ######### WRITE THE SWW FILE #############
3698    # NetCDF file definition
3699    outfile = NetCDFFile(sww_file, 'w')
3700
3701    #Create new file
3702    outfile.institution = 'Geoscience Australia'
3703    outfile.description = 'Converted from XXX'
3704
3705
3706    #For sww compatibility
3707    outfile.smoothing = 'Yes'
3708    outfile.order = 1
3709
3710    #Start time in seconds since the epoch (midnight 1/1/1970)
3711    outfile.starttime = starttime = times[0]
3712
3713
3714    # dimension definitions
3715    outfile.createDimension('number_of_volumes', number_of_volumes)
3716
3717    outfile.createDimension('number_of_vertices', 3)
3718    outfile.createDimension('number_of_points', number_of_points)
3719    outfile.createDimension('number_of_timesteps', number_of_times)
3720
3721    # variable definitions
3722    outfile.createVariable('x', precision, ('number_of_points',))
3723    outfile.createVariable('y', precision, ('number_of_points',))
3724    outfile.createVariable('elevation', precision, ('number_of_points',))
3725
3726    #FIXME: Backwards compatibility
3727    outfile.createVariable('z', precision, ('number_of_points',))
3728    #################################
3729
3730    outfile.createVariable('volumes', Int, ('number_of_volumes',
3731                                            'number_of_vertices'))
3732
3733    outfile.createVariable('time', precision,
3734                           ('number_of_timesteps',))
3735
3736    outfile.createVariable('stage', precision,
3737                           ('number_of_timesteps',
3738                            'number_of_points'))
3739
3740    outfile.createVariable('xmomentum', precision,
3741                           ('number_of_timesteps',
3742                            'number_of_points'))
3743
3744    outfile.createVariable('ymomentum', precision,
3745                           ('number_of_timesteps',
3746                            'number_of_points'))
3747
3748    #Store
3749    from anuga.coordinate_transforms.redfearn import redfearn
3750    x = zeros(number_of_points, Float)  #Easting
3751    y = zeros(number_of_points, Float)  #Northing
3752
3753    if verbose: print 'Making triangular grid'
3754    #Get zone of 1st point.
3755    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
3756
3757    vertices = {}
3758    i = 0
3759    for k, lat in enumerate(latitudes):
3760        for l, lon in enumerate(longitudes):
3761
3762            vertices[l,k] = i
3763
3764            zone, easting, northing = redfearn(lat,lon)
3765
3766            msg = 'Zone boundary crossed at longitude =', lon
3767            #assert zone == refzone, msg
3768            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
3769            x[i] = easting
3770            y[i] = northing
3771            i += 1
3772
3773
3774    #Construct 2 triangles per 'rectangular' element
3775    volumes = []
3776    for l in range(number_of_longitudes-1):    #X direction
3777        for k in range(number_of_latitudes-1): #Y direction
3778            v1 = vertices[l,k+1]
3779            v2 = vertices[l,k]
3780            v3 = vertices[l+1,k+1]
3781            v4 = vertices[l+1,k]
3782
3783            #Note, this is different to the ferrit2sww code
3784            #since the order of the lats is reversed.
3785            volumes.append([v1,v3,v2]) #Upper element
3786            volumes.append([v4,v2,v3]) #Lower element
3787
3788    volumes = array(volumes)
3789
3790    geo_ref = Geo_reference(refzone,min(x),min(y))
3791    geo_ref.write_NetCDF(outfile)
3792
3793    # This will put the geo ref in the middle
3794    #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
3795
3796
3797    if verbose:
3798        print '------------------------------------------------'
3799        print 'More Statistics:'
3800        print '  Extent (/lon):'
3801        print '    x in [%f, %f], len(lat) == %d'\
3802              %(min(x), max(x),
3803                len(x))
3804        print '    y in [%f, %f], len(lon) == %d'\
3805              %(min(y), max(y),
3806                len(y))
3807        print 'geo_ref: ',geo_ref
3808
3809    z = resize(bath_grid,outfile.variables['z'][:].shape)
3810    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
3811    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
3812    outfile.variables['z'][:] = z
3813    outfile.variables['elevation'][:] = z  #FIXME HACK
3814    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
3815
3816    stage = outfile.variables['stage']
3817    xmomentum = outfile.variables['xmomentum']
3818    ymomentum = outfile.variables['ymomentum']
3819
3820    outfile.variables['time'][:] = times   #Store time relative
3821
3822    if verbose: print 'Converting quantities'
3823    n = number_of_times
3824    for j in range(number_of_times):
3825        # load in files
3826        elevation_meta, elevation_grid = \
3827           _read_asc(elevation_dir + os.sep + elevation_files[j])
3828
3829        _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
3830        _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
3831
3832        #cut matrix to desired size
3833        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
3834        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
3835        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
3836       
3837        # handle missing values
3838        missing = (elevation_grid == elevation_meta['NODATA_value'])
3839        if sometrue (missing):
3840            if fail_on_NaN:
3841                msg = 'File %s contains missing values'\
3842                      %(elevation_files[j])
3843                raise DataMissingValuesError, msg
3844            else:
3845                elevation_grid = elevation_grid*(missing==0) + \
3846                                 missing*elevation_NaN_filler
3847
3848
3849        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
3850        i = 0
3851        for k in range(number_of_latitudes):      #Y direction
3852            for l in range(number_of_longitudes): #X direction
3853                w = zscale*elevation_grid[k,l] + mean_stage
3854                stage[j,i] = w
3855                h = w - z[i]
3856                xmomentum[j,i] = u_momentum_grid[k,l]*h
3857                ymomentum[j,i] = v_momentum_grid[k,l]*h
3858                i += 1
3859    outfile.close()
3860
3861def _get_min_max_indexes(latitudes_ref,longitudes_ref,
3862                        minlat=None, maxlat=None,
3863                        minlon=None, maxlon=None):
3864    """
3865    return max, min indexes (for slicing) of the lat and long arrays to cover the area
3866    specified with min/max lat/long
3867
3868    Think of the latitudes and longitudes describing a 2d surface.
3869    The area returned is, if possible, just big enough to cover the
3870    inputed max/min area. (This will not be possible if the max/min area
3871    has a section outside of the latitudes/longitudes area.)
3872
3873    asset  longitudes are sorted,
3874    long - from low to high (west to east, eg 148 - 151)
3875    assert latitudes are sorted, ascending or decending
3876    """
3877    latitudes = latitudes_ref[:]
3878    longitudes = longitudes_ref[:]
3879
3880    latitudes = ensure_numeric(latitudes)
3881    longitudes = ensure_numeric(longitudes)
3882   
3883    assert allclose(sort(longitudes), longitudes)
3884   
3885    lat_ascending = True
3886    if not allclose(sort(latitudes), latitudes):
3887        lat_ascending = False
3888        # reverse order of lat, so it's in ascending order         
3889        latitudes = latitudes[::-1]
3890        assert allclose(sort(latitudes), latitudes)
3891    #print "latitudes  in funct", latitudes
3892   
3893    largest_lat_index = len(latitudes)-1
3894    #Cut out a smaller extent.
3895    if minlat == None:
3896        lat_min_index = 0
3897    else:
3898        lat_min_index = searchsorted(latitudes, minlat)-1
3899        if lat_min_index <0:
3900            lat_min_index = 0
3901
3902
3903    if maxlat == None:
3904        lat_max_index = largest_lat_index #len(latitudes)
3905    else:
3906        lat_max_index = searchsorted(latitudes, maxlat)
3907        if lat_max_index > largest_lat_index:
3908            lat_max_index = largest_lat_index
3909
3910    if minlon == None:
3911        lon_min_index = 0
3912    else:
3913        lon_min_index = searchsorted(longitudes, minlon)-1
3914        if lon_min_index <0:
3915            lon_min_index = 0
3916
3917    if maxlon == None:
3918        lon_max_index = len(longitudes)
3919    else:
3920        lon_max_index = searchsorted(longitudes, maxlon)
3921
3922    # Reversing the indexes, if the lat array is decending
3923    if lat_ascending is False:
3924        lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
3925                                       largest_lat_index - lat_min_index
3926    lat_max_index = lat_max_index + 1 # taking into account how slicing works
3927    lon_max_index = lon_max_index + 1 # taking into account how slicing works
3928
3929    return lat_min_index, lat_max_index, lon_min_index, lon_max_index
3930
3931
3932def _read_asc(filename, verbose=False):
3933    """Read esri file from the following ASCII format (.asc)
3934
3935    Example:
3936
3937    ncols         3121
3938    nrows         1800
3939    xllcorner     722000
3940    yllcorner     5893000
3941    cellsize      25
3942    NODATA_value  -9999
3943    138.3698 137.4194 136.5062 135.5558 ..........
3944
3945    """
3946
3947    datafile = open(filename)
3948
3949    if verbose: print 'Reading DEM from %s' %(filename)
3950    lines = datafile.readlines()
3951    datafile.close()
3952
3953    if verbose: print 'Got', len(lines), ' lines'
3954
3955    ncols = int(lines.pop(0).split()[1].strip())
3956    nrows = int(lines.pop(0).split()[1].strip())
3957    xllcorner = float(lines.pop(0).split()[1].strip())
3958    yllcorner = float(lines.pop(0).split()[1].strip())
3959    cellsize = float(lines.pop(0).split()[1].strip())
3960    NODATA_value = float(lines.pop(0).split()[1].strip())
3961
3962    assert len(lines) == nrows
3963
3964    #Store data
3965    grid = []
3966
3967    n = len(lines)
3968    for i, line in enumerate(lines):
3969        cells = line.split()
3970        assert len(cells) == ncols
3971        grid.append(array([float(x) for x in cells]))
3972    grid = array(grid)
3973
3974    return {'xllcorner':xllcorner,
3975            'yllcorner':yllcorner,
3976            'cellsize':cellsize,
3977            'NODATA_value':NODATA_value}, grid
3978
3979
3980
3981    ####  URS 2 SWW  ###
3982
3983lon_name = 'LON'
3984lat_name = 'LAT'
3985time_name = 'TIME'
3986precision = Float # So if we want to change the precision its done here       
3987class Write_nc:
3988    """
3989    Write an nc file.
3990   
3991    Note, this should be checked to meet cdc netcdf conventions for gridded
3992    data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
3993   
3994    """
3995    def __init__(self,
3996                 quantity_name,
3997                 file_name,
3998                 time_step_count,
3999                 time_step,
4000                 lon,
4001                 lat):
4002        """
4003        time_step_count is the number of time steps.
4004        time_step is the time step size
4005       
4006        pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
4007        """
4008        self.quantity_name = quantity_name
4009        quantity_units = {'HA':'CENTIMETERS',
4010                              'UA':'CENTIMETERS/SECOND',
4011                              'VA':'CENTIMETERS/SECOND'}       
4012       
4013        multiplier_dic = {'HA':100.0, # To convert from m to cm
4014                              'UA':100.0,  #  m/s to cm/sec
4015                              'VA':-100.0}  # MUX files have positve x in the
4016        # Southern direction.  This corrects for it, when writing nc files.
4017       
4018        self.quantity_multiplier =  multiplier_dic[self.quantity_name]
4019       
4020        #self.file_name = file_name
4021        self.time_step_count = time_step_count
4022        self.time_step = time_step
4023
4024        # NetCDF file definition
4025        self.outfile = NetCDFFile(file_name, 'w')
4026        outfile = self.outfile       
4027
4028        #Create new file
4029        nc_lon_lat_header(outfile, lon, lat)
4030   
4031        # TIME
4032        outfile.createDimension(time_name, None)
4033        outfile.createVariable(time_name, precision, (time_name,))
4034
4035        #QUANTITY
4036        outfile.createVariable(self.quantity_name, precision,
4037                               (time_name, lat_name, lon_name))
4038        outfile.variables[self.quantity_name].missing_value=-1.e+034
4039        outfile.variables[self.quantity_name].units= \
4040                                 quantity_units[self.quantity_name]
4041        outfile.variables[lon_name][:]= ensure_numeric(lon)
4042        outfile.variables[lat_name][:]= ensure_numeric(lat)
4043
4044        #Assume no one will be wanting to read this, while we are writing
4045        #outfile.close()
4046       
4047    def store_timestep(self, quantity_slice):
4048        """
4049        Write a time slice of quantity info
4050        quantity_slice is the data to be stored at this time step
4051        """
4052       
4053        outfile = self.outfile
4054       
4055        # Get the variables
4056        time = outfile.variables[time_name]
4057        quantity = outfile.variables[self.quantity_name]
4058           
4059        i = len(time)
4060
4061        #Store time
4062        time[i] = i*self.time_step #self.domain.time
4063        quantity[i,:] = quantity_slice* self.quantity_multiplier
4064       
4065    def close(self):
4066        self.outfile.close()
4067
4068def urs2sww(basename_in='o', basename_out=None, verbose=False,
4069            remove_nc_files=True,
4070            minlat=None, maxlat=None,
4071            minlon= None, maxlon=None,
4072            mint=None, maxt=None,
4073            mean_stage=0,
4074            origin = None,
4075            zscale=1,
4076            fail_on_NaN=True,
4077            NaN_filler=0,
4078            elevation=None):
4079    """
4080    Convert URS C binary format for wave propagation to
4081    sww format native to abstract_2d_finite_volumes.
4082
4083    Specify only basename_in and read files of the form
4084    basefilename_velocity-z-mux, basefilename_velocity-e-mux and
4085    basefilename_waveheight-n-mux containing relative height,
4086    x-velocity and y-velocity, respectively.
4087
4088    Also convert latitude and longitude to UTM. All coordinates are
4089    assumed to be given in the GDA94 datum. The latitude and longitude
4090    information is for  a grid.
4091
4092    min's and max's: If omitted - full extend is used.
4093    To include a value min may equal it, while max must exceed it.
4094    Lat and lon are assumed to be in decimal degrees.
4095    NOTE: minlon is the most east boundary.
4096   
4097    origin is a 3-tuple with geo referenced
4098    UTM coordinates (zone, easting, northing)
4099    It will be the origin of the sww file. This shouldn't be used,
4100    since all of anuga should be able to handle an arbitary origin.
4101
4102
4103    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
4104    which means that latitude is the fastest
4105    varying dimension (row major order, so to speak)
4106
4107    In URS C binary the latitudes and longitudes are in assending order.
4108    """
4109    if basename_out == None:
4110        basename_out = basename_in
4111    files_out = urs2nc(basename_in, basename_out)
4112    ferret2sww(basename_out,
4113               minlat=minlat,
4114               maxlat=maxlat,
4115               minlon=minlon,
4116               maxlon=maxlon,
4117               mint=mint,
4118               maxt=maxt,
4119               mean_stage=mean_stage,
4120               origin=origin,
4121               zscale=zscale,
4122               fail_on_NaN=fail_on_NaN,
4123               NaN_filler=NaN_filler,
4124               inverted_bathymetry=True,
4125               verbose=verbose)
4126    #print "files_out",files_out
4127    if remove_nc_files:
4128        for file_out in files_out:
4129            os.remove(file_out)
4130   
4131def urs2nc(basename_in = 'o', basename_out = 'urs'):
4132    """
4133    Convert the 3 urs files to 4 nc files.
4134
4135    The name of the urs file names must be;
4136    [basename_in]_velocity-z-mux
4137    [basename_in]_velocity-e-mux
4138    [basename_in]_waveheight-n-mux
4139   
4140    """
4141   
4142    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4143                basename_in + EAST_VELOCITY_LABEL,
4144                basename_in + NORTH_VELOCITY_LABEL]
4145    files_out = [basename_out+'_ha.nc',
4146                 basename_out+'_ua.nc',
4147                 basename_out+'_va.nc']
4148    quantities = ['HA','UA','VA']
4149
4150    for file_name in files_in:
4151        if os.access(file_name, os.F_OK) == 0 :
4152            msg = 'File %s does not exist or is not accessible' %file_name
4153            raise IOError, msg
4154       
4155    hashed_elevation = None
4156    for file_in, file_out, quantity in map(None, files_in,
4157                                           files_out,
4158                                           quantities):
4159        lonlatdep, lon, lat, depth = _binary_c2nc(file_in,
4160                                         file_out,
4161                                         quantity)
4162        #print "lonlatdep", lonlatdep
4163        if hashed_elevation == None:
4164            elevation_file = basename_out+'_e.nc'
4165            write_elevation_nc(elevation_file,
4166                                lon,
4167                                lat,
4168                                depth)
4169            hashed_elevation = myhash(lonlatdep)
4170        else:
4171            msg = "The elevation information in the mux files is inconsistent"
4172            assert hashed_elevation == myhash(lonlatdep), msg
4173    files_out.append(elevation_file)
4174    return files_out
4175   
4176def _binary_c2nc(file_in, file_out, quantity):
4177    """
4178    Reads in a quantity urs file and writes a quantity nc file.
4179    additionally, returns the depth and lat, long info,
4180    so it can be written to a file.
4181    """
4182    columns = 3 # long, lat , depth
4183    mux_file = open(file_in, 'rb')
4184   
4185    # Number of points/stations
4186    (points_num,)= unpack('i',mux_file.read(4))
4187
4188    # nt, int - Number of time steps
4189    (time_step_count,)= unpack('i',mux_file.read(4))
4190
4191    #dt, float - time step, seconds
4192    (time_step,) = unpack('f', mux_file.read(4))
4193   
4194    msg = "Bad data in the mux file."
4195    if points_num < 0:
4196        mux_file.close()
4197        raise ANUGAError, msg
4198    if time_step_count < 0:
4199        mux_file.close()
4200        raise ANUGAError, msg
4201    if time_step < 0:
4202        mux_file.close()
4203        raise ANUGAError, msg
4204   
4205    lonlatdep = p_array.array('f')
4206    lonlatdep.read(mux_file, columns * points_num)
4207    lonlatdep = array(lonlatdep, typecode=Float)   
4208    lonlatdep = reshape(lonlatdep, (points_num, columns))
4209   
4210    lon, lat, depth = lon_lat2grid(lonlatdep)
4211    lon_sorted = list(lon)
4212    lon_sorted.sort()
4213
4214    if not lon == lon_sorted:
4215        msg = "Longitudes in mux file are not in ascending order"
4216        raise IOError, msg
4217    lat_sorted = list(lat)
4218    lat_sorted.sort()
4219
4220    if not lat == lat_sorted:
4221        msg = "Latitudes in mux file are not in ascending order"
4222   
4223    nc_file = Write_nc(quantity,
4224                       file_out,
4225                       time_step_count,
4226                       time_step,
4227                       lon,
4228                       lat)
4229
4230    for i in range(time_step_count):
4231        #Read in a time slice  from mux file 
4232        hz_p_array = p_array.array('f')
4233        hz_p_array.read(mux_file, points_num)
4234        hz_p = array(hz_p_array, typecode=Float)
4235        hz_p = reshape(hz_p, (len(lon), len(lat)))
4236        hz_p = transpose(hz_p) #mux has lat varying fastest, nc has long v.f.
4237
4238        #write time slice to nc file
4239        nc_file.store_timestep(hz_p)
4240    mux_file.close()
4241    nc_file.close()
4242
4243    return lonlatdep, lon, lat, depth
4244   
4245
4246def write_elevation_nc(file_out, lon, lat, depth_vector):
4247    """
4248    Write an nc elevation file.
4249    """
4250   
4251    # NetCDF file definition
4252    outfile = NetCDFFile(file_out, 'w')
4253
4254    #Create new file
4255    nc_lon_lat_header(outfile, lon, lat)
4256   
4257    # ELEVATION
4258    zname = 'ELEVATION'
4259    outfile.createVariable(zname, precision, (lat_name, lon_name))
4260    outfile.variables[zname].units='CENTIMETERS'
4261    outfile.variables[zname].missing_value=-1.e+034
4262
4263    outfile.variables[lon_name][:]= ensure_numeric(lon)
4264    outfile.variables[lat_name][:]= ensure_numeric(lat)
4265
4266    depth = reshape(depth_vector, ( len(lat), len(lon)))
4267    outfile.variables[zname][:]= depth
4268   
4269    outfile.close()
4270   
4271def nc_lon_lat_header(outfile, lon, lat):
4272    """
4273    outfile is the netcdf file handle.
4274    lon - a list/array of the longitudes
4275    lat - a list/array of the latitudes
4276    """
4277   
4278    outfile.institution = 'Geoscience Australia'
4279    outfile.description = 'Converted from URS binary C'
4280   
4281    # Longitude
4282    outfile.createDimension(lon_name, len(lon))
4283    outfile.createVariable(lon_name, precision, (lon_name,))
4284    outfile.variables[lon_name].point_spacing='uneven'
4285    outfile.variables[lon_name].units='degrees_east'
4286    outfile.variables[lon_name].assignValue(lon)
4287
4288
4289    # Latitude
4290    outfile.createDimension(lat_name, len(lat))
4291    outfile.createVariable(lat_name, precision, (lat_name,))
4292    outfile.variables[lat_name].point_spacing='uneven'
4293    outfile.variables[lat_name].units='degrees_north'
4294    outfile.variables[lat_name].assignValue(lat)
4295
4296
4297   
4298def lon_lat2grid(long_lat_dep):
4299    """
4300    given a list of points that are assumed to be an a grid,
4301    return the long's and lat's of the grid.
4302    long_lat_dep is an array where each row is a position.
4303    The first column is longitudes.
4304    The second column is latitudes.
4305
4306    The latitude is the fastest varying dimension - in mux files
4307    """
4308    LONG = 0
4309    LAT = 1
4310    QUANTITY = 2
4311
4312    long_lat_dep = ensure_numeric(long_lat_dep, Float)
4313   
4314    num_points = long_lat_dep.shape[0]
4315    this_rows_long = long_lat_dep[0,LONG]
4316
4317    # Count the length of unique latitudes
4318    i = 0
4319    while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
4320        i += 1
4321    # determine the lats and longsfrom the grid
4322    lat = long_lat_dep[:i, LAT]       
4323    long = long_lat_dep[::i, LONG]
4324   
4325    lenlong = len(long)
4326    lenlat = len(lat)
4327    #print 'len lat', lat, len(lat)
4328    #print 'len long', long, len(long)
4329         
4330    msg = 'Input data is not gridded'     
4331    assert num_points % lenlat == 0, msg
4332    assert num_points % lenlong == 0, msg
4333         
4334    # Test that data is gridded       
4335    for i in range(lenlong):
4336        msg = 'Data is not gridded.  It must be for this operation'
4337        first = i*lenlat
4338        last = first + lenlat
4339               
4340        assert allclose(long_lat_dep[first:last,LAT], lat), msg
4341        assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
4342   
4343   
4344#    print 'range long', min(long), max(long)
4345#    print 'range lat', min(lat), max(lat)
4346#    print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0])
4347#    print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1])
4348   
4349   
4350   
4351    msg = 'Out of range latitudes/longitudes'
4352    for l in lat:assert -90 < l < 90 , msg
4353    for l in long:assert -180 < l < 180 , msg
4354
4355    # Changing quantity from lat being the fastest varying dimension to
4356    # long being the fastest varying dimension
4357    # FIXME - make this faster/do this a better way
4358    # use numeric transpose, after reshaping the quantity vector
4359#    quantity = zeros(len(long_lat_dep), Float)
4360    quantity = zeros(num_points, Float)
4361   
4362#    print 'num',num_points
4363    for lat_i, _ in enumerate(lat):
4364        for long_i, _ in enumerate(long):
4365            q_index = lat_i*lenlong+long_i
4366            lld_index = long_i*lenlat+lat_i
4367#            print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index
4368            temp = long_lat_dep[lld_index, QUANTITY]
4369            quantity[q_index] = temp
4370           
4371    return long, lat, quantity
4372
4373    ####  END URS 2 SWW  ###
4374
4375    #### URS UNGRIDDED 2 SWW ###
4376
4377    ### PRODUCING THE POINTS NEEDED FILE ###
4378LL_LAT = -50.0
4379LL_LONG = 80.0
4380GRID_SPACING = 1.0/60.0
4381LAT_AMOUNT = 4800
4382LONG_AMOUNT = 3600
4383def URS_points_needed_to_file(file_name, boundary_polygon, zone,
4384                              ll_lat=LL_LAT, ll_long=LL_LONG,
4385                              grid_spacing=GRID_SPACING, 
4386                              lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4387                              export_csv=False, use_cache=False,
4388                              verbose=False):
4389    """
4390    file_name - name of the urs file produced for David.
4391    boundary_polygon - a list of points that describes a polygon.
4392                      The last point is assumed ot join the first point.
4393                      This is in UTM (lat long would be better though)
4394
4395    ll_lat - lower left latitude, in decimal degrees
4396    ll-long - lower left longitude, in decimal degrees
4397    grid_spacing - in deciamal degrees
4398
4399
4400    Don't add the file extension.  It will be added.
4401    """
4402    geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long,
4403                            grid_spacing, 
4404                      lat_amount, long_amount,use_cache, verbose)
4405    if not file_name[-4:] == ".urs":
4406        file_name += ".urs"
4407    geo.export_points_file(file_name)
4408    if export_csv:
4409        if file_name[-4:] == ".urs":
4410            file_name = file_name[:-4] + ".csv"
4411        geo.export_points_file(file_name)
4412
4413def URS_points_needed(boundary_polygon, zone, ll_lat=LL_LAT,
4414                      ll_long=LL_LONG, grid_spacing=GRID_SPACING, 
4415                      lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4416                      use_cache=False, verbose=False):
4417    args = (boundary_polygon,
4418                      zone)
4419    kwargs = {'ll_lat': ll_lat,
4420              'll_long': ll_long,
4421              'grid_spacing': grid_spacing,
4422              'lat_amount': lat_amount,
4423              'long_amount': long_amount} 
4424    if use_cache is True:
4425        try:
4426            from anuga.caching import cache
4427        except:
4428            msg = 'Caching was requested, but caching module'+\
4429                  'could not be imported'
4430            raise msg
4431
4432
4433        geo = cache(_URS_points_needed,
4434                  args, kwargs,
4435                  verbose=verbose,
4436                  compression=False)
4437    else:
4438        #I was getting 'got multiple values for keyword argument' errors
4439        #geo = apply(_URS_points_needed, args, kwargs)
4440        geo = _URS_points_needed(boundary_polygon,
4441                      zone, ll_lat,
4442                      ll_long, grid_spacing, 
4443                      lat_amount, long_amount)
4444
4445    return geo   
4446def _URS_points_needed(boundary_polygon,
4447                      zone, ll_lat=LL_LAT,
4448                      ll_long=LL_LONG, grid_spacing=GRID_SPACING, 
4449                      lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT):
4450    """
4451
4452    boundary_polygon - a list of points that describes a polygon.
4453                      The last point is assumed ot join the first point.
4454                      This is in UTM (lat long would be better though)
4455
4456    ll_lat - lower left latitude, in decimal degrees
4457    ll-long - lower left longitude, in decimal degrees
4458    grid_spacing - in deciamal degrees
4459
4460    """
4461   
4462    from sets import ImmutableSet
4463   
4464    msg = "grid_spacing can not be zero"
4465    assert not grid_spacing ==0, msg
4466    a = boundary_polygon
4467    # List of segments.  Each segment is two points.
4468    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
4469
4470    # convert the segs to Lat's and longs.
4471   
4472    # Don't assume the zone of the segments is the same as the lower left
4473    # corner of the lat long data!!  They can easily be in different zones
4474   
4475    lat_long_set = ImmutableSet()
4476    for seg in segs:
4477        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, 
4478                      lat_amount, long_amount, zone)
4479        lat_long_set |= ImmutableSet(points_lat_long)
4480    #print "lat_long_set",lat_long_set
4481    geo = Geospatial_data(data_points=list(lat_long_set),
4482                              points_are_lats_longs=True)
4483    return geo
4484   
4485def points_needed(seg, ll_lat, ll_long, grid_spacing, 
4486                  lat_amount, long_amount, zone):
4487    """
4488    return a list of the points, in lats and longs that are needed to
4489    interpolate any point on the segment.
4490    """
4491    from math import sqrt
4492    #print "zone",zone
4493    geo_reference = Geo_reference(zone=zone)
4494    #print "seg",seg
4495    geo = Geospatial_data(seg,geo_reference=geo_reference)
4496    seg_lat_long = geo.get_data_points(as_lat_long=True)
4497    #print "seg_lat_long", seg_lat_long
4498    # 1.415 = 2^0.5, rounded up....
4499    sqrt_2_rounded_up = 1.415
4500    buffer = sqrt_2_rounded_up * grid_spacing
4501   
4502    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
4503    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
4504    min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
4505    min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
4506
4507    first_row = (min_long - ll_long)/grid_spacing
4508    # To round up
4509    first_row_long = int(round(first_row + 0.5))
4510    #print "first_row", first_row_long
4511
4512    last_row = (max_long - ll_long)/grid_spacing # round down
4513    last_row_long = int(round(last_row))
4514    #print "last_row",last_row _long
4515   
4516    first_row = (min_lat - ll_lat)/grid_spacing
4517    # To round up
4518    first_row_lat = int(round(first_row + 0.5))
4519    #print "first_row", first_row_lat
4520
4521    last_row = (max_lat - ll_lat)/grid_spacing # round down
4522    last_row_lat = int(round(last_row))
4523    #print "last_row",last_row_lat
4524
4525    # to work out the max distance -
4526    # 111120 - horizontal distance between 1 deg latitude.
4527    #max_distance = sqrt_2_rounded_up * 111120 * grid_spacing
4528    max_distance = 157147.4112 * grid_spacing
4529    #print "max_distance", max_distance #2619.12 m for 1 minute
4530    points_lat_long = []
4531    # Create a list of the lat long points to include.
4532    for index_lat in range(first_row_lat, last_row_lat + 1):
4533        for index_long in range(first_row_long, last_row_long + 1):
4534            lat = ll_lat + index_lat*grid_spacing
4535            long = ll_long + index_long*grid_spacing
4536
4537            #filter here to keep good points
4538            if keep_point(lat, long, seg, max_distance):
4539                points_lat_long.append((lat, long)) #must be hashable
4540    #print "points_lat_long", points_lat_long
4541
4542    # Now that we have these points, lets throw ones out that are too far away
4543    return points_lat_long
4544
4545def keep_point(lat, long, seg, max_distance):
4546    """
4547    seg is two points, UTM
4548    """
4549    from math import sqrt
4550    _ , x0, y0 = redfearn(lat, long)
4551    x1 = seg[0][0]
4552    y1 = seg[0][1]
4553    x2 = seg[1][0]
4554    y2 = seg[1][1]
4555
4556    x2_1 = x2-x1
4557    y2_1 = y2-y1
4558    d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1))
4559    if d <= max_distance:
4560        return True
4561    else:
4562        return False
4563   
4564    #### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE ####
4565   
4566WAVEHEIGHT_MUX_LABEL = '_waveheight-z-mux'
4567EAST_VELOCITY_LABEL =  '_velocity-e-mux'
4568NORTH_VELOCITY_LABEL =  '_velocity-n-mux' 
4569def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
4570                      mint=None, maxt=None,
4571                      mean_stage=0,
4572                      origin=None,
4573                      hole_points_UTM=None,
4574                      zscale=1):
4575    """   
4576    Convert URS C binary format for wave propagation to
4577    sww format native to abstract_2d_finite_volumes.
4578
4579
4580    Specify only basename_in and read files of the form
4581    basefilename_velocity-z-mux, basefilename_velocity-e-mux and
4582    basefilename_waveheight-n-mux containing relative height,
4583    x-velocity and y-velocity, respectively.
4584
4585    Also convert latitude and longitude to UTM. All coordinates are
4586    assumed to be given in the GDA94 datum. The latitude and longitude
4587    information is assumed ungridded grid.
4588
4589    min's and max's: If omitted - full extend is used.
4590    To include a value min ans max may equal it.
4591    Lat and lon are assumed to be in decimal degrees.
4592   
4593    origin is a 3-tuple with geo referenced
4594    UTM coordinates (zone, easting, northing)
4595    It will be the origin of the sww file. This shouldn't be used,
4596    since all of anuga should be able to handle an arbitary origin.
4597    The mux point info is NOT relative to this origin.
4598
4599
4600    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
4601    which means that latitude is the fastest
4602    varying dimension (row major order, so to speak)
4603
4604    In URS C binary the latitudes and longitudes are in assending order.
4605
4606    Note, interpolations of the resulting sww file will be different
4607    from results of urs2sww.  This is due to the interpolation
4608    function used, and the different grid structure between urs2sww
4609    and this function.
4610   
4611    Interpolating data that has an underlying gridded source can
4612    easily end up with different values, depending on the underlying
4613    mesh.
4614
4615    consider these 4 points
4616    50  -50
4617
4618    0     0
4619
4620    The grid can be
4621     -
4622    |\|    A
4623     -
4624     or;
4625      -
4626     |/|   B
4627      -
4628      If a point is just below the center of the midpoint, it will have a
4629      +ve value in grid A and a -ve value in grid B.
4630    """ 
4631    from anuga.pmesh.mesh import Mesh, NoTrianglesError
4632
4633    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4634                basename_in + EAST_VELOCITY_LABEL,
4635                basename_in + NORTH_VELOCITY_LABEL]
4636    quantities = ['HA','UA','VA']
4637
4638    # instanciate urs_points of the three mux files.
4639    mux = {}
4640    for quantity, file in map(None, quantities, files_in):
4641        mux[quantity] = Urs_points(file)
4642       
4643    # Could check that the depth is the same. (hashing)
4644
4645    # handle to a mux file to do depth stuff
4646    a_mux = mux[quantities[0]]
4647   
4648    # Convert to utm
4649    lat = a_mux.lonlatdep[:,1]
4650    long = a_mux.lonlatdep[:,0]
4651    points_utm, zone = convert_from_latlon_to_utm( \
4652        latitudes=lat, longitudes=long)
4653    #print "points_utm", points_utm
4654    #print "zone", zone
4655
4656    elevation = a_mux.lonlatdep[:,2] * -1 #
4657   
4658    # grid ( create a mesh from the selected points)
4659    # This mesh has a problem.  Triangles are streched over ungridded areas.
4660    #  If these areas could be described as holes in pmesh, that would be great
4661
4662    # I can't just get the user to selection a point in the middle.
4663    # A boundary is needed around these points.
4664    # But if the zone of points is obvious enough auto-segment should do
4665    # a good boundary.
4666    mesh = Mesh()
4667    mesh.add_vertices(points_utm)
4668    mesh.auto_segment(smooth_indents=True, expand_pinch=True)
4669    # To try and avoid alpha shape 'hugging' too much
4670    mesh.auto_segment( mesh.shape.get_alpha()*1.1 )
4671    if hole_points_UTM is not None:
4672        point = ensure_absolute(hole_points_UTM)
4673        mesh.add_hole(point[0], point[1])
4674    try:
4675        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4676    except NoTrianglesError:
4677        # This is a bit of a hack, going in and changing the
4678        # data structure.
4679        mesh.holes = []
4680        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4681    mesh_dic = mesh.Mesh2MeshList()
4682
4683    #mesh.export_mesh_file(basename_in + '_168.tsh')
4684    #import sys; sys.exit()
4685    # These are the times of the mux file
4686    mux_times = []
4687    for i in range(a_mux.time_step_count):
4688        mux_times.append(a_mux.time_step * i) 
4689    mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt)
4690    times = mux_times[mux_times_start_i:mux_times_fin_i]
4691   
4692    if mux_times_start_i == mux_times_fin_i:
4693        # Close the mux files
4694        for quantity, file in map(None, quantities, files_in):
4695            mux[quantity].close()
4696        msg="Due to mint and maxt there's no time info in the boundary SWW."
4697        raise Exception, msg
4698       
4699    # If this raise is removed there is currently no downstream errors
4700           
4701    points_utm=ensure_numeric(points_utm)
4702    #print "mesh_dic['generatedpointlist']", mesh_dic['generatedpointlist']
4703    #print "points_utm", points_utm
4704    assert ensure_numeric(mesh_dic['generatedpointlist']) == \
4705           ensure_numeric(points_utm)
4706   
4707    volumes = mesh_dic['generatedtrianglelist']
4708   
4709    # write sww intro and grid stuff.   
4710    if basename_out is None:
4711        swwname = basename_in + '.sww'
4712    else:
4713        swwname = basename_out + '.sww'
4714
4715    if verbose: print 'Output to ', swwname
4716    outfile = NetCDFFile(swwname, 'w')
4717    # For a different way of doing this, check out tsh2sww
4718    # work out sww_times and the index range this covers
4719    sww = Write_sww()
4720    sww.store_header(outfile, times, len(volumes), len(points_utm),
4721                     verbose=verbose)
4722    outfile.mean_stage = mean_stage
4723    outfile.zscale = zscale
4724
4725    sww.store_triangulation(outfile, points_utm, volumes,
4726                            elevation, zone,  new_origin=origin,
4727                            verbose=verbose)
4728   
4729    if verbose: print 'Converting quantities'
4730    j = 0
4731    # Read in a time slice from each mux file and write it to the sww file
4732    for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
4733        if j >= mux_times_start_i and j < mux_times_fin_i:
4734            stage = zscale*ha + mean_stage
4735            h = stage - elevation
4736            xmomentum = ua*h
4737            ymomentum = -1*va*h # -1 since in mux files south is positive.
4738            sww.store_quantities(outfile, 
4739                                 slice_index=j - mux_times_start_i,
4740                                 verbose=verbose,
4741                                 stage=stage,
4742                                 xmomentum=xmomentum,
4743                                 ymomentum=ymomentum)
4744        j += 1
4745    if verbose: sww.verbose_quantities(outfile)
4746    outfile.close()
4747    #
4748    # Do some conversions while writing the sww file
4749
4750   
4751def mux2sww_time(mux_times, mint, maxt):
4752    """
4753    """
4754
4755    if mint == None:
4756        mux_times_start_i = 0
4757    else:
4758        mux_times_start_i = searchsorted(mux_times, mint)
4759       
4760    if maxt == None:
4761        mux_times_fin_i = len(mux_times)
4762    else:
4763        maxt += 0.5 # so if you specify a time where there is
4764                    # data that time will be included
4765        mux_times_fin_i = searchsorted(mux_times, maxt)
4766
4767    return mux_times_start_i, mux_times_fin_i
4768
4769
4770class Write_sww:
4771    from anuga.shallow_water.shallow_water_domain import Domain
4772
4773    # FIXME (Ole): Hardwiring the conserved quantities like
4774    # this could be a problem. I would prefer taking them from
4775    # the instantiation of Domain.
4776    sww_quantities = Domain.conserved_quantities
4777
4778
4779    RANGE = '_range'
4780    EXTREMA = ':extrema'
4781
4782    def __init__(self):
4783        pass
4784   
4785    def store_header(self,
4786                     outfile,
4787                     times,
4788                     number_of_volumes,
4789                     number_of_points,
4790                     description='Converted from XXX',
4791                     smoothing=True,
4792                     order=1, verbose=False):
4793        """
4794        outfile - the name of the file that will be written
4795        times - A list of the time slice times OR a start time
4796        Note, if a list is given the info will be made relative.
4797        number_of_volumes - the number of triangles
4798        """
4799   
4800        outfile.institution = 'Geoscience Australia'
4801        outfile.description = description
4802
4803        # For sww compatibility
4804        if smoothing is True:
4805            # Smoothing to be depreciated
4806            outfile.smoothing = 'Yes'
4807            outfile.vertices_are_stored_uniquely = 'False'
4808        else:
4809            # Smoothing to be depreciated
4810            outfile.smoothing = 'No'
4811            outfile.vertices_are_stored_uniquely = 'True'
4812        outfile.order = order
4813
4814        try:
4815            revision_number = get_revision_number()
4816        except:
4817            revision_number = None
4818        # Allow None to be stored as a string               
4819        outfile.revision_number = str(revision_number) 
4820
4821
4822       
4823        # times - A list or array of the time slice times OR a start time
4824        # times = ensure_numeric(times)
4825        # Start time in seconds since the epoch (midnight 1/1/1970)
4826
4827        # This is being used to seperate one number from a list.
4828        # what it is actually doing is sorting lists from numeric arrays.
4829        if type(times) is list or type(times) is ArrayType: 
4830            number_of_times = len(times)
4831            times = ensure_numeric(times) 
4832            if number_of_times == 0:
4833                starttime = 0
4834            else:
4835                starttime = times[0]
4836                times = times - starttime  #Store relative times
4837        else:
4838            number_of_times = 0
4839            starttime = times
4840            #times = ensure_numeric([])
4841        outfile.starttime = starttime
4842        # dimension definitions
4843        outfile.createDimension('number_of_volumes', number_of_volumes)
4844        outfile.createDimension('number_of_vertices', 3)
4845        outfile.createDimension('numbers_in_range', 2)
4846   
4847        if smoothing is True:
4848            outfile.createDimension('number_of_points', number_of_points)
4849       
4850            # FIXME(Ole): This will cause sww files for paralle domains to
4851            # have ghost nodes stored (but not used by triangles).
4852            # To clean this up, we have to change get_vertex_values and
4853            # friends in quantity.py (but I can't be bothered right now)
4854        else:
4855            outfile.createDimension('number_of_points', 3*number_of_volumes)
4856        outfile.createDimension('number_of_timesteps', number_of_times)
4857
4858        # variable definitions
4859        outfile.createVariable('x', precision, ('number_of_points',))
4860        outfile.createVariable('y', precision, ('number_of_points',))
4861        outfile.createVariable('elevation', precision, ('number_of_points',))
4862        q = 'elevation'
4863        outfile.createVariable(q+Write_sww.RANGE, precision,
4864                               ('numbers_in_range',))
4865
4866
4867        # Initialise ranges with small and large sentinels.
4868        # If this was in pure Python we could have used None sensibly
4869        outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min               
4870        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
4871
4872        # FIXME: Backwards compatibility
4873        outfile.createVariable('z', precision, ('number_of_points',))
4874        #################################
4875
4876        outfile.createVariable('volumes', Int, ('number_of_volumes',
4877                                                'number_of_vertices'))
4878       
4879        outfile.createVariable('time', precision,
4880                               ('number_of_timesteps',))
4881       
4882        for q in Write_sww.sww_quantities:
4883            outfile.createVariable(q, precision,
4884                                   ('number_of_timesteps',
4885                                    'number_of_points')) 
4886            outfile.createVariable(q+Write_sww.RANGE, precision,
4887                                   ('numbers_in_range',))
4888
4889            # Initialise ranges with small and large sentinels.
4890            # If this was in pure Python we could have used None sensibly
4891            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min               
4892            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
4893           
4894        if type(times) is list or type(times) is ArrayType: 
4895            outfile.variables['time'][:] = times    #Store time relative
4896           
4897        if verbose:
4898            print '------------------------------------------------'
4899            print 'Statistics:'
4900            print '    t in [%f, %f], len(t) == %d'\
4901                  %(min(times.flat), max(times.flat), len(times.flat))
4902
4903       
4904    def store_triangulation(self,
4905                            outfile,
4906                            points_utm,
4907                            volumes,
4908                            elevation, zone=None, new_origin=None, 
4909                            points_georeference=None, verbose=False):
4910        """
4911       
4912        new_origin - qa georeference that the points can be set to. (Maybe
4913        do this before calling this function.)
4914       
4915        points_utm - currently a list or array of the points in UTM.
4916        points_georeference - the georeference of the points_utm
4917       
4918        How about passing new_origin and current_origin.
4919        If you get both, do a convertion from the old to the new.
4920       
4921        If you only get new_origin, the points are absolute,
4922        convert to relative
4923       
4924        if you only get the current_origin the points are relative, store
4925        as relative.
4926       
4927        if you get no georefs create a new georef based on the minimums of
4928        points_utm.  (Another option would be to default to absolute)
4929       
4930        Yes, and this is done in another part of the code.
4931        Probably geospatial.
4932       
4933        If you don't supply either geo_refs, then supply a zone. If not
4934        the default zone will be used.
4935       
4936       
4937        precon
4938       
4939        header has been called.
4940        """
4941       
4942        number_of_points = len(points_utm)   
4943        volumes = array(volumes) 
4944        points_utm = array(points_utm)
4945
4946        # given the two geo_refs and the points, do the stuff
4947        # described in the method header
4948        # if this is needed else where, pull out as a function
4949        points_georeference = ensure_geo_reference(points_georeference)
4950        new_origin = ensure_geo_reference(new_origin)
4951        if new_origin is None and points_georeference is not None:
4952            points = points_utm
4953            geo_ref = points_georeference
4954        else:
4955            if new_origin is None:
4956                new_origin = Geo_reference(zone,min(points_utm[:,0]),
4957                                           min(points_utm[:,1]))
4958            points = new_origin.change_points_geo_ref(points_utm,
4959                                                      points_georeference)
4960            geo_ref = new_origin
4961
4962        # At this stage I need a georef and points
4963        # the points are relative to the georef
4964        geo_ref.write_NetCDF(outfile)
4965   
4966        # This will put the geo ref in the middle
4967        #geo_ref=Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
4968       
4969        x =  points[:,0]
4970        y =  points[:,1]
4971        z = outfile.variables['z'][:]
4972   
4973        if verbose:
4974            print '------------------------------------------------'
4975            print 'More Statistics:'
4976            print '  Extent (/lon):'
4977            print '    x in [%f, %f], len(lat) == %d'\
4978                  %(min(x), max(x),
4979                    len(x))
4980            print '    y in [%f, %f], len(lon) == %d'\
4981                  %(min(y), max(y),
4982                    len(y))
4983            print '    z in [%f, %f], len(z) == %d'\
4984                  %(min(elevation), max(elevation),
4985                    len(elevation))
4986            print 'geo_ref: ',geo_ref
4987            print '------------------------------------------------'
4988           
4989        #z = resize(bath_grid,outfile.variables['z'][:].shape)
4990        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
4991        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
4992        outfile.variables['z'][:] = elevation
4993        outfile.variables['elevation'][:] = elevation  #FIXME HACK
4994        outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
4995
4996        q = 'elevation'
4997        # This updates the _range values
4998        outfile.variables[q+Write_sww.RANGE][0] = min(elevation)
4999        outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
5000
5001
5002    def store_quantities(self, outfile, precision=Float,
5003                         slice_index=None, time=None,
5004                         verbose=False, **quant):
5005        """
5006        Write the quantity info.
5007
5008        **quant is extra keyword arguments passed in. These must be
5009          the sww quantities, currently; stage, xmomentum, ymomentum.
5010       
5011        if the time array is already been built, use the slice_index
5012        to specify the index.
5013       
5014        Otherwise, use time to increase the time dimension
5015
5016        Maybe make this general, but the viewer assumes these quantities,
5017        so maybe we don't want it general - unless the viewer is general
5018       
5019        precon
5020        triangulation and
5021        header have been called.
5022        """
5023
5024        if time is not None:
5025            file_time = outfile.variables['time']
5026            slice_index = len(file_time)
5027            file_time[slice_index] = time   
5028
5029        # write the conserved quantities from Domain.
5030        # Typically stage,  xmomentum, ymomentum
5031        # other quantities will be ignored, silently.
5032        for q in Write_sww.sww_quantities:
5033            if not quant.has_key(q):
5034                msg = 'SWW file can not write quantity %s' %q
5035                raise NewQuantity, msg
5036            else:
5037                q_values = quant[q]
5038                outfile.variables[q][slice_index] = q_values.astype(precision)
5039
5040                # This updates the _range values
5041                q_range = outfile.variables[q+Write_sww.RANGE][:]
5042                q_values_min = min(q_values)
5043                if q_values_min < q_range[0]:
5044                    outfile.variables[q+Write_sww.RANGE][0] = q_values_min
5045                q_values_max = max(q_values)
5046                if q_values_max > q_range[1]:
5047                    outfile.variables[q+Write_sww.RANGE][1] = q_values_max
5048
5049    def verbose_quantities(self, outfile):
5050        print '------------------------------------------------'
5051        print 'More Statistics:'
5052        for q in Write_sww.sww_quantities:
5053            print %s in [%f, %f]' %(q,
5054                                       outfile.variables[q+Write_sww.RANGE][0],
5055                                       outfile.variables[q+Write_sww.RANGE][1])
5056        print '------------------------------------------------'
5057
5058
5059       
5060def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation,
5061                         mean_stage=0, zscale=1,
5062                         verbose=False):   
5063    #Time stepping
5064    stage = outfile.variables['stage']
5065    xmomentum = outfile.variables['xmomentum']
5066    ymomentum = outfile.variables['ymomentum']
5067
5068    n = len(has)
5069    j=0
5070    for ha, ua, va in map(None, has, uas, vas):
5071        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
5072        w = zscale*ha + mean_stage
5073        stage[j] = w
5074        h = w - elevation
5075        xmomentum[j] = ua*h
5076        ymomentum[j] = -1*va*#  -1 since in mux files south is positive.
5077        j += 1
5078   
5079def urs2txt(basename_in, location_index=None):
5080    """
5081    Not finished or tested
5082    """
5083   
5084    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
5085                basename_in + EAST_VELOCITY_LABEL,
5086                basename_in + NORTH_VELOCITY_LABEL]
5087    quantities = ['HA','UA','VA']
5088
5089    d = ","
5090   
5091    # instanciate urs_points of the three mux files.
5092    mux = {}
5093    for quantity, file in map(None, quantities, files_in):
5094        mux[quantity] = Urs_points(file)
5095       
5096    # Could check that the depth is the same. (hashing)
5097
5098    # handle to a mux file to do depth stuff
5099    a_mux = mux[quantities[0]]
5100   
5101    # Convert to utm
5102    latitudes = a_mux.lonlatdep[:,1]
5103    longitudes = a_mux.lonlatdep[:,0]
5104    points_utm, zone = convert_from_latlon_to_utm( \
5105        latitudes=latitudes, longitudes=longitudes)
5106    #print "points_utm", points_utm
5107    #print "zone", zone
5108    depths = a_mux.lonlatdep[:,2]  #
5109   
5110    fid = open(basename_in + '.txt', 'w')
5111
5112    fid.write("zone: " + str(zone) + "\n")
5113
5114    if location_index is not None:
5115        #Title
5116        li = location_index
5117        fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5118                  'Northing' + "\n")
5119        fid.write(str(li) +d+ str(latitudes[li])+d+ \
5120              str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \
5121              str(points_utm[li][01]) + "\n")
5122
5123    # the non-time dependent stuff
5124    #Title
5125    fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5126                  'Northing' +d+ 'depth m' + "\n")
5127    i = 0
5128    for depth, point_utm, lat, long in map(None, depths,
5129                                               points_utm, latitudes,
5130                                               longitudes):
5131       
5132        fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \
5133                  str(point_utm[01]) +d+ str(depth) + "\n")
5134        i +=1
5135    #Time dependent
5136    if location_index is not None:
5137        time_step = a_mux.time_step
5138        i = 0
5139        #Title
5140        fid.write('time' +d+ 'HA depth m'+d+ \
5141                 'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \
5142                      + "\n")
5143        for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
5144            fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \
5145                      str(UA[location_index]) +d+ str(VA[location_index]) \
5146                      + "\n")
5147           
5148            i +=1
5149   
5150class Urs_points:
5151    """
5152    Read the info in URS mux files.
5153
5154    for the quantities heres a correlation between the file names and
5155    what they mean;
5156    z-mux is height above sea level, m
5157    e-mux is velocity is Eastern direction, m/s
5158    n-mux is velocity is Northern direction, m/s   
5159    """
5160    def __init__(self,urs_file):
5161        self.iterated = False
5162        columns = 3 # long, lat , depth
5163        mux_file = open(urs_file, 'rb')
5164       
5165        # Number of points/stations
5166        (self.points_num,)= unpack('i',mux_file.read(4))
5167       
5168        # nt, int - Number of time steps
5169        (self.time_step_count,)= unpack('i',mux_file.read(4))
5170        #print "self.time_step_count", self.time_step_count
5171        #dt, float - time step, seconds
5172        (self.time_step,) = unpack('f', mux_file.read(4))
5173        #print "self.time_step", self.time_step
5174        msg = "Bad data in the urs file."
5175        if self.points_num < 0:
5176            mux_file.close()
5177            raise ANUGAError, msg
5178        if self.time_step_count < 0:
5179            mux_file.close()
5180            raise ANUGAError, msg
5181        if self.time_step < 0:
5182            mux_file.close()
5183            raise ANUGAError, msg
5184
5185        # the depth is in meters, and it is the distance from the ocean
5186        # to the sea bottom.
5187        lonlatdep = p_array.array('f')
5188        lonlatdep.read(mux_file, columns * self.points_num)
5189        lonlatdep = array(lonlatdep, typecode=Float)   
5190        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
5191        #print 'lonlatdep',lonlatdep
5192        self.lonlatdep = lonlatdep
5193       
5194        self.mux_file = mux_file
5195        # check this array
5196
5197    def __iter__(self):
5198        """
5199        iterate over quantity data which is with respect to time.
5200
5201        Note: You can only interate once over an object
5202       
5203        returns quantity infomation for each time slice
5204        """
5205        msg =  "You can only interate once over a urs file."
5206        assert not self.iterated, msg
5207        self.iter_time_step = 0
5208        self.iterated = True
5209        return self
5210   
5211    def next(self):
5212        if self.time_step_count == self.iter_time_step:
5213            self.close()
5214            raise StopIteration
5215        #Read in a time slice  from mux file 
5216        hz_p_array = p_array.array('f')
5217        hz_p_array.read(self.mux_file, self.points_num)
5218        hz_p = array(hz_p_array, typecode=Float)
5219        self.iter_time_step += 1
5220       
5221        return hz_p
5222
5223    def close(self):
5224        self.mux_file.close()
5225       
5226    #### END URS UNGRIDDED 2 SWW ###
5227
5228       
5229def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
5230                         print_to_screen=False, verbose=False):
5231    """
5232    Used to store screen output and errors to file, if run on multiple
5233    processes eachprocessor will have its own output and error file.
5234   
5235    extra_info - is used as a string that can identify outputs with another
5236    string eg. '_other'
5237    """
5238    import sys
5239    dir_name = dir_name
5240    if access(dir_name,W_OK) == 0:
5241        if verbose: print 'Make directory %s' %dir_name
5242        if verbose: print "myid", myid
5243        mkdir (dir_name,0777)
5244    if myid <>'':
5245        myid = '_'+str(myid)
5246    if numprocs <>'':
5247        numprocs = '_'+str(numprocs)
5248    if extra_info <>'':
5249        extra_info = '_'+str(extra_info)
5250    screen_output_name = dir_name + "screen_output%s%s%s.txt" %(myid,
5251                                                                numprocs,
5252                                                                extra_info)
5253    screen_error_name = dir_name + "screen_error%s%s%s.txt" %(myid,
5254                                                              numprocs,
5255                                                              extra_info)
5256    print screen_output_name
5257    #used to catch screen output to file
5258    sys.stdout = Screen_Catcher(screen_output_name)
5259    sys.stderr = Screen_Catcher(screen_error_name)
5260
5261class Screen_Catcher:
5262    """this simply catches the screen output and stores it to file defined by
5263    start_screen_catcher (above)
5264    """
5265   
5266    def __init__(self, filename):
5267        self.filename = filename
5268 
5269        if exists(self.filename)is True:
5270            print'Old existing file "%s" has been deleted' %(self.filename)
5271            remove(self.filename)
5272
5273    def write(self, stuff):
5274        fid = open(self.filename, 'a')
5275        fid.write(stuff)
5276#        if print_to_screen: print stuff
5277
5278def copy_code_files(dir_name, filename1, filename2):
5279    """Copies "filename1" and "filename2" to "dir_name". Very useful for
5280    information management
5281    filename1 and filename2 are both absolute pathnames   
5282    """
5283
5284    if access(dir_name,F_OK) == 0:
5285        print 'Make directory %s' %dir_name
5286        mkdir (dir_name,0777)
5287    shutil.copy(filename1, dir_name + sep + basename(filename1))
5288    shutil.copy(filename2, dir_name + sep + basename(filename2))
5289#    copy (__file__, project.output_run_time_dir + basename(__file__))
5290    print 'Files %s and %s copied' %(filename1, filename2)
5291
5292def get_data_from_file(filename,separator_value = ','):
5293    """
5294    Read in data information from file and
5295   
5296    Returns:
5297        header_fields, a string? of the first line separated
5298        by the 'separator_value'
5299       
5300        data, a array (N data columns X M lines) in the file
5301        excluding the header
5302       
5303    NOTE: wont deal with columns with different lenghts and there must be
5304    no blank lines at the end.
5305    """
5306   
5307    fid = open(filename)
5308    lines = fid.readlines()
5309   
5310    fid.close()
5311   
5312    header_line = lines[0]
5313    header_fields = header_line.split(separator_value)
5314
5315    #array to store data, number in there is to allow float...
5316    #i'm sure there is a better way!
5317    data=array([],typecode=Float)
5318    data=resize(data,((len(lines)-1),len(header_fields)))
5319#    print 'number of fields',range(len(header_fields))
5320#    print 'number of lines',len(lines), shape(data)
5321#    print'data',data[1,1],header_line
5322
5323    array_number = 0
5324    line_number = 1
5325    while line_number < (len(lines)):
5326        for i in range(len(header_fields)): 
5327            #this get line below the header, explaining the +1
5328            #and also the line_number can be used as the array index
5329            fields = lines[line_number].split(separator_value)
5330            #assign to array
5331            data[array_number,i] = float(fields[i])
5332           
5333        line_number = line_number +1
5334        array_number = array_number +1
5335       
5336    return header_fields, data
5337
5338def store_parameters(verbose=False,**kwargs):
5339    """
5340    Store "kwargs" into a temp csv file, if "completed" is a kwargs csv file is
5341    kwargs[file_name] else it is kwargs[output_dir] + details_temp.csv
5342   
5343    Must have a file_name keyword arg, this is what is writing to.
5344    might be a better way to do this using CSV module Writer and writeDict
5345   
5346    writes file to "output_dir" unless "completed" is in kwargs, then
5347    it writes to "file_name" kwargs
5348
5349    """
5350    import types
5351#    import os
5352   
5353    # Check that kwargs is a dictionary
5354    if type(kwargs) != types.DictType:
5355        raise TypeError
5356   
5357    #is completed is kwargs?
5358    try:
5359        kwargs['completed']
5360        completed=True
5361    except:
5362        completed=False
5363 
5364    #get file name and removes from dict and assert that a file_name exists
5365    if completed:
5366        try:
5367            file = str(kwargs['file_name'])
5368        except:
5369            raise 'kwargs must have file_name'
5370    else:
5371        #write temp file in output directory
5372        try:
5373            file = str(kwargs['output_dir'])+'detail_temp.csv'
5374        except:
5375            raise 'kwargs must have output_dir'
5376       
5377    #extracts the header info and the new line info
5378    line=''
5379    header=''
5380    count=0
5381    keys = kwargs.keys()
5382    keys.sort()
5383   
5384    #used the sorted keys to create the header and line data
5385    for k in keys:
5386#        print "%s = %s" %(k, kwargs[k])
5387        header = header+str(k)
5388        line = line+str(kwargs[k])
5389        count+=1
5390        if count <len(kwargs):
5391            header = header+','
5392            line = line+','
5393    header+='\n'
5394    line+='\n'
5395
5396    # checks the header info, if the same, then write, if not create a new file
5397    #try to open!
5398    try:
5399        fid = open(file,"r")
5400        file_header=fid.readline()
5401        fid.close()
5402        if verbose: print 'read file header %s' %file_header
5403       
5404    except:
5405        msg = 'try to create new file',file
5406        if verbose: print msg
5407        #tries to open file, maybe directory is bad
5408        try:
5409            fid = open(file,"w")
5410            fid.write(header)
5411            fid.close()
5412            file_header=header
5413        except:
5414            msg = 'cannot create new file',file
5415            raise msg
5416           
5417    #if header is same or this is a new file
5418    if file_header==str(header):
5419        fid=open(file,"a")
5420        #write new line
5421        fid.write(line)
5422        fid.close()
5423    else:
5424        #backup plan,
5425        # if header is different and has completed will append info to
5426        #end of details_temp.cvs file in output directory
5427        file = str(kwargs['output_dir'])+'detail_temp.csv'
5428        fid=open(file,"a")
5429        fid.write(header)
5430        fid.write(line)
5431        fid.close()
5432        if verbose: print 'file',file_header.strip('\n')
5433        if verbose: print 'head',header.strip('\n')
5434        if file_header.strip('\n')==str(header): print 'they equal'
5435        msg = 'WARNING: File header does not match input info, the input variables have changed, suggest to change file name'
5436        print msg
5437
5438
5439
5440
5441
5442def get_maximum_inundation_elevation(filename,
5443                                     polygon=None,
5444                                     time_interval=None,
5445                                     verbose=False):
5446   
5447    """Return highest elevation where depth > 0
5448   
5449    Usage:
5450    max_runup = get_maximum_inundation_elevation(filename,
5451                                                 polygon=None,
5452                                                 time_interval=None,
5453                                                 verbose=False)
5454
5455    filename is a NetCDF sww file containing ANUGA model output.   
5456    Optional arguments polygon and time_interval restricts the maximum
5457    runup calculation
5458    to a points that lie within the specified polygon and time interval.
5459
5460    If no inundation is found within polygon and time_interval the return value
5461    is None signifying "No Runup" or "Everything is dry".
5462
5463    See general function get_maximum_inundation_data for details.
5464   
5465    """
5466   
5467    runup, _ = get_maximum_inundation_data(filename,
5468                                           polygon=polygon,
5469                                           time_interval=time_interval,
5470                                           verbose=verbose)
5471    return runup
5472
5473
5474
5475
5476def get_maximum_inundation_location(filename,
5477                                    polygon=None,
5478                                    time_interval=None,
5479                                    verbose=False):
5480    """Return location of highest elevation where h > 0
5481   
5482   
5483    Usage:
5484    max_runup_location = get_maximum_inundation_location(filename,
5485                                                         polygon=None,
5486                                                         time_interval=None,
5487                                                         verbose=False)
5488
5489    filename is a NetCDF sww file containing ANUGA model output.
5490    Optional arguments polygon and time_interval restricts the maximum
5491    runup calculation
5492    to a points that lie within the specified polygon and time interval.
5493
5494    If no inundation is found within polygon and time_interval the return value
5495    is None signifying "No Runup" or "Everything is dry".
5496
5497    See general function get_maximum_inundation_data for details.
5498    """
5499   
5500    _, max_loc = get_maximum_inundation_data(filename,
5501                                             polygon=polygon,
5502                                             time_interval=time_interval,
5503                                             verbose=verbose)
5504    return max_loc
5505   
5506
5507
5508def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
5509                                use_centroid_values=False,
5510                                verbose=False):
5511    """Compute maximum run up height from sww file.
5512
5513
5514    Usage:
5515    runup, location = get_maximum_inundation_data(filename,
5516                                                  polygon=None,
5517                                                  time_interval=None,
5518                                                  verbose=False)
5519   
5520
5521    Algorithm is as in get_maximum_inundation_elevation from
5522    shallow_water_domain
5523    except that this function works with the sww file and computes the maximal
5524    runup height over multiple timesteps.
5525   
5526    Optional arguments polygon and time_interval restricts the
5527    maximum runup calculation
5528    to a points that lie within the specified polygon and time interval.
5529    Polygon is
5530    assumed to be in (absolute) UTM coordinates in the same zone as domain.
5531
5532    If no inundation is found within polygon and time_interval the return value
5533    is None signifying "No Runup" or "Everything is dry".
5534    """
5535
5536    # We are using nodal values here as that is what is stored in sww files.
5537
5538    # Water depth below which it is considered to be 0 in the model
5539    # FIXME (Ole): Allow this to be specified as a keyword argument as well
5540
5541    from anuga.utilities.polygon import inside_polygon   
5542    from anuga.config import minimum_allowed_height
5543    from Scientific.IO.NetCDF import NetCDFFile
5544
5545    dir, base = os.path.split(filename)
5546           
5547    iterate_over = get_all_swwfiles(dir,base)
5548   
5549    # Read sww file
5550    if verbose: 
5551        print 'Reading from %s' %filename
5552        # FIXME: Use general swwstats (when done)
5553   
5554    maximal_runup = None
5555    maximal_runup_location = None
5556   
5557    for file, swwfile in enumerate (iterate_over):
5558       
5559        # Read sww file
5560        filename = join(dir,swwfile+'.sww')
5561       
5562        if verbose: 
5563            print 'Reading from %s' %filename
5564            # FIXME: Use general swwstats (when done)
5565               
5566        fid = NetCDFFile(filename)
5567   
5568        # Get geo_reference
5569        # sww files don't have to have a geo_ref
5570        try:
5571            geo_reference = Geo_reference(NetCDFObject=fid)
5572        except AttributeError, e:
5573            geo_reference = Geo_reference() # Default georef object
5574           
5575        xllcorner = geo_reference.get_xllcorner()
5576        yllcorner = geo_reference.get_yllcorner()
5577        zone = geo_reference.get_zone()
5578       
5579        # Get extent
5580        volumes = fid.variables['volumes'][:]   
5581        x = fid.variables['x'][:] + xllcorner
5582        y = fid.variables['y'][:] + yllcorner
5583   
5584   
5585        # Get the relevant quantities
5586        elevation = fid.variables['elevation'][:] 
5587        stage = fid.variables['stage'][:]
5588   
5589   
5590        # Here's where one could convert nodal information to centroid
5591        # information
5592        # but is probably something we need to write in C.
5593        # Here's a Python thought which is NOT finished!!!
5594        if use_centroid_values is True:
5595            x = get_centroid_values(x, volumes)
5596            y = get_centroid_values(y, volumes)   
5597            elevation = get_centroid_values(elevation, volumes)   
5598   
5599   
5600        # Spatial restriction
5601        if polygon is not None:
5602            msg = 'polygon must be a sequence of points.'
5603            assert len(polygon[0]) == 2, msg
5604            # FIXME (Ole): Make a generic polygon input check in polygon.py
5605            # and call it here
5606           
5607            points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1)
5608   
5609            point_indices = inside_polygon(points, polygon)
5610   
5611            # Restrict quantities to polygon
5612            elevation = take(elevation, point_indices)
5613            stage = take(stage, point_indices, axis=1)
5614   
5615            # Get info for location of maximal runup
5616            points_in_polygon = take(points, point_indices)
5617            x = points_in_polygon[:,0]
5618            y = points_in_polygon[:,1]       
5619        else:
5620            # Take all points
5621            point_indices = arange(len(x))
5622           
5623   
5624        # Temporal restriction
5625        time = fid.variables['time'][:]
5626        all_timeindices = arange(len(time))       
5627        if time_interval is not None:
5628           
5629            msg = 'time_interval must be a sequence of length 2.'
5630            assert len(time_interval) == 2, msg
5631            msg = 'time_interval %s must not be decreasing.' %(time_interval)
5632            assert time_interval[1] >= time_interval[0], msg
5633           
5634            msg = 'Specified time interval [%.8f:%.8f]' %tuple(time_interval)
5635            msg += ' must does not match model time interval: [%.8f, %.8f]\n'\
5636                   %(time[0], time[-1])
5637            if time_interval[1] < time[0]: raise ValueError(msg)
5638            if time_interval[0] > time[-1]: raise ValueError(msg)
5639   
5640            # Take time indices corresponding to interval (& is bitwise AND)
5641            timesteps = compress((time_interval[0] <= time) & (time <= time_interval[1]),
5642                                 all_timeindices)
5643   
5644   
5645            msg = 'time_interval %s did not include any model timesteps.' %(time_interval)       
5646            assert not alltrue(timesteps == 0), msg
5647   
5648   
5649        else:
5650            # Take them all
5651            timesteps = all_timeindices
5652       
5653   
5654        fid.close()
5655   
5656        # Compute maximal runup for each timestep
5657        #maximal_runup = None
5658        #maximal_runup_location = None
5659        #maximal_runups = [None]
5660        #maximal_runup_locations = [None]
5661       
5662        for i in timesteps:
5663            if use_centroid_values is True:
5664                stage_i = get_centroid_values(stage[i,:], volumes)   
5665            else:
5666                stage_i = stage[i,:]
5667               
5668            depth = stage_i  - elevation
5669       
5670            # Get wet nodes i.e. nodes with depth>0 within given region and timesteps
5671            wet_nodes = compress(depth > minimum_allowed_height, arange(len(depth)))
5672   
5673            if alltrue(wet_nodes == 0):
5674                runup = None
5675            else:   
5676                # Find maximum elevation among wet nodes
5677                wet_elevation = take(elevation, wet_nodes)
5678   
5679                runup_index = argmax(wet_elevation)
5680                runup = max(wet_elevation)
5681                assert wet_elevation[runup_index] == runup # Must always be True
5682            if runup > maximal_runup:
5683                maximal_runup = runup      # This works even if maximal_runups is None
5684                #print "NEW RUNUP",runup
5685   
5686                # Record location
5687                wet_x = take(x, wet_nodes)
5688                wet_y = take(y, wet_nodes)           
5689                maximal_runup_location = [wet_x[runup_index], wet_y[runup_index]]
5690   
5691    #print 'maximal_runup, maximal_runup_location',maximal_runup, maximal_runup_location
5692    return maximal_runup, maximal_runup_location
5693
5694def get_all_swwfiles(look_in_dir='',base_name='',verbose=False):
5695    '''
5696    Finds all the sww files in a "look_in_dir" which contains a "base_name".
5697    will accept base_name with or without the extension ".sww"
5698   
5699    Returns: a list of strings
5700       
5701    Usage:     iterate_over = get_all_swwfiles(dir, name)
5702    then
5703               for swwfile in iterate_over:
5704                   do stuff
5705                   
5706    Check "export_grids" and "get_maximum_inundation_data" for examples
5707    '''
5708   
5709    #plus tests the extension
5710    name, extension = os.path.splitext(base_name)
5711
5712    if extension <>'' and extension <> '.sww':
5713        msg = msg='file %s %s must be an NetCDF sww file!'%(base_name,extension)
5714        raise IOError, msg
5715
5716    if look_in_dir == "":
5717        look_in_dir = "." # Unix compatibility
5718   
5719    dir_ls = os.listdir(look_in_dir)
5720    #print 'dir_ls',dir_ls, base
5721    iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == '.sww']
5722    if len(iterate_over) == 0:
5723        msg = 'No files of the base name %s'\
5724              %(name)
5725        raise IOError, msg
5726    if verbose: print 'iterate over %s' %(iterate_over)
5727
5728    return iterate_over
5729
5730def get_all_files_with_extension(look_in_dir='',base_name='',extension='.sww',verbose=False):
5731    '''
5732    Finds all the sww files in a "look_in_dir" which contains a "base_name".
5733   
5734   
5735    Returns: a list of strings
5736       
5737    Usage:     iterate_over = get_all_swwfiles(dir, name)
5738    then
5739               for swwfile in iterate_over:
5740                   do stuff
5741                   
5742    Check "export_grids" and "get_maximum_inundation_data" for examples
5743    '''
5744   
5745    #plus tests the extension
5746    name, ext = os.path.splitext(base_name)
5747#    print 'look_in_dir',look_in_dir
5748
5749    if ext <>'' and ext <> extension:
5750        msg = msg='base_name %s must be an file with %s extension!'%(base_name,extension)
5751        raise IOError, msg
5752
5753    if look_in_dir == "":
5754        look_in_dir = "." # Unix compatibility
5755#    print 'look_in_dir',look_in_dir, getcwd()
5756    dir_ls = os.listdir(look_in_dir)
5757    #print 'dir_ls',dir_ls, base_name
5758    iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == extension]
5759    if len(iterate_over) == 0:
5760        msg = 'No files of the base name %s in %s'\
5761              %(name, look_in_dir)
5762        raise IOError, msg
5763    if verbose: print 'iterate over %s' %(iterate_over)
5764
5765    return iterate_over
5766
5767def get_all_directories_with_name(look_in_dir='',base_name='',verbose=False):
5768    '''
5769    Finds all the sww files in a "look_in_dir" which contains a "base_name".
5770   
5771   
5772    Returns: a list of strings
5773       
5774    Usage:     iterate_over = get_all_swwfiles(dir, name)
5775    then
5776               for swwfile in iterate_over:
5777                   do stuff
5778                   
5779    Check "export_grids" and "get_maximum_inundation_data" for examples
5780    '''
5781   
5782    #plus tests the extension
5783
5784    if look_in_dir == "":
5785        look_in_dir = "." # Unix compatibility
5786#    print 'look_in_dir',look_in_dir
5787    dir_ls = os.listdir(look_in_dir)
5788#    print 'dir_ls',dir_ls
5789    iterate_over = [x for x in dir_ls if base_name in x]
5790    if len(iterate_over) == 0:
5791        msg = 'No files of the base name %s'\
5792              %(name)
5793        raise IOError, msg
5794    if verbose: print 'iterate over %s' %(iterate_over)
5795
5796    return iterate_over
5797
5798
5799#-------------------------------------------------------------
5800if __name__ == "__main__":
5801    #setting umask from config to force permissions for all files and directories
5802    # created to the same. (it was noticed the "mpirun" doesn't honour the umask
5803    # set in your .bashrc etc file)
5804    from config import umask
5805    import os 
5806    os.umask(umask)
5807
Note: See TracBrowser for help on using the repository browser.