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

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

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

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