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

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

Work on ticket:175

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