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

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

Added functionality for clipping buildings by a list of polygons.
This is not finished and test is only rudimentary.

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