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

Last change on this file since 5470 was 5470, checked in by ole, 14 years ago

Added verbose flag to c code

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