source: branches/numpy/anuga/shallow_water/data_manager.py @ 7035

Last change on this file since 7035 was 7035, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk.

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