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

Last change on this file since 4554 was 4554, checked in by ole, 17 years ago

Implemented get_maximum_inundation_elevation and
get_maximum_inundation_location for use with sww files, added tests and
updated user manual.

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