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

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

Allow urs_ext to store only a subset of gauges in mux file according to a list of ordered points

File size: 221.9 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,permutation=None,verbose=False):
4841    """
4842    Access the mux files specified in the filenames list. Combine the
4843    data found therin as a weighted linear sum as specifed by the weights.
4844    If permutation is None extract timeseries data for all gauges within the
4845    files.
4846   
4847    Input:
4848    filenames:   List of filenames specifiying the file containing the
4849                 timeseries data (mux2 format) for each source
4850    weights:     Weights associated with each source
4851    permutation: Specifies the gauge numbers that for which data is to be
4852                 extracted
4853    """
4854    from Numeric import ones,Float,compress,zeros,arange
4855    from urs_ext import read_mux2
4856
4857    numSrc=len(filenames)
4858
4859    file_params=-1*ones(3,Float)#[nsta,dt,nt]
4860   
4861    # Convert verbose to int C flag
4862    if verbose:
4863        verbose=1
4864    else:
4865        verbose=0
4866       
4867    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
4868
4869    msg='File parameter values were not read in correctly from c file'
4870    assert len(compress(file_params>0,file_params))!=0,msg
4871    msg='The number of stations specifed in the c array and in the file are inconsitent'
4872    assert file_params[0]==data.shape[0],msg
4873   
4874    nsta=int(file_params[0])
4875    msg='Must have at least one station'
4876    assert nsta>0,msg
4877    dt=file_params[1]
4878    msg='Must have a postive timestep'
4879    assert dt>0,msg
4880    nt=int(file_params[2])
4881    msg='Must have at least one gauge value'
4882    assert nt>0,msg
4883   
4884    OFFSET=5 #number of site parameters p passed back with data
4885    #p=[geolat,geolon,depth,start_tstep,finish_tstep]
4886
4887    times=dt*arange(1,(data.shape[1]-OFFSET)+1)
4888    latitudes=zeros(data.shape[0],Float)
4889    longitudes=zeros(data.shape[0],Float)
4890    elevation=zeros(data.shape[0],Float)
4891    quantity=zeros((data.shape[0],data.shape[1]-OFFSET),Float)
4892    starttime=1e16
4893    for i in range(0,data.shape[0]):
4894        latitudes[i]=data[i][data.shape[1]-OFFSET]
4895        longitudes[i]=data[i][data.shape[1]-OFFSET+1]
4896        elevation[i]=-data[i][data.shape[1]-OFFSET+2]
4897        quantity[i]=data[i][:-OFFSET]
4898        starttime=min(dt*data[i][data.shape[1]-OFFSET+3],starttime)
4899       
4900    return times, latitudes, longitudes, elevation, quantity, starttime
4901
4902def mux2sww_time(mux_times, mint, maxt):
4903    """
4904    """
4905
4906    if mint == None:
4907        mux_times_start_i = 0
4908    else:
4909        mux_times_start_i = searchsorted(mux_times, mint)
4910       
4911    if maxt == None:
4912        mux_times_fin_i = len(mux_times)
4913    else:
4914        maxt += 0.5 # so if you specify a time where there is
4915                    # data that time will be included
4916        mux_times_fin_i = searchsorted(mux_times, maxt)
4917
4918    return mux_times_start_i, mux_times_fin_i
4919
4920
4921def urs2sts(basename_in, basename_out = None, weights=None,
4922            verbose = False, origin = None,
4923            mean_stage=0.0, zscale=1.0,
4924            ordering_filename = None,
4925            minlat = None, maxlat = None,
4926            minlon = None, maxlon = None):
4927    """Convert URS mux2 format for wave propagation to sts format
4928
4929    Also convert latitude and longitude to UTM. All coordinates are
4930    assumed to be given in the GDA94 datum
4931
4932    origin is a 3-tuple with geo referenced
4933    UTM coordinates (zone, easting, northing)
4934   
4935    input:
4936    basename_in: list of source file prefixes
4937   
4938    These are combined with the extensions:
4939    WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage
4940    EAST_VELOCITY_MUX2_LABEL =  '-e-mux2' xmomentum
4941    NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2' and ymomentum
4942   
4943    to create a 2D list of mux2 file. The rows are associated with each quantity and must have the above extensions
4944    the columns are the list of file prefixes.
4945   
4946    ordering: a .txt file name specifying which mux2 gauge points are sto be stored. This is indicated by the index of
4947    the gauge in the ordering file.
4948   
4949    ordering file format:
4950    header
4951    index,longitude,latitude
4952   
4953    header='index,longitude,latitude\n'
4954    If ordering is None all points are taken in the order they appear in the mux2 file subject to rectangular clipping box (see below).
4955   
4956   
4957    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
4958    information at gauge points not on the boundary.
4959    if ordering_filename is specified clipping box cannot be used
4960   
4961   
4962    output:
4963    baename_out: name of sts file in which mux2 data is stored.
4964   
4965   
4966    """
4967    import os
4968    from Scientific.IO.NetCDF import NetCDFFile
4969    from Numeric import Float, Int, Int32, searchsorted, zeros, array
4970    from Numeric import allclose, around,ones,Float
4971    from types import ListType,StringType
4972    from operator import __and__
4973   
4974    if not isinstance(basename_in, ListType):
4975        if verbose: print 'Reading single source'
4976        basename_in=[basename_in]
4977
4978    # Check that basename is a list of strings
4979    if not reduce(__and__, map(lambda z:isinstance(z,StringType),basename_in)):
4980        msg= 'basename_in must be a string or list of strings'
4981        raise Exception, msg
4982
4983    # Find the number of sources to be used
4984    numSrc=len(basename_in)
4985
4986    # A weight must be specified for each source
4987    if weights is None:
4988        weights=ones(numSrc,Float)/numSrc
4989    else:
4990        weights = ensure_numeric(weights)
4991        msg = 'When combining multiple sources must specify a weight for '+\
4992              'mux2 source file'
4993        assert len(weights)== numSrc, msg
4994
4995    precision = Float
4996
4997    if ordering_filename is not None:
4998        msg = 'If ordering_filename is specified,'
4999        msg += ' rectangular clipping box cannot be specified as well. '
5000        assert minlat is None and maxlat is None and\
5001            minlon is None and maxlon is None, msg
5002   
5003   
5004    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
5005    if minlat != None:
5006        assert -90 < minlat < 90 , msg
5007    if maxlat != None:
5008        assert -90 < maxlat < 90 , msg
5009        if minlat != None:
5010            assert maxlat > minlat
5011    if minlon != None:
5012        assert -180 < minlon < 180 , msg
5013    if maxlon != None:
5014        assert -180 < maxlon < 180 , msg
5015        if minlon != None:
5016            assert maxlon > minlon
5017
5018
5019    # Check output filename   
5020    if basename_out is None:
5021        msg = 'STS filename must be specified'
5022        raise Exception, msg
5023   
5024    stsname = basename_out + '.sts'
5025
5026    files_in=[[],[],[]]
5027    for files in basename_in:
5028        files_in[0].append(files + WAVEHEIGHT_MUX2_LABEL),
5029        files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL)
5030        files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL)
5031       
5032    #files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
5033    #            basename_in + EAST_VELOCITY_MUX2_LABEL,
5034    #            basename_in + NORTH_VELOCITY_MUX2_LABEL]
5035   
5036    quantities = ['HA','UA','VA']
5037
5038    # For each source file check that there exists three files ending with:
5039    # WAVEHEIGHT_MUX2_LABEL,
5040    # EAST_VELOCITY_MUX2_LABEL, and
5041    # NORTH_VELOCITY_MUX2_LABEL
5042    for i in range(len(quantities)): 
5043        for file_in in files_in[i]:
5044            if (os.access(file_in, os.F_OK) == 0):
5045                msg = 'File %s does not exist or is not accessible' %file_in
5046                raise IOError, msg
5047
5048    if ordering_filename is not None:
5049        # Read ordering file
5050       
5051        try:
5052            fid=open(ordering_filename, 'r')
5053            file_header=fid.readline().split(',')
5054            ordering_lines=fid.readlines()
5055            fid.close()
5056        except:
5057            msg = 'Cannot open %s'%ordering_filename
5058            raise Exception, msg
5059
5060        reference_header = 'index, longitude, latitude\n'
5061        reference_header_split = reference_header.split(',')
5062        for i in range(3):
5063            if not file_header[i] == reference_header_split[i]:
5064                msg = 'File must contain header: '+reference_header+'\n'
5065                raise Exception, msg
5066
5067        if len(ordering_lines)<2:
5068            msg = 'File must contain at least two points'
5069            raise Exception, msg
5070
5071   
5072       
5073        permutation = [int(line.split(',')[0]) for line in ordering_lines]
5074    else:
5075        permutation = None 
5076
5077    if (verbose): print 'reading mux2 file'
5078    mux={}
5079    for i,quantity in enumerate(quantities):
5080        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
5081        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
5082        if quantity!=quantities[0]:
5083            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
5084            assert allclose(times,times_old),msg
5085            assert allclose(latitudes_urs,latitudes_urs_old),msg
5086            assert allclose(longitudes_urs,longitudes_urs_old),msg
5087            assert allclose(elevation,elevation_old),msg
5088            assert allclose(starttime,starttime_old)
5089        times_old=times
5090        latitudes_urs_old=latitudes_urs
5091        longitudes_urs_old=longitudes_urs
5092        elevation_old=elevation
5093        starttime_old=starttime
5094       
5095    if ordering_filename is not None:
5096        latitudes = latitudes_urs
5097        longitudes = longitudes_urs
5098       
5099        # Self check - can be removed to improve speed
5100        ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]               
5101        ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]       
5102       
5103        msg = 'Longitudes specified in ordering file do not match those found in mux files'
5104        msg += 'I got %s instead of %s (only beginning shown)'\
5105            %(str(longitudes[:10]) + '...',
5106              str(ref_longitudes[:10]) + '...')       
5107        assert allclose(longitudes, ref_longitudes), msg       
5108       
5109        msg = 'Latitudes specified in ordering file do not match those found in mux files. '
5110        msg += 'I got %s instead of %s (only beginning shown)'\
5111            %(str(latitudes[:10]) + '...',
5112              str(ref_latitudes[:10]) + '...')               
5113        assert allclose(latitudes, ref_latitudes), msg
5114       
5115    elif (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
5116        #FIXME - this branch is probably not working   
5117        if verbose: print 'Cliping urs data'
5118        latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)
5119        longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)
5120       
5121        permutation = range(latitudes.shape[0])       
5122       
5123        msg = 'Clipping is not done yet'
5124        raise Exception, msg
5125       
5126    else:
5127        latitudes=latitudes_urs
5128        longitudes=longitudes_urs
5129        permutation = range(latitudes.shape[0])               
5130
5131    msg='File is empty and or clipped region not in file region'
5132    assert len(latitudes>0),msg
5133
5134    number_of_points = latitudes.shape[0]
5135    number_of_times = times.shape[0]
5136    number_of_latitudes = latitudes.shape[0]
5137    number_of_longitudes = longitudes.shape[0]
5138
5139    # NetCDF file definition
5140    outfile = NetCDFFile(stsname, 'w')
5141
5142    description = 'Converted from URS mux2 files: %s'\
5143                  %(basename_in)
5144   
5145    # Create new file
5146    #starttime = times[0]
5147    sts = Write_sts()
5148    sts.store_header(outfile, 
5149                     times+starttime,
5150                     number_of_points, 
5151                     description=description,
5152                     verbose=verbose,
5153                     sts_precision=Float)
5154   
5155    # Store
5156    from anuga.coordinate_transforms.redfearn import redfearn
5157    x = zeros(number_of_points, Float)  #Easting
5158    y = zeros(number_of_points, Float)  #Northing
5159
5160    # Check zone boundaries
5161    refzone, _, _ = redfearn(latitudes[0],longitudes[0]) 
5162
5163    for i in range(number_of_points):
5164        zone, easting, northing = redfearn(latitudes[i],longitudes[i])
5165        x[i] = easting
5166        y[i] = northing
5167        msg='all sts gauges need to be in the same zone'
5168        assert zone==refzone,msg
5169
5170    if origin is None:
5171        origin = Geo_reference(refzone,min(x),min(y))
5172    geo_ref = write_NetCDF_georeference(origin, outfile)
5173
5174    z = elevation
5175   
5176    #print geo_ref.get_xllcorner()
5177    #print geo_ref.get_yllcorner()
5178
5179    z = resize(z,outfile.variables['z'][:].shape)
5180    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
5181    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
5182    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
5183    outfile.variables['elevation'][:] = z
5184
5185    stage = outfile.variables['stage']
5186    xmomentum = outfile.variables['xmomentum']
5187    ymomentum = outfile.variables['ymomentum']
5188
5189    if verbose: print 'Converting quantities'
5190    for j in range(len(times)):
5191        for i in range(number_of_points):
5192            w = zscale*mux['HA'][i,j] + mean_stage
5193            h=w-elevation[i]
5194            stage[j,i] = w
5195
5196            xmomentum[j,i] = mux['UA'][i,j]*h
5197            ymomentum[j,i] = mux['VA'][i,j]*h
5198
5199    outfile.close()
5200
5201def create_sts_boundary(stsname):
5202    """
5203        Create boundary segments from .sts file. Points can be stored in
5204    arbitrary order within the .sts file. The order in which the .sts points
5205    make up the boundary are given in order.txt file
5206   
5207    FIXME:
5208    Point coordinates are stored in relative eastings and northings.
5209    But boundary is produced in absolute coordinates
5210   
5211    """
5212    try:
5213        fid = NetCDFFile(stsname+'.sts', 'r')
5214    except:
5215        msg = 'Cannot open %s'%filename+'.sts'
5216
5217
5218    xllcorner = fid.xllcorner[0]
5219    yllcorner = fid.yllcorner[0]
5220    #Points stored in sts file are normalised to [xllcorner,yllcorner] but
5221    #we cannot assume that boundary polygon will be. At least the
5222    #additional points specified by the user after this function is called
5223    x = fid.variables['x'][:]+xllcorner
5224    y = fid.variables['y'][:]+yllcorner
5225
5226    x = reshape(x, (len(x),1))
5227    y = reshape(y, (len(y),1))
5228    sts_points = concatenate((x,y), axis=1)
5229
5230    return sts_points.tolist()
5231
5232class Write_sww:
5233    from anuga.shallow_water.shallow_water_domain import Domain
5234
5235    # FIXME (Ole): Hardwiring the conserved quantities like
5236    # this could be a problem. I would prefer taking them from
5237    # the instantiation of Domain.
5238    #
5239    # (DSG) There is not always a Domain instance when Write_sww is used.
5240    # Check to see if this is the same level of hardwiring as is in
5241    # shallow water doamain.
5242   
5243    sww_quantities = Domain.conserved_quantities
5244
5245
5246    RANGE = '_range'
5247    EXTREMA = ':extrema'
5248
5249    def __init__(self):
5250        pass
5251   
5252    def store_header(self,
5253                     outfile,
5254                     times,
5255                     number_of_volumes,
5256                     number_of_points,
5257                     description='Converted from XXX',
5258                     smoothing=True,
5259                     order=1,
5260                     sww_precision=Float32,
5261                     verbose=False):
5262        """
5263        outfile - the name of the file that will be written
5264        times - A list of the time slice times OR a start time
5265        Note, if a list is given the info will be made relative.
5266        number_of_volumes - the number of triangles
5267        """
5268   
5269        outfile.institution = 'Geoscience Australia'
5270        outfile.description = description
5271
5272        # For sww compatibility
5273        if smoothing is True:
5274            # Smoothing to be depreciated
5275            outfile.smoothing = 'Yes'
5276            outfile.vertices_are_stored_uniquely = 'False'
5277        else:
5278            # Smoothing to be depreciated
5279            outfile.smoothing = 'No'
5280            outfile.vertices_are_stored_uniquely = 'True'
5281        outfile.order = order
5282
5283        try:
5284            revision_number = get_revision_number()
5285        except:
5286            revision_number = None
5287        # Allow None to be stored as a string               
5288        outfile.revision_number = str(revision_number) 
5289
5290
5291       
5292        # times - A list or array of the time slice times OR a start time
5293        # times = ensure_numeric(times)
5294        # Start time in seconds since the epoch (midnight 1/1/1970)
5295
5296        # This is being used to seperate one number from a list.
5297        # what it is actually doing is sorting lists from numeric arrays.
5298        if type(times) is list or type(times) is ArrayType: 
5299            number_of_times = len(times)
5300            times = ensure_numeric(times) 
5301            if number_of_times == 0:
5302                starttime = 0
5303            else:
5304                starttime = times[0]
5305                times = times - starttime  #Store relative times
5306        else:
5307            number_of_times = 0
5308            starttime = times
5309            #times = ensure_numeric([])
5310        outfile.starttime = starttime
5311        # dimension definitions
5312        outfile.createDimension('number_of_volumes', number_of_volumes)
5313        outfile.createDimension('number_of_vertices', 3)
5314        outfile.createDimension('numbers_in_range', 2)
5315   
5316        if smoothing is True:
5317            outfile.createDimension('number_of_points', number_of_points)
5318       
5319            # FIXME(Ole): This will cause sww files for paralle domains to
5320            # have ghost nodes stored (but not used by triangles).
5321            # To clean this up, we have to change get_vertex_values and
5322            # friends in quantity.py (but I can't be bothered right now)
5323        else:
5324            outfile.createDimension('number_of_points', 3*number_of_volumes)
5325        outfile.createDimension('number_of_timesteps', number_of_times)
5326
5327        # variable definitions
5328        outfile.createVariable('x', sww_precision, ('number_of_points',))
5329        outfile.createVariable('y', sww_precision, ('number_of_points',))
5330        outfile.createVariable('elevation', sww_precision, ('number_of_points',))
5331        q = 'elevation'
5332        outfile.createVariable(q+Write_sww.RANGE, sww_precision,
5333                               ('numbers_in_range',))
5334
5335
5336        # Initialise ranges with small and large sentinels.
5337        # If this was in pure Python we could have used None sensibly
5338        outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
5339        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
5340
5341        # FIXME: Backwards compatibility
5342        outfile.createVariable('z', sww_precision, ('number_of_points',))
5343        #################################
5344
5345        outfile.createVariable('volumes', Int, ('number_of_volumes',
5346                                                'number_of_vertices'))
5347        # Doing sww_precision instead of Float gives cast errors.
5348        outfile.createVariable('time', Float,
5349                               ('number_of_timesteps',))
5350       
5351        for q in Write_sww.sww_quantities:
5352            outfile.createVariable(q, sww_precision,
5353                                   ('number_of_timesteps',
5354                                    'number_of_points')) 
5355            outfile.createVariable(q+Write_sww.RANGE, sww_precision,
5356                                   ('numbers_in_range',))
5357
5358            # Initialise ranges with small and large sentinels.
5359            # If this was in pure Python we could have used None sensibly
5360            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
5361            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
5362           
5363        if type(times) is list or type(times) is ArrayType: 
5364            outfile.variables['time'][:] = times    #Store time relative
5365           
5366        if verbose:
5367            print '------------------------------------------------'
5368            print 'Statistics:'
5369            print '    t in [%f, %f], len(t) == %d'\
5370                  %(min(times.flat), max(times.flat), len(times.flat))
5371
5372       
5373    def store_triangulation(self,
5374                            outfile,
5375                            points_utm,
5376                            volumes,
5377                            elevation, zone=None, new_origin=None, 
5378                            points_georeference=None, verbose=False):
5379        """
5380       
5381        new_origin - qa georeference that the points can be set to. (Maybe
5382        do this before calling this function.)
5383       
5384        points_utm - currently a list or array of the points in UTM.
5385        points_georeference - the georeference of the points_utm
5386       
5387        How about passing new_origin and current_origin.
5388        If you get both, do a convertion from the old to the new.
5389       
5390        If you only get new_origin, the points are absolute,
5391        convert to relative
5392       
5393        if you only get the current_origin the points are relative, store
5394        as relative.
5395       
5396        if you get no georefs create a new georef based on the minimums of
5397        points_utm.  (Another option would be to default to absolute)
5398       
5399        Yes, and this is done in another part of the code.
5400        Probably geospatial.
5401       
5402        If you don't supply either geo_refs, then supply a zone. If not
5403        the default zone will be used.
5404       
5405       
5406        precon
5407       
5408        header has been called.
5409        """
5410       
5411        number_of_points = len(points_utm)   
5412        volumes = array(volumes) 
5413        points_utm = array(points_utm)
5414
5415        # given the two geo_refs and the points, do the stuff
5416        # described in the method header
5417        # if this is needed else where, pull out as a function
5418        points_georeference = ensure_geo_reference(points_georeference)
5419        new_origin = ensure_geo_reference(new_origin)
5420        if new_origin is None and points_georeference is not None:
5421            points = points_utm
5422            geo_ref = points_georeference
5423        else:
5424            if new_origin is None:
5425                new_origin = Geo_reference(zone,min(points_utm[:,0]),
5426                                           min(points_utm[:,1]))
5427            points = new_origin.change_points_geo_ref(points_utm,
5428                                                      points_georeference)
5429            geo_ref = new_origin
5430
5431        # At this stage I need a georef and points
5432        # the points are relative to the georef
5433        geo_ref.write_NetCDF(outfile)
5434   
5435        # This will put the geo ref in the middle
5436        #geo_ref=Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
5437       
5438        x =  points[:,0]
5439        y =  points[:,1]
5440        z = outfile.variables['z'][:]
5441   
5442        if verbose:
5443            print '------------------------------------------------'
5444            print 'More Statistics:'
5445            print '  Extent (/lon):'
5446            print '    x in [%f, %f], len(lat) == %d'\
5447                  %(min(x), max(x),
5448                    len(x))
5449            print '    y in [%f, %f], len(lon) == %d'\
5450                  %(min(y), max(y),
5451                    len(y))
5452            print '    z in [%f, %f], len(z) == %d'\
5453                  %(min(elevation), max(elevation),
5454                    len(elevation))
5455            print 'geo_ref: ',geo_ref
5456            print '------------------------------------------------'
5457           
5458        #z = resize(bath_grid,outfile.variables['z'][:].shape)
5459        #print "points[:,0]", points[:,0]
5460        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
5461        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
5462        outfile.variables['z'][:] = elevation
5463        outfile.variables['elevation'][:] = elevation  #FIXME HACK
5464        outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
5465
5466        q = 'elevation'
5467        # This updates the _range values
5468        outfile.variables[q+Write_sww.RANGE][0] = min(elevation)
5469        outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
5470
5471
5472    def store_quantities(self, outfile, sww_precision=Float32,
5473                         slice_index=None, time=None,
5474                         verbose=False, **quant):
5475        """
5476        Write the quantity info.
5477
5478        **quant is extra keyword arguments passed in. These must be
5479          the sww quantities, currently; stage, xmomentum, ymomentum.
5480       
5481        if the time array is already been built, use the slice_index
5482        to specify the index.
5483       
5484        Otherwise, use time to increase the time dimension
5485
5486        Maybe make this general, but the viewer assumes these quantities,
5487        so maybe we don't want it general - unless the viewer is general
5488       
5489        precon
5490        triangulation and
5491        header have been called.
5492        """
5493
5494        if time is not None:
5495            file_time = outfile.variables['time']
5496            slice_index = len(file_time)
5497            file_time[slice_index] = time   
5498
5499        # Write the conserved quantities from Domain.
5500        # Typically stage,  xmomentum, ymomentum
5501        # other quantities will be ignored, silently.
5502        # Also write the ranges: stage_range,
5503        # xmomentum_range and ymomentum_range
5504        for q in Write_sww.sww_quantities:
5505            if not quant.has_key(q):
5506                msg = 'SWW file can not write quantity %s' %q
5507                raise NewQuantity, msg
5508            else:
5509                q_values = quant[q]
5510                outfile.variables[q][slice_index] = \
5511                                q_values.astype(sww_precision)
5512
5513                # This updates the _range values
5514                q_range = outfile.variables[q+Write_sww.RANGE][:]
5515                q_values_min = min(q_values)
5516                if q_values_min < q_range[0]:
5517                    outfile.variables[q+Write_sww.RANGE][0] = q_values_min
5518                q_values_max = max(q_values)
5519                if q_values_max > q_range[1]:
5520                    outfile.variables[q+Write_sww.RANGE][1] = q_values_max
5521
5522    def verbose_quantities(self, outfile):
5523        print '------------------------------------------------'
5524        print 'More Statistics:'
5525        for q in Write_sww.sww_quantities:
5526            print %s in [%f, %f]' %(q,
5527                                       outfile.variables[q+Write_sww.RANGE][0],
5528                                       outfile.variables[q+Write_sww.RANGE][1])
5529        print '------------------------------------------------'
5530
5531
5532       
5533def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation,
5534                         mean_stage=0, zscale=1,
5535                         verbose=False):   
5536    #Time stepping
5537    stage = outfile.variables['stage']
5538    xmomentum = outfile.variables['xmomentum']
5539    ymomentum = outfile.variables['ymomentum']
5540
5541    n = len(has)
5542    j=0
5543    for ha, ua, va in map(None, has, uas, vas):
5544        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
5545        w = zscale*ha + mean_stage
5546        stage[j] = w
5547        h = w - elevation
5548        xmomentum[j] = ua*h
5549        ymomentum[j] = -1*va*#  -1 since in mux files south is positive.
5550        j += 1
5551   
5552def urs2txt(basename_in, location_index=None):
5553    """
5554    Not finished or tested
5555    """
5556   
5557    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
5558                basename_in + EAST_VELOCITY_LABEL,
5559                basename_in + NORTH_VELOCITY_LABEL]
5560    quantities = ['HA','UA','VA']
5561
5562    d = ","
5563   
5564    # instanciate urs_points of the three mux files.
5565    mux = {}
5566    for quantity, file in map(None, quantities, files_in):
5567        mux[quantity] = Urs_points(file)
5568       
5569    # Could check that the depth is the same. (hashing)
5570
5571    # handle to a mux file to do depth stuff
5572    a_mux = mux[quantities[0]]
5573   
5574    # Convert to utm
5575    latitudes = a_mux.lonlatdep[:,1]
5576    longitudes = a_mux.lonlatdep[:,0]
5577    points_utm, zone = convert_from_latlon_to_utm( \
5578        latitudes=latitudes, longitudes=longitudes)
5579    #print "points_utm", points_utm
5580    #print "zone", zone
5581    depths = a_mux.lonlatdep[:,2]  #
5582   
5583    fid = open(basename_in + '.txt', 'w')
5584
5585    fid.write("zone: " + str(zone) + "\n")
5586
5587    if location_index is not None:
5588        #Title
5589        li = location_index
5590        fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5591                  'Northing' + "\n")
5592        fid.write(str(li) +d+ str(latitudes[li])+d+ \
5593              str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \
5594              str(points_utm[li][01]) + "\n")
5595
5596    # the non-time dependent stuff
5597    #Title
5598    fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
5599              'Northing' +d+ 'depth m' + "\n")
5600    i = 0
5601    for depth, point_utm, lat, long in map(None, depths,
5602                                               points_utm, latitudes,
5603                                               longitudes):
5604       
5605        fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \
5606                  str(point_utm[01]) +d+ str(depth) + "\n")
5607        i +=1
5608    #Time dependent
5609    if location_index is not None:
5610        time_step = a_mux.time_step
5611        i = 0
5612        #Title
5613        fid.write('time' +d+ 'HA depth m'+d+ \
5614                 'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \
5615                      + "\n")
5616        for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
5617            fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \
5618                      str(UA[location_index]) +d+ str(VA[location_index]) \
5619                      + "\n")
5620           
5621            i +=1
5622
5623class Write_sts:
5624
5625    sts_quantities = ['stage','xmomentum','ymomentum']
5626
5627
5628    RANGE = '_range'
5629    EXTREMA = ':extrema'
5630   
5631    def __init__(self):
5632        pass
5633
5634    def store_header(self,
5635                     outfile,
5636                     times,
5637                     number_of_points,
5638                     description='Converted from URS mux2 format',
5639                     sts_precision=Float32,
5640                     verbose=False):
5641        """
5642        outfile - the name of the file that will be written
5643        times - A list of the time slice times OR a start time
5644        Note, if a list is given the info will be made relative.
5645        number_of_points - the number of urs gauges sites
5646        """
5647
5648        outfile.institution = 'Geoscience Australia'
5649        outfile.description = description
5650       
5651        try:
5652            revision_number = get_revision_number()
5653        except:
5654            revision_number = None
5655        # Allow None to be stored as a string               
5656        outfile.revision_number = str(revision_number) 
5657       
5658        # times - A list or array of the time slice times OR a start time
5659        # Start time in seconds since the epoch (midnight 1/1/1970)
5660
5661        # This is being used to seperate one number from a list.
5662        # what it is actually doing is sorting lists from numeric arrays.
5663        if type(times) is list or type(times) is ArrayType: 
5664            number_of_times = len(times)
5665            times = ensure_numeric(times) 
5666            if number_of_times == 0:
5667                starttime = 0
5668            else:
5669                starttime = times[0]
5670                times = times - starttime  #Store relative times
5671        else:
5672            number_of_times = 0
5673            starttime = times
5674
5675        outfile.starttime = starttime
5676
5677        # Dimension definitions
5678        outfile.createDimension('number_of_points', number_of_points)
5679        outfile.createDimension('number_of_timesteps', number_of_times)
5680        outfile.createDimension('numbers_in_range', 2)
5681
5682        # Variable definitions
5683        outfile.createVariable('x', sts_precision, ('number_of_points',))
5684        outfile.createVariable('y', sts_precision, ('number_of_points',))
5685        outfile.createVariable('elevation', sts_precision,('number_of_points',))
5686
5687        q = 'elevation'
5688        outfile.createVariable(q+Write_sts.RANGE, sts_precision,
5689                               ('numbers_in_range',))
5690
5691        # Initialise ranges with small and large sentinels.
5692        # If this was in pure Python we could have used None sensibly
5693        outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
5694        outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
5695
5696        outfile.createVariable('z', sts_precision, ('number_of_points',))
5697        # Doing sts_precision instead of Float gives cast errors.
5698        outfile.createVariable('time', Float, ('number_of_timesteps',))
5699
5700        for q in Write_sts.sts_quantities:
5701            outfile.createVariable(q, sts_precision,
5702                                   ('number_of_timesteps',
5703                                    'number_of_points'))
5704            outfile.createVariable(q+Write_sts.RANGE, sts_precision,
5705                                   ('numbers_in_range',))
5706            # Initialise ranges with small and large sentinels.
5707            # If this was in pure Python we could have used None sensibly
5708            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
5709            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
5710
5711        if type(times) is list or type(times) is ArrayType: 
5712            outfile.variables['time'][:] = times    #Store time relative
5713
5714        if verbose:
5715            print '------------------------------------------------'
5716            print 'Statistics:'
5717            print '    t in [%f, %f], len(t) == %d'\
5718                  %(min(times.flat), max(times.flat), len(times.flat))
5719
5720    def store_points(self,
5721                     outfile,
5722                     points_utm,
5723                     elevation, zone=None, new_origin=None, 
5724                     points_georeference=None, verbose=False):
5725
5726        """
5727        points_utm - currently a list or array of the points in UTM.
5728        points_georeference - the georeference of the points_utm
5729
5730        How about passing new_origin and current_origin.
5731        If you get both, do a convertion from the old to the new.
5732       
5733        If you only get new_origin, the points are absolute,
5734        convert to relative
5735       
5736        if you only get the current_origin the points are relative, store
5737        as relative.
5738       
5739        if you get no georefs create a new georef based on the minimums of
5740        points_utm.  (Another option would be to default to absolute)
5741       
5742        Yes, and this is done in another part of the code.
5743        Probably geospatial.
5744       
5745        If you don't supply either geo_refs, then supply a zone. If not
5746        the default zone will be used.
5747
5748        precondition:
5749             header has been called.
5750        """
5751
5752        number_of_points = len(points_utm)
5753        points_utm = array(points_utm)
5754
5755        # given the two geo_refs and the points, do the stuff
5756        # described in the method header
5757        points_georeference = ensure_geo_reference(points_georeference)
5758        new_origin = ensure_geo_reference(new_origin)
5759       
5760        if new_origin is None and points_georeference is not None:
5761            points = points_utm
5762            geo_ref = points_georeference
5763        else:
5764            if new_origin is None:
5765                new_origin = Geo_reference(zone,min(points_utm[:,0]),
5766                                           min(points_utm[:,1]))
5767            points = new_origin.change_points_geo_ref(points_utm,
5768                                                      points_georeference)
5769            geo_ref = new_origin
5770
5771        # At this stage I need a georef and points
5772        # the points are relative to the georef
5773        geo_ref.write_NetCDF(outfile)
5774
5775        x =  points[:,0]
5776        y =  points[:,1]
5777        z = outfile.variables['z'][:]
5778   
5779        if verbose:
5780            print '------------------------------------------------'
5781            print 'More Statistics:'
5782            print '  Extent (/lon):'
5783            print '    x in [%f, %f], len(lat) == %d'\
5784                  %(min(x), max(x),
5785                    len(x))
5786            print '    y in [%f, %f], len(lon) == %d'\
5787                  %(min(y), max(y),
5788                    len(y))
5789            print '    z in [%f, %f], len(z) == %d'\
5790                  %(min(elevation), max(elevation),
5791                    len(elevation))
5792            print 'geo_ref: ',geo_ref
5793            print '------------------------------------------------'
5794           
5795        #z = resize(bath_grid,outfile.variables['z'][:].shape)
5796        #print "points[:,0]", points[:,0]
5797        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
5798        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
5799        outfile.variables['z'][:] = elevation
5800        outfile.variables['elevation'][:] = elevation  #FIXME HACK4
5801
5802        q = 'elevation'
5803        # This updates the _range values
5804        outfile.variables[q+Write_sts.RANGE][0] = min(elevation)
5805        outfile.variables[q+Write_sts.RANGE][1] = max(elevation)
5806
5807    def store_quantities(self, outfile, sts_precision=Float32,
5808                         slice_index=None, time=None,
5809                         verbose=False, **quant):
5810       
5811        """
5812        Write the quantity info.
5813
5814        **quant is extra keyword arguments passed in. These must be
5815          the sts quantities, currently; stage.
5816       
5817        if the time array is already been built, use the slice_index
5818        to specify the index.
5819       
5820        Otherwise, use time to increase the time dimension
5821
5822        Maybe make this general, but the viewer assumes these quantities,
5823        so maybe we don't want it general - unless the viewer is general
5824       
5825        precondition:
5826            triangulation and
5827            header have been called.
5828        """
5829        if time is not None:
5830            file_time = outfile.variables['time']
5831            slice_index = len(file_time)
5832            file_time[slice_index] = time   
5833
5834        # Write the conserved quantities from Domain.
5835        # Typically stage,  xmomentum, ymomentum
5836        # other quantities will be ignored, silently.
5837        # Also write the ranges: stage_range
5838        for q in Write_sts.sts_quantities:
5839            if not quant.has_key(q):
5840                msg = 'STS file can not write quantity %s' %q
5841                raise NewQuantity, msg
5842            else:
5843                q_values = quant[q]
5844                outfile.variables[q][slice_index] = \
5845                                q_values.astype(sts_precision)
5846
5847                # This updates the _range values
5848                q_range = outfile.variables[q+Write_sts.RANGE][:]
5849                q_values_min = min(q_values)
5850                if q_values_min < q_range[0]:
5851                    outfile.variables[q+Write_sts.RANGE][0] = q_values_min
5852                q_values_max = max(q_values)
5853                if q_values_max > q_range[1]:
5854                    outfile.variables[q+Write_sts.RANGE][1] = q_values_max
5855
5856
5857   
5858class Urs_points:
5859    """
5860    Read the info in URS mux files.
5861
5862    for the quantities heres a correlation between the file names and
5863    what they mean;
5864    z-mux is height above sea level, m
5865    e-mux is velocity is Eastern direction, m/s
5866    n-mux is velocity is Northern direction, m/s   
5867    """
5868    def __init__(self,urs_file):
5869        self.iterated = False
5870        columns = 3 # long, lat , depth
5871        mux_file = open(urs_file, 'rb')
5872       
5873        # Number of points/stations
5874        (self.points_num,)= unpack('i',mux_file.read(4))
5875       
5876        # nt, int - Number of time steps
5877        (self.time_step_count,)= unpack('i',mux_file.read(4))
5878        #print "self.time_step_count", self.time_step_count
5879        #dt, float - time step, seconds
5880        (self.time_step,) = unpack('f', mux_file.read(4))
5881        #print "self.time_step", self.time_step
5882        msg = "Bad data in the urs file."
5883        if self.points_num < 0:
5884            mux_file.close()
5885            raise ANUGAError, msg
5886        if self.time_step_count < 0:
5887            mux_file.close()
5888            raise ANUGAError, msg
5889        if self.time_step < 0:
5890            mux_file.close()
5891            raise ANUGAError, msg
5892
5893        # the depth is in meters, and it is the distance from the ocean
5894        # to the sea bottom.
5895        lonlatdep = p_array.array('f')
5896        lonlatdep.read(mux_file, columns * self.points_num)
5897        lonlatdep = array(lonlatdep, typecode=Float)   
5898        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
5899        #print 'lonlatdep',lonlatdep
5900        self.lonlatdep = lonlatdep
5901       
5902        self.mux_file = mux_file
5903        # check this array
5904
5905    def __iter__(self):
5906        """
5907        iterate over quantity data which is with respect to time.
5908
5909        Note: You can only interate once over an object
5910       
5911        returns quantity infomation for each time slice
5912        """
5913        msg =  "You can only interate once over a urs file."
5914        assert not self.iterated, msg
5915        self.iter_time_step = 0
5916        self.iterated = True
5917        return self
5918   
5919    def next(self):
5920        if self.time_step_count == self.iter_time_step:
5921            self.close()
5922            raise StopIteration
5923        #Read in a time slice  from mux file 
5924        hz_p_array = p_array.array('f')
5925        hz_p_array.read(self.mux_file, self.points_num)
5926        hz_p = array(hz_p_array, typecode=Float)
5927        self.iter_time_step += 1
5928       
5929        return hz_p
5930
5931    def close(self):
5932        self.mux_file.close()
5933       
5934    #### END URS UNGRIDDED 2 SWW ###
5935
5936       
5937def start_screen_catcher(dir_name=None, myid='', numprocs='', extra_info='',
5938                         verbose=True):
5939    """
5940    Used to store screen output and errors to file, if run on multiple
5941    processes eachprocessor will have its own output and error file.
5942   
5943    extra_info - is used as a string that can identify outputs with another
5944    string eg. '_other'
5945   
5946    FIXME: Would be good if you could suppress all the screen output and
5947    only save it to file... however it seems a bit tricky as this capture
5948    techique response to sys.stdout and by this time it is already printed out.
5949    """
5950   
5951    import sys
5952#    dir_name = dir_name
5953    if dir_name == None:
5954        dir_name=getcwd()
5955       
5956    if access(dir_name,W_OK) == 0:
5957        if verbose: print 'Making directory %s' %dir_name
5958      #  if verbose: print "myid", myid
5959        mkdir (dir_name,0777)
5960
5961    if myid <>'':
5962        myid = '_'+str(myid)
5963    if numprocs <>'':
5964        numprocs = '_'+str(numprocs)
5965    if extra_info <>'':
5966        extra_info = '_'+str(extra_info)
5967#    print 'hello1'
5968    screen_output_name = join(dir_name, "screen_output%s%s%s.txt" %(myid,
5969                                                                numprocs,
5970                                                                extra_info))
5971    screen_error_name = join(dir_name,  "screen_error%s%s%s.txt" %(myid,
5972                                                              numprocs,
5973                                                              extra_info))
5974
5975    if verbose: print 'Starting ScreenCatcher, all output will be stored in %s' \
5976                                     %(screen_output_name)
5977    #used to catch screen output to file
5978    sys.stdout = Screen_Catcher(screen_output_name)
5979    sys.stderr = Screen_Catcher(screen_error_name)
5980
5981class Screen_Catcher:
5982    """this simply catches the screen output and stores it to file defined by
5983    start_screen_catcher (above)
5984    """
5985   
5986    def __init__(self, filename):
5987        self.filename = filename
5988#        print 'init'
5989        if exists(self.filename)is True:
5990            print'Old existing file "%s" has been deleted' %(self.filename)
5991            remove(self.filename)
5992
5993    def write(self, stuff):
5994        fid = open(self.filename, 'a')
5995        fid.write(stuff)
5996        fid.close()
5997       
5998# FIXME (DSG): Add unit test, make general, not just 2 files,
5999# but any number of files.
6000def copy_code_files(dir_name, filename1, filename2=None):
6001    """Copies "filename1" and "filename2" to "dir_name". Very useful for
6002    information management
6003    filename1 and filename2 are both absolute pathnames   
6004    """
6005
6006    if access(dir_name,F_OK) == 0:
6007        print 'Make directory %s' %dir_name
6008        mkdir (dir_name,0777)
6009    shutil.copy(filename1, dir_name + sep + basename(filename1))
6010    if filename2!=None:
6011        shutil.copy(filename2, dir_name + sep + basename(filename2))
6012        print 'Files %s and %s copied' %(filename1, filename2)
6013    else:
6014        print 'File %s copied' %(filename1)
6015
6016def get_data_from_file(filename,separator_value = ','):
6017    """
6018    Read in data information from file and
6019   
6020    Returns:
6021        header_fields, a string? of the first line separated
6022        by the 'separator_value'
6023       
6024        data, a array (N data columns X M lines) in the file
6025        excluding the header
6026       
6027    NOTE: wont deal with columns with different lenghts and there must be
6028    no blank lines at the end.
6029    """
6030   
6031    fid = open(filename)
6032    lines = fid.readlines()
6033   
6034    fid.close()
6035   
6036    header_line = lines[0]
6037    header_fields = header_line.split(separator_value)
6038
6039    #array to store data, number in there is to allow float...
6040    #i'm sure there is a better way!
6041    data=array([],typecode=Float)
6042    data=resize(data,((len(lines)-1),len(header_fields)))
6043#    print 'number of fields',range(len(header_fields))
6044#    print 'number of lines',len(lines), shape(data)
6045#    print'data',data[1,1],header_line
6046
6047    array_number = 0
6048    line_number = 1
6049    while line_number < (len(lines)):
6050        for i in range(len(header_fields)): 
6051            #this get line below the header, explaining the +1
6052            #and also the line_number can be used as the array index
6053            fields = lines[line_number].split(separator_value)
6054            #assign to array
6055            data[array_number,i] = float(fields[i])
6056           
6057        line_number = line_number +1
6058        array_number = array_number +1
6059       
6060    return header_fields, data
6061
6062def store_parameters(verbose=False,**kwargs):
6063    """
6064    Store "kwargs" into a temp csv file, if "completed" is a kwargs csv file is
6065    kwargs[file_name] else it is kwargs[output_dir] + details_temp.csv
6066   
6067    Must have a file_name keyword arg, this is what is writing to.
6068    might be a better way to do this using CSV module Writer and writeDict
6069   
6070    writes file to "output_dir" unless "completed" is in kwargs, then
6071    it writes to "file_name" kwargs
6072
6073    """
6074    import types
6075#    import os
6076   
6077    # Check that kwargs is a dictionary
6078    if type(kwargs) != types.DictType:
6079        raise TypeError
6080   
6081    #is completed is kwargs?
6082    try:
6083        kwargs['completed']
6084        completed=True
6085    except:
6086        completed=False
6087 
6088    #get file name and removes from dict and assert that a file_name exists
6089    if completed:
6090        try:
6091            file = str(kwargs['file_name'])
6092        except:
6093            raise 'kwargs must have file_name'
6094    else:
6095        #write temp file in output directory
6096        try:
6097            file = str(kwargs['output_dir'])+'detail_temp.csv'
6098        except:
6099            raise 'kwargs must have output_dir'
6100       
6101    #extracts the header info and the new line info
6102    line=''
6103    header=''
6104    count=0
6105    keys = kwargs.keys()
6106    keys.sort()
6107   
6108    #used the sorted keys to create the header and line data
6109    for k in keys:
6110#        print "%s = %s" %(k, kwargs[k])
6111        header = header+str(k)
6112        line = line+str(kwargs[k])
6113        count+=1
6114        if count <len(kwargs):
6115            header = header+','
6116            line = line+','
6117    header+='\n'
6118    line+='\n'
6119
6120    # checks the header info, if the same, then write, if not create a new file
6121    #try to open!
6122    try:
6123        fid = open(file,"r")
6124        file_header=fid.readline()
6125        fid.close()
6126        if verbose: print 'read file header %s' %file_header
6127       
6128    except:
6129        msg = 'try to create new file',file
6130        if verbose: print msg
6131        #tries to open file, maybe directory is bad
6132        try:
6133            fid = open(file,"w")
6134            fid.write(header)
6135            fid.close()
6136            file_header=header
6137        except:
6138            msg = 'cannot create new file',file
6139            raise msg
6140           
6141    #if header is same or this is a new file
6142    if file_header==str(header):
6143        fid=open(file,"a")
6144        #write new line
6145        fid.write(line)
6146        fid.close()
6147    else:
6148        #backup plan,
6149        # if header is different and has completed will append info to
6150        #end of details_temp.cvs file in output directory
6151        file = str(kwargs['output_dir'])+'detail_temp.csv'
6152        fid=open(file,"a")
6153        fid.write(header)
6154        fid.write(line)
6155        fid.close()
6156        if verbose: print 'file',file_header.strip('\n')
6157        if verbose: print 'head',header.strip('\n')
6158        if file_header.strip('\n')==str(header): print 'they equal'
6159        msg = 'WARNING: File header does not match input info, the input variables have changed, suggest to change file name'
6160        print msg
6161
6162
6163
6164# ----------------------------------------------
6165# Functions to obtain diagnostics from sww files
6166#-----------------------------------------------
6167
6168def get_mesh_and_quantities_from_file(filename,
6169                                      quantities=None,
6170                                      verbose=False):
6171    """Get and rebuild mesh structure and the associated quantities from sww file
6172
6173    Input:
6174        filename - Name os sww file
6175        quantities - Names of quantities to load
6176
6177    Output:
6178        mesh - instance of class Interpolate
6179               (including mesh and interpolation functionality)
6180        quantities - arrays with quantity values at each mesh node
6181        time - vector of stored timesteps
6182    """
6183 
6184    # FIXME (Ole): Maybe refactor filefunction using this more fundamental code.
6185   
6186    import types
6187    from Scientific.IO.NetCDF import NetCDFFile
6188    from shallow_water import Domain
6189    from Numeric import asarray, transpose, resize
6190    from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
6191
6192    if verbose: print 'Reading from ', filename
6193    fid = NetCDFFile(filename, 'r')    # Open existing file for read
6194    time = fid.variables['time'][:]    # Time vector
6195    time += fid.starttime[0]
6196   
6197    # Get the variables as Numeric arrays
6198    x = fid.variables['x'][:]                   # x-coordinates of nodes
6199    y = fid.variables['y'][:]                   # y-coordinates of nodes
6200    elevation = fid.variables['elevation'][:]   # Elevation
6201    stage = fid.variables['stage'][:]           # Water level
6202    xmomentum = fid.variables['xmomentum'][:]   # Momentum in the x-direction
6203    ymomentum = fid.variables['ymomentum'][:]   # Momentum in the y-direction
6204
6205
6206    # Mesh (nodes (Mx2), triangles (Nx3))
6207    nodes = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
6208    triangles = fid.variables['volumes'][:]
6209   
6210    # Get geo_reference
6211    try:
6212        geo_reference = Geo_reference(NetCDFObject=fid)
6213    except: #AttributeError, e:
6214        # Sww files don't have to have a geo_ref
6215        geo_reference = None
6216
6217    if verbose: print '    building mesh from sww file %s' %filename
6218    boundary = None
6219
6220    #FIXME (Peter Row): Should this be in mesh?
6221    if fid.smoothing != 'Yes':
6222        nodes = nodes.tolist()
6223        triangles = triangles.tolist()
6224        nodes, triangles, boundary=weed(nodes, triangles, boundary)
6225
6226
6227    try:
6228        mesh = Mesh(nodes, triangles, boundary,
6229                    geo_reference=geo_reference)
6230    except AssertionError, e:
6231        fid.close()
6232        msg = 'Domain could not be created: %s. "' %e
6233        raise DataDomainError, msg
6234
6235
6236    quantities = {}
6237    quantities['elevation'] = elevation
6238    quantities['stage'] = stage   
6239    quantities['xmomentum'] = xmomentum
6240    quantities['ymomentum'] = ymomentum   
6241
6242    fid.close()
6243    return mesh, quantities, time
6244
6245
6246
6247def get_flow_through_cross_section(filename,
6248                                   polyline,
6249                                   verbose=False):