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

Last change on this file since 4558 was 4558, checked in by duncan, 17 years ago

removing stuff that made the sww files absolute. see ticket#173

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