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

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

Initial commit of numpy changes. Still a long way to go.

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