source: anuga_core/source/anuga/shallow_water/data_manager_joaquims_patch.py @ 7340

Last change on this file since 7340 was 7340, checked in by ole, 15 years ago

More refactoring in preparation for ticket:191
These are mostly simplifications of sww file creation

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