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

Last change on this file since 6215 was 6215, checked in by ole, 15 years ago

Removed references to _range in sww files in preparation for new viewer and ticket:192
Current extrema still needs to be removed and replaced by vertex based extrema.

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