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

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

Added polygon values example file

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