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

Last change on this file since 5589 was 5589, checked in by sexton, 15 years ago

unit tests for comparing urs gauges to sts outputs

File size: 221.6 KB
RevLine 
[2852]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
[4663]13.csv: ASCII format for storing arbitrary points and associated attributes
[2852]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)
[3535]33TSH:          Triangular meshes (e.g. created from anuga.pmesh)
[2852]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
[3560]47TSH + Boundary SWW -> SWW: Simluation using abstract_2d_finite_volumes
[2852]48
49"""
50
[5485]51# This file was reverted from changeset:5484 to changeset:5470 on 10th July
52# by Ole.
53
[2852]54import exceptions
[3292]55class TitleValueError(exceptions.Exception): pass
[2852]56class DataMissingValuesError(exceptions.Exception): pass
57class DataFileNotOpenError(exceptions.Exception): pass
58class DataTimeError(exceptions.Exception): pass
59class DataDomainError(exceptions.Exception): pass
[4455]60class NewQuantity(exceptions.Exception): pass
[2852]61
62
63
[3292]64import csv
[4598]65import os, sys
[4500]66import shutil
[3720]67from struct import unpack
68import array as p_array
[4271]69#import time, os
[4595]70from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
[2852]71
[4595]72
[4868]73from Numeric import concatenate, array, Float, Int, Int32, resize, \
74     sometrue, searchsorted, zeros, allclose, around, reshape, \
75     transpose, sort, NewAxis, ArrayType, compress, take, arange, \
76     argmax, alltrue, shape, Float32, size
[4595]77
78import string
79
[3720]80from Scientific.IO.NetCDF import NetCDFFile
[4500]81#from shutil import copy
[4595]82from os.path import exists, basename, join
[4636]83from os import getcwd
[2852]84
85
[4271]86from anuga.coordinate_transforms.redfearn import redfearn, \
87     convert_from_latlon_to_utm
[4387]88from anuga.coordinate_transforms.geo_reference import Geo_reference, \
[4455]89     write_NetCDF_georeference, ensure_geo_reference
[4382]90from anuga.geospatial_data.geospatial_data import Geospatial_data,\
91     ensure_absolute
[3642]92from anuga.config import minimum_storable_height as default_minimum_storable_height
[4699]93from anuga.config import max_float
[4271]94from anuga.utilities.numerical_tools import ensure_numeric,  mean
[3720]95from anuga.caching.caching import myhash
96from anuga.utilities.anuga_exceptions import ANUGAError
[4271]97from anuga.shallow_water import Domain
98from anuga.abstract_2d_finite_volumes.pmesh2domain import \
99     pmesh_to_domain_instance
[4480]100from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
[4567]101     remove_lone_verts, sww2timeseries, get_centroid_values
[4497]102from anuga.load_mesh.loadASCII import export_mesh_file
[5226]103from anuga.utilities.polygon import intersection
104
105
[4462]106# formula mappings
107
108quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
109                    'depth':'stage-elevation',
110                    'speed': \
111 '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
112
113
[4271]114   
[2852]115def make_filename(s):
[4612]116    """Transform argument string into a Sexsuitable filename
[2852]117    """
118
119    s = s.strip()
120    s = s.replace(' ', '_')
121    s = s.replace('(', '')
122    s = s.replace(')', '')
123    s = s.replace('__', '_')
124
125    return s
126
127
128def check_dir(path, verbose=None):
129    """Check that specified path exists.
130    If path does not exist it will be created if possible
131
132    USAGE:
133       checkdir(path, verbose):
134
135    ARGUMENTS:
136        path -- Directory
137        verbose -- Flag verbose output (default: None)
138
139    RETURN VALUE:
140        Verified path including trailing separator
141
142    """
143
144    import os.path
145
146    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
147        unix = 0
148    else:
149        unix = 1
150
151
152    if path[-1] != os.sep:
153        path = path + os.sep  # Add separator for directories
154
155    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
156    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
157        try:
158            exitcode=os.mkdir(path)
159
160            # Change access rights if possible
161            #
162            if unix:
163                exitcode=os.system('chmod 775 '+path)
164            else:
165                pass  # FIXME: What about acces rights under Windows?
166
167            if verbose: print 'MESSAGE: Directory', path, 'created.'
168
169        except:
170            print 'WARNING: Directory', path, 'could not be created.'
171            if unix:
172                path = '/tmp/'
173            else:
174                path = 'C:'
175
176            print 'Using directory %s instead' %path
177
178    return(path)
179
180
181
182def del_dir(path):
183    """Recursively delete directory path and all its contents
184    """
185
186    import os
187
188    if os.path.isdir(path):
189        for file in os.listdir(path):
190            X = os.path.join(path, file)
191
192
193            if os.path.isdir(X) and not os.path.islink(X):
194                del_dir(X)
195            else:
196                try:
197                    os.remove(X)
198                except:
199                    print "Could not remove file %s" %X
200
201        os.rmdir(path)
[4598]202       
203       
204# ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007   
[4603]205def rmgeneric(path, __func__,verbose=False):
206    ERROR_STR= """Error removing %(path)s, %(error)s """
[2852]207
[4603]208    try:
209        __func__(path)
210        if verbose: print 'Removed ', path
211    except OSError, (errno, strerror):
212        print ERROR_STR % {'path' : path, 'error': strerror }
213           
214def removeall(path,verbose=False):
[2852]215
[4603]216    if not os.path.isdir(path):
217        return
218   
219    files=os.listdir(path)
[2852]220
[4603]221    for x in files:
222        fullpath=os.path.join(path, x)
223        if os.path.isfile(fullpath):
224            f=os.remove
225            rmgeneric(fullpath, f)
226        elif os.path.isdir(fullpath):
227            removeall(fullpath)
228            f=os.rmdir
229            rmgeneric(fullpath, f,verbose)
230
231
232
[2852]233def create_filename(datadir, filename, format, size=None, time=None):
234
235    import os
[3514]236    #from anuga.config import data_dir
[2852]237
238    FN = check_dir(datadir) + filename
239
240    if size is not None:
241        FN += '_size%d' %size
242
243    if time is not None:
244        FN += '_time%.2f' %time
245
246    FN += '.' + format
247    return FN
248
249
250def get_files(datadir, filename, format, size):
[4595]251    """Get all file (names) with given name, size and format
[2852]252    """
253
254    import glob
255
256    import os
[3514]257    #from anuga.config import data_dir
[2852]258
259    dir = check_dir(datadir)
260
261    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
262    return glob.glob(pattern)
263
264
265
266#Generic class for storing output to e.g. visualisation or checkpointing
267class Data_format:
268    """Generic interface to data formats
269    """
270
271
272    def __init__(self, domain, extension, mode = 'w'):
273        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
274                                        '''   'w' (write)'''+\
275                                        '''   'r' (read)''' +\
276                                        '''   'a' (append)'''
277
278        #Create filename
279        self.filename = create_filename(domain.get_datadir(),
[3928]280                                        domain.get_name(), extension)
[2852]281
282        #print 'F', self.filename
283        self.timestep = 0
284        self.domain = domain
[3928]285       
[2852]286
287
[3928]288        # Exclude ghosts in case this is a parallel domain
289        self.number_of_nodes = domain.number_of_full_nodes       
290        self.number_of_volumes = domain.number_of_full_triangles
291        #self.number_of_volumes = len(domain)       
292
293
294
295
[2852]296        #FIXME: Should we have a general set_precision function?
297
298
299
300#Class for storing output to e.g. visualisation
301class Data_format_sww(Data_format):
302    """Interface to native NetCDF format (.sww) for storing model output
303
304    There are two kinds of data
305
306    1: Constant data: Vertex coordinates and field values. Stored once
307    2: Variable data: Conserved quantities. Stored once per timestep.
308
309    All data is assumed to reside at vertex locations.
310    """
311
312
313    def __init__(self, domain, mode = 'w',\
314                 max_size = 2000000000,
315                 recursion = False):
316        from Scientific.IO.NetCDF import NetCDFFile
317        from Numeric import Int, Float, Float32
318
319        self.precision = Float32 #Use single precision for quantities
320        if hasattr(domain, 'max_size'):
321            self.max_size = domain.max_size #file size max is 2Gig
322        else:
323            self.max_size = max_size
324        self.recursion = recursion
325        self.mode = mode
326
327        Data_format.__init__(self, domain, 'sww', mode)
328
[3642]329        if hasattr(domain, 'minimum_storable_height'):
[4704]330            self.minimum_storable_height = domain.minimum_storable_height
[3529]331        else:
[3642]332            self.minimum_storable_height = default_minimum_storable_height
[2852]333
334        # NetCDF file definition
335        fid = NetCDFFile(self.filename, mode)
336
337        if mode == 'w':
[4455]338            description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
339            self.writer = Write_sww()
[4704]340            self.writer.store_header(fid,
341                                     domain.starttime,
342                                     self.number_of_volumes,
343                                     self.domain.number_of_full_nodes,
344                                     description=description,
345                                     smoothing=domain.smooth,
[4862]346                                     order=domain.default_order,
347                                     sww_precision=self.precision)
[4704]348
349            # Extra optional information
[4455]350            if hasattr(domain, 'texture'):
[4704]351                fid.texture = domain.texture
[2852]352
[4704]353            if domain.quantities_to_be_monitored is not None:
354                fid.createDimension('singleton', 1)
[4705]355                fid.createDimension('two', 2)               
356
357                poly = domain.monitor_polygon
358                if poly is not None:
359                    N = len(poly)
360                    fid.createDimension('polygon_length', N)
[4706]361                    fid.createVariable('extrema.polygon',
[4705]362                                       self.precision,
363                                       ('polygon_length',
364                                        'two'))
[4706]365                    fid.variables['extrema.polygon'][:] = poly                                   
366
[4705]367                   
368                interval = domain.monitor_time_interval
369                if interval is not None:
[4706]370                    fid.createVariable('extrema.time_interval',
[4705]371                                       self.precision,
372                                       ('two',))
[4706]373                    fid.variables['extrema.time_interval'][:] = interval
[4705]374
375               
[4704]376                for q in domain.quantities_to_be_monitored:
377                    #print 'doing', q
[4706]378                    fid.createVariable(q+'.extrema', self.precision,
[4704]379                                       ('numbers_in_range',))
[4706]380                    fid.createVariable(q+'.min_location', self.precision,
[4704]381                                       ('numbers_in_range',))
[4706]382                    fid.createVariable(q+'.max_location', self.precision,
[4704]383                                       ('numbers_in_range',))
[4706]384                    fid.createVariable(q+'.min_time', self.precision,
[4704]385                                       ('singleton',))
[4706]386                    fid.createVariable(q+'.max_time', self.precision,
[4704]387                                       ('singleton',))
[2852]388
[4704]389                   
[2852]390        fid.close()
391
392
393    def store_connectivity(self):
394        """Specialisation of store_connectivity for net CDF format
395
396        Writes x,y,z coordinates of triangles constituting
397        the bed elevation.
398        """
399
400        from Scientific.IO.NetCDF import NetCDFFile
401
402        from Numeric import concatenate, Int
403
404        domain = self.domain
405
406        #Get NetCDF
407        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
408
409        # Get the variables
410        x = fid.variables['x']
411        y = fid.variables['y']
412        z = fid.variables['elevation']
413
414        volumes = fid.variables['volumes']
415
416        # Get X, Y and bed elevation Z
417        Q = domain.quantities['elevation']
418        X,Y,Z,V = Q.get_vertex_values(xy=True,
[3926]419                                      precision=self.precision)
[2852]420
[4455]421        #
422        points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
[4704]423        self.writer.store_triangulation(fid,
424                                        points,
425                                        V.astype(volumes.typecode()),
426                                        Z,
427                                        points_georeference= \
428                                        domain.geo_reference)
[2852]429
[4704]430        # Close
431        fid.close()
[2852]432
433
[4868]434    def store_timestep(self, names=None):
[2852]435        """Store time and named quantities to file
436        """
[4868]437       
[2852]438        from Scientific.IO.NetCDF import NetCDFFile
439        import types
440        from time import sleep
441        from os import stat
442
443        from Numeric import choose
[4868]444
445
446        if names is None:
447            # Standard shallow water wave equation quantitites in ANUGA
448            names = ['stage', 'xmomentum', 'ymomentum']
[4455]449       
[4704]450        # Get NetCDF       
[2852]451        retries = 0
452        file_open = False
453        while not file_open and retries < 10:
454            try:
[4704]455                fid = NetCDFFile(self.filename, 'a') # Open existing file
[2852]456            except IOError:
[4704]457                # This could happen if someone was reading the file.
458                # In that case, wait a while and try again
[2852]459                msg = 'Warning (store_timestep): File %s could not be opened'\
460                      %self.filename
461                msg += ' - trying step %s again' %self.domain.time
462                print msg
463                retries += 1
464                sleep(1)
465            else:
466                file_open = True
467
468        if not file_open:
469            msg = 'File %s could not be opened for append' %self.filename
470            raise DataFileNotOpenError, msg
471
472
473
[4704]474        # Check to see if the file is already too big:
[2852]475        time = fid.variables['time']
476        i = len(time)+1
477        file_size = stat(self.filename)[6]
478        file_size_increase =  file_size/i
479        if file_size + file_size_increase > self.max_size*(2**self.recursion):
[4704]480            # In order to get the file name and start time correct,
481            # I change the domain.filename and domain.starttime.
482            # This is the only way to do this without changing
483            # other modules (I think).
[2852]484
[4704]485            # Write a filename addon that won't break swollens reader
486            # (10.sww is bad)
[2852]487            filename_ext = '_time_%s'%self.domain.time
488            filename_ext = filename_ext.replace('.', '_')
[4704]489           
490            # Remember the old filename, then give domain a
491            # name with the extension
[3846]492            old_domain_filename = self.domain.get_name()
[2852]493            if not self.recursion:
[3846]494                self.domain.set_name(old_domain_filename+filename_ext)
[2852]495
[3846]496
[4704]497            # Change the domain starttime to the current time
[2852]498            old_domain_starttime = self.domain.starttime
499            self.domain.starttime = self.domain.time
500
[4704]501            # Build a new data_structure.
[2852]502            next_data_structure=\
503                Data_format_sww(self.domain, mode=self.mode,\
504                                max_size = self.max_size,\
505                                recursion = self.recursion+1)
506            if not self.recursion:
507                print '    file_size = %s'%file_size
508                print '    saving file to %s'%next_data_structure.filename
509            #set up the new data_structure
510            self.domain.writer = next_data_structure
511
512            #FIXME - could be cleaner to use domain.store_timestep etc.
513            next_data_structure.store_connectivity()
514            next_data_structure.store_timestep(names)
515            fid.sync()
516            fid.close()
517
518            #restore the old starttime and filename
519            self.domain.starttime = old_domain_starttime
[3846]520            self.domain.set_name(old_domain_filename)           
[2852]521        else:
522            self.recursion = False
523            domain = self.domain
524
525            # Get the variables
526            time = fid.variables['time']
527            stage = fid.variables['stage']
528            xmomentum = fid.variables['xmomentum']
529            ymomentum = fid.variables['ymomentum']
530            i = len(time)
531            if type(names) not in [types.ListType, types.TupleType]:
532                names = [names]
533
[4455]534            if 'stage' in names and 'xmomentum' in names and \
[4704]535               'ymomentum' in names:
[4455]536
[4868]537                # Get stage, elevation, depth and select only those
538                # values where minimum_storable_height is exceeded
[4455]539                Q = domain.quantities['stage']
[4868]540                A, _ = Q.get_vertex_values(xy = False,
541                                           precision = self.precision)
[4455]542                z = fid.variables['elevation']
[4868]543
544                storable_indices = A-z[:] >= self.minimum_storable_height
545                stage = choose(storable_indices, (z[:], A))
[4455]546               
[4868]547                # Define a zero vector of same size and type as A
548                # for use with momenta
549                null = zeros(size(A), A.typecode())
550               
551                # Get xmomentum where depth exceeds minimum_storable_height
[4455]552                Q = domain.quantities['xmomentum']
[4868]553                xmom, _ = Q.get_vertex_values(xy = False,
554                                              precision = self.precision)
555                xmomentum = choose(storable_indices, (null, xmom))
[4455]556               
[4868]557
558                # Get ymomentum where depth exceeds minimum_storable_height
[4455]559                Q = domain.quantities['ymomentum']
[4868]560                ymom, _ = Q.get_vertex_values(xy = False,
561                                              precision = self.precision)
562                ymomentum = choose(storable_indices, (null, ymom))               
[4704]563               
564                # Write quantities to NetCDF
565                self.writer.store_quantities(fid, 
566                                             time=self.domain.time,
[4862]567                                             sww_precision=self.precision,
[4704]568                                             stage=stage,
569                                             xmomentum=xmomentum,
570                                             ymomentum=ymomentum)
[4455]571            else:
[4868]572                msg = 'Quantities stored must be: stage, xmomentum, ymomentum.'
573                msg += ' Instead I got: ' + str(names)
574                raise Exception, msg
575           
[2852]576
[3946]577
[4704]578            # Update extrema if requested
579            domain = self.domain
580            if domain.quantities_to_be_monitored is not None:
581                for q, info in domain.quantities_to_be_monitored.items():
582
583                    if info['min'] is not None:
[4706]584                        fid.variables[q + '.extrema'][0] = info['min']
585                        fid.variables[q + '.min_location'][:] =\
[4704]586                                        info['min_location']
[4706]587                        fid.variables[q + '.min_time'][0] = info['min_time']
[4704]588                       
589                    if info['max'] is not None:
[4706]590                        fid.variables[q + '.extrema'][1] = info['max']
591                        fid.variables[q + '.max_location'][:] =\
[4704]592                                        info['max_location']
[4706]593                        fid.variables[q + '.max_time'][0] = info['max_time']
[4704]594
595           
596
[4868]597            # Flush and close
[2852]598            fid.sync()
599            fid.close()
600
601
602
[4704]603# Class for handling checkpoints data
[2852]604class Data_format_cpt(Data_format):
605    """Interface to native NetCDF format (.cpt)
606    """
607
608
609    def __init__(self, domain, mode = 'w'):
610        from Scientific.IO.NetCDF import NetCDFFile
611        from Numeric import Int, Float, Float
612
613        self.precision = Float #Use full precision
614
615        Data_format.__init__(self, domain, 'sww', mode)
616
617
618        # NetCDF file definition
619        fid = NetCDFFile(self.filename, mode)
620
621        if mode == 'w':
622            #Create new file
623            fid.institution = 'Geoscience Australia'
624            fid.description = 'Checkpoint data'
625            #fid.smooth = domain.smooth
626            fid.order = domain.default_order
627
628            # dimension definitions
629            fid.createDimension('number_of_volumes', self.number_of_volumes)
630            fid.createDimension('number_of_vertices', 3)
631
632            #Store info at all vertices (no smoothing)
633            fid.createDimension('number_of_points', 3*self.number_of_volumes)
634            fid.createDimension('number_of_timesteps', None) #extensible
635
636            # variable definitions
637
638            #Mesh
639            fid.createVariable('x', self.precision, ('number_of_points',))
640            fid.createVariable('y', self.precision, ('number_of_points',))
641
642
643            fid.createVariable('volumes', Int, ('number_of_volumes',
644                                                'number_of_vertices'))
645
646            fid.createVariable('time', self.precision,
647                               ('number_of_timesteps',))
648
649            #Allocate space for all quantities
650            for name in domain.quantities.keys():
651                fid.createVariable(name, self.precision,
652                                   ('number_of_timesteps',
653                                    'number_of_points'))
654
655        #Close
656        fid.close()
657
658
659    def store_checkpoint(self):
660        """
661        Write x,y coordinates of triangles.
662        Write connectivity (
663        constituting
664        the bed elevation.
665        """
666
667        from Scientific.IO.NetCDF import NetCDFFile
668
669        from Numeric import concatenate
670
671        domain = self.domain
672
673        #Get NetCDF
674        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
675
676        # Get the variables
677        x = fid.variables['x']
678        y = fid.variables['y']
679
680        volumes = fid.variables['volumes']
681
682        # Get X, Y and bed elevation Z
683        Q = domain.quantities['elevation']
684        X,Y,Z,V = Q.get_vertex_values(xy=True,
685                      precision = self.precision)
686
687
688
689        x[:] = X.astype(self.precision)
690        y[:] = Y.astype(self.precision)
691        z[:] = Z.astype(self.precision)
692
693        volumes[:] = V
694
695        #Close
696        fid.close()
697
698
699    def store_timestep(self, name):
700        """Store time and named quantity to file
701        """
702        from Scientific.IO.NetCDF import NetCDFFile
703        from time import sleep
704
705        #Get NetCDF
706        retries = 0
707        file_open = False
708        while not file_open and retries < 10:
709            try:
710                fid = NetCDFFile(self.filename, 'a') #Open existing file
711            except IOError:
712                #This could happen if someone was reading the file.
713                #In that case, wait a while and try again
714                msg = 'Warning (store_timestep): File %s could not be opened'\
715                  %self.filename
716                msg += ' - trying again'
717                print msg
718                retries += 1
719                sleep(1)
720            else:
721                file_open = True
722
723        if not file_open:
724            msg = 'File %s could not be opened for append' %self.filename
725            raise DataFileNotOPenError, msg
726
727
728        domain = self.domain
729
730        # Get the variables
731        time = fid.variables['time']
732        stage = fid.variables['stage']
733        i = len(time)
734
735        #Store stage
736        time[i] = self.domain.time
737
738        # Get quantity
739        Q = domain.quantities[name]
740        A,V = Q.get_vertex_values(xy=False,
[4704]741                                  precision = self.precision)
[2852]742
743        stage[i,:] = A.astype(self.precision)
744
745        #Flush and close
746        fid.sync()
747        fid.close()
748
749
[4326]750#### NED is national exposure database (name changed to NEXIS)
[3292]751   
752LAT_TITLE = 'LATITUDE'
753LONG_TITLE = 'LONGITUDE'
[3336]754X_TITLE = 'x'
755Y_TITLE = 'y'
[3292]756class Exposure_csv:
757    def __init__(self,file_name, latitude_title=LAT_TITLE,
[3398]758                 longitude_title=LONG_TITLE, is_x_y_locations=None,
[3336]759                 x_title=X_TITLE, y_title=Y_TITLE,
[4326]760                 refine_polygon=None, title_check_list=None):
[3292]761        """
[3296]762        This class is for handling the exposure csv file.
763        It reads the file in and converts the lats and longs to a geospatial
764        data object.
765        Use the methods to read and write columns.
766
767        The format of the csv files it reads is;
768           The first row is a title row.
769           comma's are the delimiters
770           each column is a 'set' of data
771
772        Feel free to use/expand it to read other csv files.
773           
774           
775        It is not for adding and deleting rows
776       
[3292]777        Can geospatial handle string attributes? It's not made for them.
778        Currently it can't load and save string att's.
779
780        So just use geospatial to hold the x, y and georef? Bad, since
781        different att's are in diferent structures.  Not so bad, the info
782        to write if the .csv file is saved is in attribute_dic
783
784        The location info is in the geospatial attribute.
785       
786       
787        """
788        self._file_name = file_name
789        self._geospatial = None #
790
791        # self._attribute_dic is a dictionary.
792        #The keys are the column titles.
793        #The values are lists of column data
794       
795        # self._title_index_dic is a dictionary.
796        #The keys are the column titles.
797        #The values are the index positions of file columns.
798        self._attribute_dic, self._title_index_dic = \
[4612]799            csv2dict(self._file_name, title_check_list=title_check_list)
[3292]800        try:
[4059]801            #Have code here that handles caps or lower
[3292]802            lats = self._attribute_dic[latitude_title]
803            longs = self._attribute_dic[longitude_title]
804           
805        except KeyError:
806            # maybe a warning..
[3398]807            #Let's see if this works..
808            if False != is_x_y_locations:
809                is_x_y_locations = True
[3292]810            pass
811        else:
812            self._geospatial = Geospatial_data(latitudes = lats,
813                 longitudes = longs)
[3336]814
[3398]815        if is_x_y_locations is True:
[3336]816            if self._geospatial is not None:
817                pass #fixme throw an error
818            try:
819                xs = self._attribute_dic[x_title]
820                ys = self._attribute_dic[y_title]
821                points = [[float(i),float(j)] for i,j in map(None,xs,ys)]
822            except KeyError:
823                # maybe a warning..
[3664]824                msg = "Could not find location information."
[3336]825                raise TitleValueError, msg
826            else:
827                self._geospatial = Geospatial_data(data_points=points)
[3292]828           
829        # create a list of points that are in the refining_polygon
830        # described by a list of indexes representing the points
831
832    def __cmp__(self, other):
833        #print "self._attribute_dic",self._attribute_dic
834        #print "other._attribute_dic",other._attribute_dic
835        #print "self._title_index_dic", self._title_index_dic
836        #print "other._title_index_dic", other._title_index_dic
837       
838        #check that a is an instance of this class
839        if isinstance(self, type(other)):
840            result = cmp(self._attribute_dic, other._attribute_dic)
841            if result <>0:
842                return result
843            # The order of the columns is important. Therefore..
844            result = cmp(self._title_index_dic, other._title_index_dic)
845            if result <>0:
846                return result
847            for self_ls, other_ls in map(None,self._attribute_dic, \
848                    other._attribute_dic):
849                result = cmp(self._attribute_dic[self_ls],
850                             other._attribute_dic[other_ls])
851                if result <>0:
852                    return result
853            return 0
854        else:
855            return 1
856   
857
858    def get_column(self, column_name, use_refind_polygon=False):
859        """
860        Given a column name return a list of the column values
861
862        Note, the type of the values will be String!
[3437]863        do this to change a list of strings to a list of floats
864        time = [float(x) for x in time]
[3292]865       
866        Not implemented:
867        if use_refind_polygon is True, only return values in the
868        refined polygon
869        """
870        if not self._attribute_dic.has_key(column_name):
871            msg = 'Therer is  no column called %s!' %column_name
872            raise TitleValueError, msg
873        return self._attribute_dic[column_name]
874
[3437]875
876    def get_value(self, value_column_name,
877                  known_column_name,
878                  known_values,
879                  use_refind_polygon=False):
880        """
881        Do linear interpolation on the known_colum, using the known_value,
882        to return a value of the column_value_name.
883        """
884        pass
885
886
[3292]887    def get_location(self, use_refind_polygon=False):
888        """
889        Return a geospatial object which describes the
890        locations of the location file.
891
892        Note, if there is not location info, this returns None.
893       
894        Not implemented:
895        if use_refind_polygon is True, only return values in the
896        refined polygon
897        """
898        return self._geospatial
899
900    def set_column(self, column_name, column_values, overwrite=False):
901        """
902        Add a column to the 'end' (with the right most column being the end)
903        of the csv file.
904
905        Set overwrite to True if you want to overwrite a column.
906       
907        Note, in column_name white space is removed and case is not checked.
908        Precondition
909        The column_name and column_values cannot have comma's in it.
910        """
911       
912        value_row_count = \
913            len(self._attribute_dic[self._title_index_dic.keys()[0]])
914        if len(column_values) <> value_row_count: 
915            msg = 'The number of column values must equal the number of rows.'
916            raise DataMissingValuesError, msg
917       
918        if self._attribute_dic.has_key(column_name):
919            if not overwrite:
920                msg = 'Column name %s already in use!' %column_name
921                raise TitleValueError, msg
922        else:
923            # New title.  Add it to the title index.
924            self._title_index_dic[column_name] = len(self._title_index_dic)
925        self._attribute_dic[column_name] = column_values
926        #print "self._title_index_dic[column_name]",self._title_index_dic
927
928    def save(self, file_name=None):
929        """
930        Save the exposure csv file
931        """
932        if file_name is None:
933            file_name = self._file_name
934       
935        fd = open(file_name,'wb')
936        writer = csv.writer(fd)
937       
938        #Write the title to a cvs file
939        line = [None]* len(self._title_index_dic)
940        for title in self._title_index_dic.iterkeys():
941            line[self._title_index_dic[title]]= title
942        writer.writerow(line)
943       
944        # Write the values to a cvs file
945        value_row_count = \
946            len(self._attribute_dic[self._title_index_dic.keys()[0]])
947        for row_i in range(value_row_count):
948            line = [None]* len(self._title_index_dic)
949            for title in self._title_index_dic.iterkeys():
950                line[self._title_index_dic[title]]= \
951                     self._attribute_dic[title][row_i]
952            writer.writerow(line)
953
954
[5586]955def csv2array(file_name):
956    """Convert CSV files of the form
957   
958    time, discharge, velocity
959    0.0,  1.2,       0.0
960    0.1,  3.2,       1.1
961    ...
962   
963    to a dictionary of numeric arrays.
964   
965   
966    See underlying function csv2dict for more details.
967   
968    """
969   
970   
971    X, _ = csv2dict(file_name)
972   
973    Y = {}
974    for key in X.keys():
975        Y[key] = array([float(x) for x in X[key]])
976       
977    return Y   
978   
979           
[4612]980def csv2dict(file_name, title_check_list=None):
981    """
982    Load in the csv as a dic, title as key and column info as value, .
983    Also, create a dic, title as key and column index as value,
984    to keep track of the column order.
[4775]985
986    Two dictionaries are returned.
[4612]987   
[5586]988    WARNING: Values are returned as strings.
[4612]989    do this to change a list of strings to a list of floats
990        time = [float(x) for x in time]
[4775]991
992       
[4612]993    """
994   
995    #
996    attribute_dic = {}
997    title_index_dic = {}
998    titles_stripped = [] # list of titles
999    reader = csv.reader(file(file_name))
1000
1001    # Read in and manipulate the title info
1002    titles = reader.next()
1003    for i,title in enumerate(titles):
1004        titles_stripped.append(title.strip())
1005        title_index_dic[title.strip()] = i
1006    title_count = len(titles_stripped)       
1007    #print "title_index_dic",title_index_dic
1008    if title_check_list is not None:
1009        for title_check in title_check_list:
1010            #msg = "Reading error.  This row is not present ", title_check
1011            #assert title_index_dic.has_key(title_check), msg
1012            if not title_index_dic.has_key(title_check):
1013                #reader.close()
1014                msg = "Reading error.  This row is not present ", \
1015                      title_check                     
1016                raise IOError, msg
1017               
1018   
1019   
1020    #create a dic of colum values, indexed by column title
1021    for line in reader:
1022        if len(line) <> title_count:
1023            raise IOError #FIXME make this nicer
1024        for i, value in enumerate(line):
1025            attribute_dic.setdefault(titles_stripped[i],[]).append(value)
1026       
1027    return attribute_dic, title_index_dic
1028
1029
[2852]1030#Auxiliary
1031def write_obj(filename,x,y,z):
1032    """Store x,y,z vectors into filename (obj format)
1033       Vectors are assumed to have dimension (M,3) where
1034       M corresponds to the number elements.
1035       triangles are assumed to be disconnected
1036
1037       The three numbers in each vector correspond to three vertices,
1038
1039       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
1040
1041    """
1042    #print 'Writing obj to %s' % filename
1043
1044    import os.path
1045
1046    root, ext = os.path.splitext(filename)
1047    if ext == '.obj':
1048        FN = filename
1049    else:
1050        FN = filename + '.obj'
1051
1052
1053    outfile = open(FN, 'wb')
1054    outfile.write("# Triangulation as an obj file\n")
1055
1056    M, N = x.shape
1057    assert N==3  #Assuming three vertices per element
1058
1059    for i in range(M):
1060        for j in range(N):
1061            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
1062
1063    for i in range(M):
1064        base = i*N
1065        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
1066
1067    outfile.close()
1068
1069
1070#########################################################
1071#Conversion routines
1072########################################################
1073
1074def sww2obj(basefilename, size):
1075    """Convert netcdf based data output to obj
1076    """
1077    from Scientific.IO.NetCDF import NetCDFFile
1078
1079    from Numeric import Float, zeros
1080
1081    #Get NetCDF
1082    FN = create_filename('.', basefilename, 'sww', size)
1083    print 'Reading from ', FN
1084    fid = NetCDFFile(FN, 'r')  #Open existing file for read
1085
1086
1087    # Get the variables
1088    x = fid.variables['x']
1089    y = fid.variables['y']
1090    z = fid.variables['elevation']
1091    time = fid.variables['time']
1092    stage = fid.variables['stage']
1093
1094    M = size  #Number of lines
1095    xx = zeros((M,3), Float)
1096    yy = zeros((M,3), Float)
1097    zz = zeros((M,3), Float)
1098
1099    for i in range(M):
1100        for j in range(3):
1101            xx[i,j] = x[i+j*M]
1102            yy[i,j] = y[i+j*M]
1103            zz[i,j] = z[i+j*M]
1104
1105    #Write obj for bathymetry
1106    FN = create_filename('.', basefilename, 'obj', size)
1107    write_obj(FN,xx,yy,zz)
1108
1109
1110    #Now read all the data with variable information, combine with
1111    #x,y info and store as obj
1112
1113    for k in range(len(time)):
1114        t = time[k]
1115        print 'Processing timestep %f' %t
1116
1117        for i in range(M):
1118            for j in range(3):
1119                zz[i,j] = stage[k,i+j*M]
1120
1121
1122        #Write obj for variable data
1123        #FN = create_filename(basefilename, 'obj', size, time=t)
1124        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
1125        write_obj(FN,xx,yy,zz)
1126
1127
1128def dat2obj(basefilename):
1129    """Convert line based data output to obj
1130    FIXME: Obsolete?
1131    """
1132
1133    import glob, os
[3514]1134    from anuga.config import data_dir
[2852]1135
1136
1137    #Get bathymetry and x,y's
1138    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
1139
1140    from Numeric import zeros, Float
1141
1142    M = len(lines)  #Number of lines
1143    x = zeros((M,3), Float)
1144    y = zeros((M,3), Float)
1145    z = zeros((M,3), Float)
1146
1147    ##i = 0
1148    for i, line in enumerate(lines):
1149        tokens = line.split()
1150        values = map(float,tokens)
1151
1152        for j in range(3):
1153            x[i,j] = values[j*3]
1154            y[i,j] = values[j*3+1]
1155            z[i,j] = values[j*3+2]
1156
1157        ##i += 1
1158
1159
1160    #Write obj for bathymetry
1161    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
1162
1163
1164    #Now read all the data files with variable information, combine with
1165    #x,y info
1166    #and store as obj
1167
1168    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
1169
1170    for filename in files:
1171        print 'Processing %s' % filename
1172
1173        lines = open(data_dir+os.sep+filename,'r').readlines()
1174        assert len(lines) == M
1175        root, ext = os.path.splitext(filename)
1176
1177        #Get time from filename
1178        i0 = filename.find('_time=')
1179        if i0 == -1:
1180            #Skip bathymetry file
1181            continue
1182
1183        i0 += 6  #Position where time starts
1184        i1 = filename.find('.dat')
1185
1186        if i1 > i0:
1187            t = float(filename[i0:i1])
1188        else:
1189            raise DataTimeError, 'Hmmmm'
1190
1191
1192
1193        ##i = 0
1194        for i, line in enumerate(lines):
1195            tokens = line.split()
1196            values = map(float,tokens)
1197
1198            for j in range(3):
1199                z[i,j] = values[j]
1200
1201            ##i += 1
1202
1203        #Write obj for variable data
1204        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
1205
1206
1207def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
1208    """Read netcdf filename1, pick timesteps first:step:last and save to
1209    nettcdf file filename2
1210    """
1211    from Scientific.IO.NetCDF import NetCDFFile
1212
1213    #Get NetCDF
1214    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
1215    outfile = NetCDFFile(filename2, 'w')  #Open new file
1216
1217
1218    #Copy dimensions
1219    for d in infile.dimensions:
1220        outfile.createDimension(d, infile.dimensions[d])
1221
1222    for name in infile.variables:
1223        var = infile.variables[name]
1224        outfile.createVariable(name, var.typecode(), var.dimensions)
1225
1226
1227    #Copy the static variables
1228    for name in infile.variables:
1229        if name == 'time' or name == 'stage':
1230            pass
1231        else:
1232            #Copy
1233            outfile.variables[name][:] = infile.variables[name][:]
1234
1235    #Copy selected timesteps
1236    time = infile.variables['time']
1237    stage = infile.variables['stage']
1238
1239    newtime = outfile.variables['time']
1240    newstage = outfile.variables['stage']
1241
1242    if last is None:
1243        last = len(time)
1244
1245    selection = range(first, last, step)
1246    for i, j in enumerate(selection):
1247        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
1248        newtime[i] = time[j]
1249        newstage[i,:] = stage[j,:]
1250
1251    #Close
1252    infile.close()
1253    outfile.close()
1254
1255
1256#Get data objects
1257def get_dataobject(domain, mode='w'):
1258    """Return instance of class of given format using filename
1259    """
1260
1261    cls = eval('Data_format_%s' %domain.format)
1262    return cls(domain, mode)
1263
1264
1265
1266
1267def dem2pts(basename_in, basename_out=None,
1268            easting_min=None, easting_max=None,
1269            northing_min=None, northing_max=None,
1270            use_cache=False, verbose=False,):
1271    """Read Digitial Elevation model from the following NetCDF format (.dem)
1272
1273    Example:
1274
1275    ncols         3121
1276    nrows         1800
1277    xllcorner     722000
1278    yllcorner     5893000
1279    cellsize      25
1280    NODATA_value  -9999
1281    138.3698 137.4194 136.5062 135.5558 ..........
1282
1283    Convert to NetCDF pts format which is
1284
1285    points:  (Nx2) Float array
1286    elevation: N Float array
1287    """
1288
1289
1290
1291    kwargs = {'basename_out': basename_out,
1292              'easting_min': easting_min,
1293              'easting_max': easting_max,
1294              'northing_min': northing_min,
1295              'northing_max': northing_max,
1296              'verbose': verbose}
1297
1298    if use_cache is True:
1299        from caching import cache
1300        result = cache(_dem2pts, basename_in, kwargs,
1301                       dependencies = [basename_in + '.dem'],
1302                       verbose = verbose)
1303
1304    else:
1305        result = apply(_dem2pts, [basename_in], kwargs)
1306
1307    return result
1308
1309
1310def _dem2pts(basename_in, basename_out=None, verbose=False,
1311            easting_min=None, easting_max=None,
1312            northing_min=None, northing_max=None):
1313    """Read Digitial Elevation model from the following NetCDF format (.dem)
1314
1315    Internal function. See public function dem2pts for details.
1316    """
1317
[4776]1318    # FIXME: Can this be written feasibly using write_pts?
[2852]1319
1320    import os
1321    from Scientific.IO.NetCDF import NetCDFFile
1322    from Numeric import Float, zeros, reshape, sum
1323
1324    root = basename_in
1325
[4776]1326    # Get NetCDF
1327    infile = NetCDFFile(root + '.dem', 'r')  # Open existing netcdf file for read
[2852]1328
1329    if verbose: print 'Reading DEM from %s' %(root + '.dem')
1330
1331    ncols = infile.ncols[0]
1332    nrows = infile.nrows[0]
[4776]1333    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
1334    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
[2852]1335    cellsize = infile.cellsize[0]
1336    NODATA_value = infile.NODATA_value[0]
1337    dem_elevation = infile.variables['elevation']
1338
1339    zone = infile.zone[0]
1340    false_easting = infile.false_easting[0]
1341    false_northing = infile.false_northing[0]
1342
[4776]1343    # Text strings
[2852]1344    projection = infile.projection
1345    datum = infile.datum
1346    units = infile.units
1347
1348
[4776]1349    # Get output file
[2852]1350    if basename_out == None:
1351        ptsname = root + '.pts'
1352    else:
1353        ptsname = basename_out + '.pts'
1354
1355    if verbose: print 'Store to NetCDF file %s' %ptsname
1356    # NetCDF file definition
1357    outfile = NetCDFFile(ptsname, 'w')
1358
[4776]1359    # Create new file
[2852]1360    outfile.institution = 'Geoscience Australia'
1361    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
1362                      'of spatial point data'
[4776]1363    # Assign default values
[2852]1364    if easting_min is None: easting_min = xllcorner
1365    if easting_max is None: easting_max = xllcorner + ncols*cellsize
1366    if northing_min is None: northing_min = yllcorner
1367    if northing_max is None: northing_max = yllcorner + nrows*cellsize
1368
[4776]1369    # Compute offsets to update georeferencing
[2852]1370    easting_offset = xllcorner - easting_min
1371    northing_offset = yllcorner - northing_min
1372
[4776]1373    # Georeferencing
[2852]1374    outfile.zone = zone
[4776]1375    outfile.xllcorner = easting_min # Easting of lower left corner
1376    outfile.yllcorner = northing_min # Northing of lower left corner
[2852]1377    outfile.false_easting = false_easting
1378    outfile.false_northing = false_northing
1379
1380    outfile.projection = projection
1381    outfile.datum = datum
1382    outfile.units = units
1383
1384
[4776]1385    # Grid info (FIXME: probably not going to be used, but heck)
[2852]1386    outfile.ncols = ncols
1387    outfile.nrows = nrows
1388
1389    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
1390    totalnopoints = nrows*ncols
1391
[4776]1392    # Calculating number of NODATA_values for each row in clipped region
1393    # FIXME: use array operations to do faster
[2852]1394    nn = 0
1395    k = 0
1396    i1_0 = 0
1397    j1_0 = 0
1398    thisj = 0
1399    thisi = 0
1400    for i in range(nrows):
1401        y = (nrows-i-1)*cellsize + yllcorner
1402        for j in range(ncols):
1403            x = j*cellsize + xllcorner
1404            if easting_min <= x <= easting_max and \
1405               northing_min <= y <= northing_max:
1406                thisj = j
1407                thisi = i
1408                if dem_elevation_r[i,j] == NODATA_value: nn += 1
1409
1410                if k == 0:
1411                    i1_0 = i
1412                    j1_0 = j
1413                k += 1
1414
1415    index1 = j1_0
1416    index2 = thisj
1417
[4776]1418    # Dimension definitions
[2852]1419    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
1420    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
1421
1422    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
1423    nopoints = clippednopoints-nn
1424
1425    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
1426
[3664]1427    if verbose:
[2852]1428        print 'There are %d values in the elevation' %totalnopoints
1429        print 'There are %d values in the clipped elevation' %clippednopoints
1430        print 'There are %d NODATA_values in the clipped elevation' %nn
1431
1432    outfile.createDimension('number_of_points', nopoints)
1433    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1434
[4776]1435    # Variable definitions
[2852]1436    outfile.createVariable('points', Float, ('number_of_points',
1437                                             'number_of_dimensions'))
1438    outfile.createVariable('elevation', Float, ('number_of_points',))
1439
1440    # Get handles to the variables
1441    points = outfile.variables['points']
1442    elevation = outfile.variables['elevation']
1443
1444    lenv = index2-index1+1
[4776]1445    # Store data
[2852]1446    global_index = 0
[4776]1447    # for i in range(nrows):
[2852]1448    for i in range(i1_0,thisi+1,1):
1449        if verbose and i%((nrows+10)/10)==0:
1450            print 'Processing row %d of %d' %(i, nrows)
1451
1452        lower_index = global_index
1453
1454        v = dem_elevation_r[i,index1:index2+1]
1455        no_NODATA = sum(v == NODATA_value)
1456        if no_NODATA > 0:
[4776]1457            newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA
[2852]1458        else:
[4776]1459            newcols = lenv # ncols_in_bounding_box
[2852]1460
1461        telev = zeros(newcols, Float)
1462        tpoints = zeros((newcols, 2), Float)
1463
1464        local_index = 0
1465
1466        y = (nrows-i-1)*cellsize + yllcorner
1467        #for j in range(ncols):
1468        for j in range(j1_0,index2+1,1):
1469
1470            x = j*cellsize + xllcorner
1471            if easting_min <= x <= easting_max and \
1472               northing_min <= y <= northing_max and \
1473               dem_elevation_r[i,j] <> NODATA_value:
1474                tpoints[local_index, :] = [x-easting_min,y-northing_min]
1475                telev[local_index] = dem_elevation_r[i, j]
1476                global_index += 1
1477                local_index += 1
1478
1479        upper_index = global_index
1480
1481        if upper_index == lower_index + newcols:
1482            points[lower_index:upper_index, :] = tpoints
1483            elevation[lower_index:upper_index] = telev
1484
1485    assert global_index == nopoints, 'index not equal to number of points'
1486
1487    infile.close()
1488    outfile.close()
1489
1490
1491
1492def _read_hecras_cross_sections(lines):
1493    """Return block of surface lines for each cross section
1494    Starts with SURFACE LINE,
1495    Ends with END CROSS-SECTION
1496    """
1497
1498    points = []
1499
1500    reading_surface = False
1501    for i, line in enumerate(lines):
1502
1503        if len(line.strip()) == 0:    #Ignore blanks
1504            continue
1505
1506        if lines[i].strip().startswith('SURFACE LINE'):
1507            reading_surface = True
1508            continue
1509
1510        if lines[i].strip().startswith('END') and reading_surface:
1511            yield points
1512            reading_surface = False
1513            points = []
1514
1515        if reading_surface:
1516            fields = line.strip().split(',')
1517            easting = float(fields[0])
1518            northing = float(fields[1])
1519            elevation = float(fields[2])
1520            points.append([easting, northing, elevation])
1521
1522
1523
1524
1525def hecras_cross_sections2pts(basename_in,
1526                              basename_out=None,
1527                              verbose=False):
1528    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
1529
1530    Example:
1531
1532
1533# RAS export file created on Mon 15Aug2005 11:42
1534# by HEC-RAS Version 3.1.1
1535
1536BEGIN HEADER:
1537  UNITS: METRIC
1538  DTM TYPE: TIN
1539  DTM: v:\1\cit\perth_topo\river_tin
1540  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
1541  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
1542  MAP PROJECTION: UTM
1543  PROJECTION ZONE: 50
1544  DATUM: AGD66
1545  VERTICAL DATUM:
1546  NUMBER OF REACHES:  19
1547  NUMBER OF CROSS-SECTIONS:  14206
1548END HEADER:
1549
1550
1551Only the SURFACE LINE data of the following form will be utilised
1552
1553  CROSS-SECTION:
1554    STREAM ID:Southern-Wungong
1555    REACH ID:Southern-Wungong
1556    STATION:19040.*
1557    CUT LINE:
1558      405548.671603161 , 6438142.7594925
1559      405734.536092045 , 6438326.10404912
1560      405745.130459356 , 6438331.48627354
1561      405813.89633823 , 6438368.6272789
1562    SURFACE LINE:
1563     405548.67,   6438142.76,   35.37
1564     405552.24,   6438146.28,   35.41
1565     405554.78,   6438148.78,   35.44
1566     405555.80,   6438149.79,   35.44
1567     405559.37,   6438153.31,   35.45
1568     405560.88,   6438154.81,   35.44
1569     405562.93,   6438156.83,   35.42
1570     405566.50,   6438160.35,   35.38
1571     405566.99,   6438160.83,   35.37
1572     ...
1573   END CROSS-SECTION
1574
1575    Convert to NetCDF pts format which is
1576
1577    points:  (Nx2) Float array
1578    elevation: N Float array
1579    """
1580
1581    import os
1582    from Scientific.IO.NetCDF import NetCDFFile
1583    from Numeric import Float, zeros, reshape
1584
1585    root = basename_in
1586
1587    #Get ASCII file
1588    infile = open(root + '.sdf', 'r')  #Open SDF file for read
1589
1590    if verbose: print 'Reading DEM from %s' %(root + '.sdf')
1591
1592    lines = infile.readlines()
1593    infile.close()
1594
1595    if verbose: print 'Converting to pts format'
1596
1597    i = 0
1598    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
1599        i += 1
1600
1601    assert lines[i].strip().upper() == 'BEGIN HEADER:'
1602    i += 1
1603
1604    assert lines[i].strip().upper().startswith('UNITS:')
1605    units = lines[i].strip().split()[1]
1606    i += 1
1607
1608    assert lines[i].strip().upper().startswith('DTM TYPE:')
1609    i += 1
1610
1611    assert lines[i].strip().upper().startswith('DTM:')
1612    i += 1
1613
1614    assert lines[i].strip().upper().startswith('STREAM')
1615    i += 1
1616
1617    assert lines[i].strip().upper().startswith('CROSS')
1618    i += 1
1619
1620    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
1621    projection = lines[i].strip().split(':')[1]
1622    i += 1
1623
1624    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
1625    zone = int(lines[i].strip().split(':')[1])
1626    i += 1
1627
1628    assert lines[i].strip().upper().startswith('DATUM:')
1629    datum = lines[i].strip().split(':')[1]
1630    i += 1
1631
1632    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
1633    i += 1
1634
1635    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
1636    i += 1
1637
1638    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
1639    number_of_cross_sections = int(lines[i].strip().split(':')[1])
1640    i += 1
1641
1642
1643    #Now read all points
1644    points = []
1645    elevation = []
1646    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
1647        for k, entry in enumerate(entries):
1648            points.append(entry[:2])
1649            elevation.append(entry[2])
1650
1651
1652    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
1653          %(j+1, number_of_cross_sections)
1654    assert j+1 == number_of_cross_sections, msg
1655
1656    #Get output file
1657    if basename_out == None:
1658        ptsname = root + '.pts'
1659    else:
1660        ptsname = basename_out + '.pts'
1661
[4455]1662    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
1663    geo = Geospatial_data(points, {"elevation":elevation},
1664                          verbose=verbose, geo_reference=geo_ref)
1665    geo.export_points_file(ptsname)
[2852]1666
[4462]1667def export_grid(basename_in, extra_name_out = None,
1668                quantities = None, # defaults to elevation
1669                timestep = None,
1670                reduction = None,
1671                cellsize = 10,
1672                NODATA_value = -9999,
1673                easting_min = None,
1674                easting_max = None,
1675                northing_min = None,
1676                northing_max = None,
1677                verbose = False,
1678                origin = None,
1679                datum = 'WGS84',
1680                format = 'ers'):
1681    """
1682   
1683    Wrapper for sww2dem. - see sww2dem to find out what most of the
1684    parameters do.
[2852]1685
[4462]1686    Quantities is a list of quantities.  Each quantity will be
1687    calculated for each sww file.
1688
1689    This returns the basenames of the files returned, which is made up
1690    of the dir and all of the file name, except the extension.
1691
1692    This function returns the names of the files produced.
[4535]1693
1694    It will also produce as many output files as there are input sww files.
[4462]1695    """
1696   
1697    if quantities is None:
1698        quantities = ['elevation']
1699       
1700    if type(quantities) is str:
1701            quantities = [quantities]
1702
1703    # How many sww files are there?
1704    dir, base = os.path.split(basename_in)
[4489]1705    #print "basename_in",basename_in
1706    #print "base",base
[4586]1707
1708    iterate_over = get_all_swwfiles(dir,base,verbose)
[4463]1709   
[4526]1710    if dir == "":
1711        dir = "." # Unix compatibility
[4462]1712   
1713    files_out = []
[4586]1714    #print 'sww_file',iterate_over
[4548]1715    for sww_file in iterate_over:
[4462]1716        for quantity in quantities:
1717            if extra_name_out is None:
1718                basename_out = sww_file + '_' + quantity
1719            else:
1720                basename_out = sww_file + '_' + quantity + '_' \
1721                               + extra_name_out
[4524]1722#            print "basename_out", basename_out
[4462]1723       
[4524]1724            file_out = sww2dem(dir+sep+sww_file, dir+sep+basename_out,
[4462]1725                               quantity, 
1726                               timestep,
1727                               reduction,
1728                               cellsize,
1729                               NODATA_value,
1730                               easting_min,
1731                               easting_max,
1732                               northing_min,
1733                               northing_max,
1734                               verbose,
1735                               origin,
1736                               datum,
1737                               format)
1738            files_out.append(file_out)
1739    #print "basenames_out after",basenames_out
1740    return files_out
[4545]1741
1742
1743def get_timeseries(production_dirs, output_dir, scenario_name, gauges_dir_name,
1744                   plot_quantity, generate_fig = False,
1745                   reportname = None, surface = False, time_min = None,
1746                   time_max = None, title_on = False, verbose = True,
1747                   nodes=None):
1748    """
1749    nodes - number of processes used.
1750
1751    warning - this function has no tests
1752    """
1753    if reportname == None:
1754        report = False
1755    else:
1756        report = True
1757       
1758    if nodes is None:
1759        is_parallel = False
1760    else:
1761        is_parallel = True
1762       
1763    # Generate figures
1764    swwfiles = {}
[4462]1765   
[4545]1766    if is_parallel is True:   
1767        for i in range(nodes):
1768            print 'Sending node %d of %d' %(i,nodes)
1769            swwfiles = {}
1770            if not reportname == None:
1771                reportname = report_name + '_%s' %(i)
1772            for label_id in production_dirs.keys():
1773                if label_id == 'boundaries':
1774                    swwfile = best_boundary_sww
1775                else:
1776                    file_loc = output_dir + label_id + sep
1777                    sww_extra = '_P%s_%s' %(i,nodes)
1778                    swwfile = file_loc + scenario_name + sww_extra + '.sww'
1779                    print 'swwfile',swwfile
1780                    swwfiles[swwfile] = label_id
1781
1782                texname, elev_output = sww2timeseries(swwfiles,
1783                                              gauges_dir_name,
1784                                              production_dirs,
1785                                              report = report,
1786                                              reportname = reportname,
1787                                              plot_quantity = plot_quantity,
1788                                              generate_fig = generate_fig,
1789                                              surface = surface,
1790                                              time_min = time_min,
1791                                              time_max = time_max,
1792                                              title_on = title_on,
1793                                              verbose = verbose)
1794    else:   
1795        for label_id in production_dirs.keys():       
1796            if label_id == 'boundaries':
1797                print 'boundaries'
1798                file_loc = project.boundaries_in_dir
1799                swwfile = project.boundaries_dir_name3 + '.sww'
1800                #  swwfile = boundary_dir_filename
1801            else:
1802                file_loc = output_dir + label_id + sep
1803                swwfile = file_loc + scenario_name + '.sww'
1804            swwfiles[swwfile] = label_id
1805       
1806        texname, elev_output = sww2timeseries(swwfiles,
1807                                              gauges_dir_name,
1808                                              production_dirs,
1809                                              report = report,
1810                                              reportname = reportname,
1811                                              plot_quantity = plot_quantity,
1812                                              generate_fig = generate_fig,
1813                                              surface = surface,
1814                                              time_min = time_min,
1815                                              time_max = time_max,
1816                                              title_on = title_on,
1817                                              verbose = verbose)
1818                                         
1819
1820   
[2852]1821def sww2dem(basename_in, basename_out = None,
[4462]1822            quantity = None, # defaults to elevation
[2852]1823            timestep = None,
1824            reduction = None,
1825            cellsize = 10,
1826            NODATA_value = -9999,
1827            easting_min = None,
1828            easting_max = None,
1829            northing_min = None,
1830            northing_max = None,
1831            verbose = False,
1832            origin = None,
1833            datum = 'WGS84',
[4462]1834            format = 'ers'):
[2852]1835
[4663]1836    """Read SWW file and convert to Digitial Elevation model format
1837    (.asc or .ers)
[2852]1838
1839    Example (ASC):
1840
1841    ncols         3121
1842    nrows         1800
1843    xllcorner     722000
1844    yllcorner     5893000
1845    cellsize      25
1846    NODATA_value  -9999
1847    138.3698 137.4194 136.5062 135.5558 ..........
1848
1849    Also write accompanying file with same basename_in but extension .prj
1850    used to fix the UTM zone, datum, false northings and eastings.
1851
1852    The prj format is assumed to be as
1853
1854    Projection    UTM
1855    Zone          56
1856    Datum         WGS84
1857    Zunits        NO
1858    Units         METERS
1859    Spheroid      WGS84
1860    Xshift        0.0000000000
1861    Yshift        10000000.0000000000
1862    Parameters
1863
1864
1865    The parameter quantity must be the name of an existing quantity or
1866    an expression involving existing quantities. The default is
[4462]1867    'elevation'. Quantity is not a list of quantities.
[2852]1868
1869    if timestep (an index) is given, output quantity at that timestep
1870
1871    if reduction is given use that to reduce quantity over all timesteps.
1872
1873    datum
1874
1875    format can be either 'asc' or 'ers'
1876    """
1877
1878    import sys
[4663]1879    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, \
1880         sometrue
[2852]1881    from Numeric import array2string
1882
[4663]1883    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
1884         separate_points_by_polygon
1885    from anuga.abstract_2d_finite_volumes.util import \
1886         apply_expression_to_dictionary
[2852]1887
1888    msg = 'Format must be either asc or ers'
1889    assert format.lower() in ['asc', 'ers'], msg
1890
1891
1892    false_easting = 500000
1893    false_northing = 10000000
1894
1895    if quantity is None:
1896        quantity = 'elevation'
[4462]1897       
[2852]1898    if reduction is None:
1899        reduction = max
1900
1901    if basename_out is None:
1902        basename_out = basename_in + '_%s' %quantity
1903
[4462]1904    if quantity_formula.has_key(quantity):
1905        quantity = quantity_formula[quantity]
1906       
[2852]1907    swwfile = basename_in + '.sww'
1908    demfile = basename_out + '.' + format
1909    # Note the use of a .ers extension is optional (write_ermapper_grid will
1910    # deal with either option
[4413]1911   
[4548]1912    #if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors
[4413]1913   
[4551]1914    # Read sww file
[4524]1915    if verbose: 
1916        print 'Reading from %s' %swwfile
1917        print 'Output directory is %s' %basename_out
1918   
[2852]1919    from Scientific.IO.NetCDF import NetCDFFile
1920    fid = NetCDFFile(swwfile)
1921
1922    #Get extent and reference
1923    x = fid.variables['x'][:]
1924    y = fid.variables['y'][:]
1925    volumes = fid.variables['volumes'][:]
[3961]1926    if timestep is not None:
1927        times = fid.variables['time'][timestep]
[3967]1928    else:
1929        times = fid.variables['time'][:]
[2852]1930
1931    number_of_timesteps = fid.dimensions['number_of_timesteps']
1932    number_of_points = fid.dimensions['number_of_points']
[3967]1933   
[2852]1934    if origin is None:
1935
[4551]1936        # Get geo_reference
1937        # sww files don't have to have a geo_ref
[2852]1938        try:
1939            geo_reference = Geo_reference(NetCDFObject=fid)
1940        except AttributeError, e:
[4551]1941            geo_reference = Geo_reference() # Default georef object
[2852]1942
1943        xllcorner = geo_reference.get_xllcorner()
1944        yllcorner = geo_reference.get_yllcorner()
1945        zone = geo_reference.get_zone()
1946    else:
1947        zone = origin[0]
1948        xllcorner = origin[1]
1949        yllcorner = origin[2]
1950
1951
1952
[4663]1953    # FIXME: Refactor using code from Interpolation_function.statistics
1954    # (in interpolate.py)
[4551]1955    # Something like print swwstats(swwname)
[2852]1956    if verbose:
1957        print '------------------------------------------------'
1958        print 'Statistics of SWW file:'
1959        print '  Name: %s' %swwfile
1960        print '  Reference:'
1961        print '    Lower left corner: [%f, %f]'\
1962              %(xllcorner, yllcorner)
[3961]1963        if timestep is not None:
1964            print '    Time: %f' %(times)
1965        else:
1966            print '    Start time: %f' %fid.starttime[0]
[2852]1967        print '  Extent:'
1968        print '    x [m] in [%f, %f], len(x) == %d'\
1969              %(min(x.flat), max(x.flat), len(x.flat))
1970        print '    y [m] in [%f, %f], len(y) == %d'\
1971              %(min(y.flat), max(y.flat), len(y.flat))
[3961]1972        if timestep is not None:
1973            print '    t [s] = %f, len(t) == %d' %(times, 1)
1974        else:
1975            print '    t [s] in [%f, %f], len(t) == %d'\
[3967]1976              %(min(times), max(times), len(times))
[2852]1977        print '  Quantities [SI units]:'
[5189]1978        # Comment out for reduced memory consumption
[3967]1979        for name in ['stage', 'xmomentum', 'ymomentum']:
[2852]1980            q = fid.variables[name][:].flat
[3967]1981            if timestep is not None:
1982                q = q[timestep*len(x):(timestep+1)*len(x)]
1983            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1984        for name in ['elevation']:
1985            q = fid.variables[name][:].flat
1986            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1987           
[2891]1988    # Get quantity and reduce if applicable
[2852]1989    if verbose: print 'Processing quantity %s' %quantity
1990
[2891]1991    # Turn NetCDF objects into Numeric arrays
[5115]1992    try:
1993        q = fid.variables[quantity][:] 
[5189]1994       
1995       
[5115]1996    except:
1997        quantity_dict = {}
1998        for name in fid.variables.keys():
1999            quantity_dict[name] = fid.variables[name][:] 
2000        #Convert quantity expression to quantities found in sww file   
2001        q = apply_expression_to_dictionary(quantity, quantity_dict)
[5189]2002    #print "q.shape",q.shape
[2852]2003    if len(q.shape) == 2:
2004        #q has a time component and needs to be reduced along
2005        #the temporal dimension
2006        if verbose: print 'Reducing quantity %s' %quantity
2007        q_reduced = zeros( number_of_points, Float )
[3967]2008       
2009        if timestep is not None:
2010            for k in range(number_of_points):
2011                q_reduced[k] = q[timestep,k] 
2012        else:
2013            for k in range(number_of_points):
2014                q_reduced[k] = reduction( q[:,k] )
[2852]2015
2016        q = q_reduced
2017
2018    #Post condition: Now q has dimension: number_of_points
2019    assert len(q.shape) == 1
2020    assert q.shape[0] == number_of_points
2021
2022    if verbose:
2023        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
2024
2025
2026    #Create grid and update xll/yll corner and x,y
2027
2028    #Relative extent
2029    if easting_min is None:
2030        xmin = min(x)
2031    else:
2032        xmin = easting_min - xllcorner
2033
2034    if easting_max is None:
2035        xmax = max(x)
2036    else:
2037        xmax = easting_max - xllcorner
2038
2039    if northing_min is None:
2040        ymin = min(y)
2041    else:
2042        ymin = northing_min - yllcorner
2043
2044    if northing_max is None:
2045        ymax = max(y)
2046    else:
2047        ymax = northing_max - yllcorner
2048
2049
2050
2051    if verbose: print 'Creating grid'
2052    ncols = int((xmax-xmin)/cellsize)+1
2053    nrows = int((ymax-ymin)/cellsize)+1
2054
2055
2056    #New absolute reference and coordinates
2057    newxllcorner = xmin+xllcorner
2058    newyllcorner = ymin+yllcorner
2059
2060    x = x+xllcorner-newxllcorner
2061    y = y+yllcorner-newyllcorner
[3961]2062   
[2852]2063    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2064    assert len(vertex_points.shape) == 2
2065
2066    grid_points = zeros ( (ncols*nrows, 2), Float )
2067
2068
2069    for i in xrange(nrows):
2070        if format.lower() == 'asc':
2071            yg = i*cellsize
2072        else:
2073        #this will flip the order of the y values for ers
2074            yg = (nrows-i)*cellsize
2075
2076        for j in xrange(ncols):
2077            xg = j*cellsize
2078            k = i*ncols + j
2079
2080            grid_points[k,0] = xg
2081            grid_points[k,1] = yg
2082
2083    #Interpolate
[3514]2084    from anuga.fit_interpolate.interpolate import Interpolate
[2852]2085
[4480]2086    # Remove loners from vertex_points, volumes here
2087    vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
[4497]2088    #export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
[4522]2089    #import sys; sys.exit()
[2852]2090    interp = Interpolate(vertex_points, volumes, verbose = verbose)
2091
2092    #Interpolate using quantity values
2093    if verbose: print 'Interpolating'
2094    grid_values = interp.interpolate(q, grid_points).flat
2095
2096
2097    if verbose:
2098        print 'Interpolated values are in [%f, %f]' %(min(grid_values),
2099                                                      max(grid_values))
2100
2101    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2102    P = interp.mesh.get_boundary_polygon()
2103    outside_indices = outside_polygon(grid_points, P, closed=True)
2104
2105    for i in outside_indices:
2106        grid_values[i] = NODATA_value
2107
2108
2109
2110
2111    if format.lower() == 'ers':
2112        # setup ERS header information
2113        grid_values = reshape(grid_values,(nrows, ncols))
2114        header = {}
2115        header['datum'] = '"' + datum + '"'
2116        # FIXME The use of hardwired UTM and zone number needs to be made optional
2117        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
2118        header['projection'] = '"UTM-' + str(zone) + '"'
2119        header['coordinatetype'] = 'EN'
2120        if header['coordinatetype'] == 'LL':
2121            header['longitude'] = str(newxllcorner)
2122            header['latitude'] = str(newyllcorner)
2123        elif header['coordinatetype'] == 'EN':
2124            header['eastings'] = str(newxllcorner)
2125            header['northings'] = str(newyllcorner)
2126        header['nullcellvalue'] = str(NODATA_value)
2127        header['xdimension'] = str(cellsize)
2128        header['ydimension'] = str(cellsize)
2129        header['value'] = '"' + quantity + '"'
2130        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
2131
2132
2133        #Write
2134        if verbose: print 'Writing %s' %demfile
2135        import ermapper_grids
2136        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
2137
2138        fid.close()
2139    else:
2140        #Write to Ascii format
2141
2142        #Write prj file
2143        prjfile = basename_out + '.prj'
2144
2145        if verbose: print 'Writing %s' %prjfile
2146        prjid = open(prjfile, 'w')
2147        prjid.write('Projection    %s\n' %'UTM')
2148        prjid.write('Zone          %d\n' %zone)
2149        prjid.write('Datum         %s\n' %datum)
2150        prjid.write('Zunits        NO\n')
2151        prjid.write('Units         METERS\n')
2152        prjid.write('Spheroid      %s\n' %datum)
2153        prjid.write('Xshift        %d\n' %false_easting)
2154        prjid.write('Yshift        %d\n' %false_northing)
2155        prjid.write('Parameters\n')
2156        prjid.close()
2157
2158
2159
2160        if verbose: print 'Writing %s' %demfile
2161
2162        ascid = open(demfile, 'w')
2163
2164        ascid.write('ncols         %d\n' %ncols)
2165        ascid.write('nrows         %d\n' %nrows)
2166        ascid.write('xllcorner     %d\n' %newxllcorner)
2167        ascid.write('yllcorner     %d\n' %newyllcorner)
2168        ascid.write('cellsize      %f\n' %cellsize)
2169        ascid.write('NODATA_value  %d\n' %NODATA_value)
2170
2171
2172        #Get bounding polygon from mesh
2173        #P = interp.mesh.get_boundary_polygon()
2174        #inside_indices = inside_polygon(grid_points, P)
2175
2176        for i in range(nrows):
2177            if verbose and i%((nrows+10)/10)==0:
2178                print 'Doing row %d of %d' %(i, nrows)
2179
2180            base_index = (nrows-i-1)*ncols
2181
2182            slice = grid_values[base_index:base_index+ncols]
2183            s = array2string(slice, max_line_width=sys.maxint)
2184            ascid.write(s[1:-1] + '\n')
2185
2186
2187            #print
2188            #for j in range(ncols):
2189            #    index = base_index+j#
2190            #    print grid_values[index],
2191            #    ascid.write('%f ' %grid_values[index])
2192            #ascid.write('\n')
2193
2194        #Close
2195        ascid.close()
2196        fid.close()
[4462]2197        return basename_out
2198
[5189]2199
[2852]2200#Backwards compatibility
2201def sww2asc(basename_in, basename_out = None,
2202            quantity = None,
2203            timestep = None,
2204            reduction = None,
2205            cellsize = 10,
2206            verbose = False,
2207            origin = None):
2208    print 'sww2asc will soon be obsoleted - please use sww2dem'
2209    sww2dem(basename_in,
2210            basename_out = basename_out,
2211            quantity = quantity,
2212            timestep = timestep,
2213            reduction = reduction,
2214            cellsize = cellsize,
2215            verbose = verbose,
2216            origin = origin,
2217        datum = 'WGS84',
2218        format = 'asc')
2219
2220def sww2ers(basename_in, basename_out = None,
2221            quantity = None,
2222            timestep = None,
2223            reduction = None,
2224            cellsize = 10,
2225            verbose = False,
2226            origin = None,
2227            datum = 'WGS84'):
2228    print 'sww2ers will soon be obsoleted - please use sww2dem'
2229    sww2dem(basename_in,
2230            basename_out = basename_out,
2231            quantity = quantity,
2232            timestep = timestep,
2233            reduction = reduction,
2234            cellsize = cellsize,
2235            verbose = verbose,
2236            origin = origin,
[2891]2237            datum = datum,
2238            format = 'ers')
[2852]2239################################# END COMPATIBILITY ##############
2240
2241
2242
[2891]2243def sww2pts(basename_in, basename_out=None,
2244            data_points=None,
2245            quantity=None,
2246            timestep=None,
2247            reduction=None,
2248            NODATA_value=-9999,
2249            verbose=False,
2250            origin=None):
2251            #datum = 'WGS84')
2252
2253
2254    """Read SWW file and convert to interpolated values at selected points
2255
2256    The parameter quantity' must be the name of an existing quantity or
2257    an expression involving existing quantities. The default is
2258    'elevation'.
2259
2260    if timestep (an index) is given, output quantity at that timestep
2261
2262    if reduction is given use that to reduce quantity over all timesteps.
2263
2264    data_points (Nx2 array) give locations of points where quantity is to be computed.
2265   
2266    """
2267
2268    import sys
2269    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
2270    from Numeric import array2string
2271
[3514]2272    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
[3560]2273    from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
[2891]2274
[3514]2275    from anuga.geospatial_data.geospatial_data import Geospatial_data
[2891]2276
2277    if quantity is None:
2278        quantity = 'elevation'
2279
2280    if reduction is None:
2281        reduction = max
2282
2283    if basename_out is None:
2284        basename_out = basename_in + '_%s' %quantity
2285
2286    swwfile = basename_in + '.sww'
2287    ptsfile = basename_out + '.pts'
2288
2289    # Read sww file
2290    if verbose: print 'Reading from %s' %swwfile
2291    from Scientific.IO.NetCDF import NetCDFFile
2292    fid = NetCDFFile(swwfile)
2293
2294    # Get extent and reference
2295    x = fid.variables['x'][:]
2296    y = fid.variables['y'][:]
2297    volumes = fid.variables['volumes'][:]
2298
2299    number_of_timesteps = fid.dimensions['number_of_timesteps']
2300    number_of_points = fid.dimensions['number_of_points']
2301    if origin is None:
2302
2303        # Get geo_reference
2304        # sww files don't have to have a geo_ref
2305        try:
2306            geo_reference = Geo_reference(NetCDFObject=fid)
2307        except AttributeError, e:
2308            geo_reference = Geo_reference() #Default georef object
2309
2310        xllcorner = geo_reference.get_xllcorner()
2311        yllcorner = geo_reference.get_yllcorner()
2312        zone = geo_reference.get_zone()
2313    else:
2314        zone = origin[0]
2315        xllcorner = origin[1]
2316        yllcorner = origin[2]
2317
2318
2319
2320    # FIXME: Refactor using code from file_function.statistics
2321    # Something like print swwstats(swwname)
2322    if verbose:
2323        x = fid.variables['x'][:]
2324        y = fid.variables['y'][:]
2325        times = fid.variables['time'][:]
2326        print '------------------------------------------------'
2327        print 'Statistics of SWW file:'
2328        print '  Name: %s' %swwfile
2329        print '  Reference:'
2330        print '    Lower left corner: [%f, %f]'\
2331              %(xllcorner, yllcorner)
2332        print '    Start time: %f' %fid.starttime[0]
2333        print '  Extent:'
2334        print '    x [m] in [%f, %f], len(x) == %d'\
2335              %(min(x.flat), max(x.flat), len(x.flat))
2336        print '    y [m] in [%f, %f], len(y) == %d'\
2337              %(min(y.flat), max(y.flat), len(y.flat))
2338        print '    t [s] in [%f, %f], len(t) == %d'\
2339              %(min(times), max(times), len(times))
2340        print '  Quantities [SI units]:'
2341        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2342            q = fid.variables[name][:].flat
2343            print '    %s in [%f, %f]' %(name, min(q), max(q))
2344
2345
2346
2347    # Get quantity and reduce if applicable
2348    if verbose: print 'Processing quantity %s' %quantity
2349
2350    # Turn NetCDF objects into Numeric arrays
2351    quantity_dict = {}
2352    for name in fid.variables.keys():
2353        quantity_dict[name] = fid.variables[name][:]
2354
2355
2356
2357    # Convert quantity expression to quantities found in sww file   
2358    q = apply_expression_to_dictionary(quantity, quantity_dict)
2359
2360
2361
2362    if len(q.shape) == 2:
2363        # q has a time component and needs to be reduced along
2364        # the temporal dimension
2365        if verbose: print 'Reducing quantity %s' %quantity
2366        q_reduced = zeros( number_of_points, Float )
2367
2368        for k in range(number_of_points):
2369            q_reduced[k] = reduction( q[:,k] )
2370
2371        q = q_reduced
2372
2373    # Post condition: Now q has dimension: number_of_points
2374    assert len(q.shape) == 1
2375    assert q.shape[0] == number_of_points
2376
2377
2378    if verbose:
2379        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
2380
2381
2382    # Create grid and update xll/yll corner and x,y
2383    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2384    assert len(vertex_points.shape) == 2
2385
2386    # Interpolate
[3514]2387    from anuga.fit_interpolate.interpolate import Interpolate
[2891]2388    interp = Interpolate(vertex_points, volumes, verbose = verbose)
2389
2390    # Interpolate using quantity values
2391    if verbose: print 'Interpolating'
2392    interpolated_values = interp.interpolate(q, data_points).flat
2393
2394
2395    if verbose:
2396        print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
2397                                                      max(interpolated_values))
2398
2399    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2400    P = interp.mesh.get_boundary_polygon()
2401    outside_indices = outside_polygon(data_points, P, closed=True)
2402
2403    for i in outside_indices:
2404        interpolated_values[i] = NODATA_value
2405
2406
2407    # Store results   
2408    G = Geospatial_data(data_points=data_points,
2409                        attributes=interpolated_values)
2410
2411    G.export_points_file(ptsfile, absolute = True)
2412
[2931]2413    fid.close()
[2891]2414
2415
[2852]2416def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2417                                  use_cache = False,
2418                                  verbose = False):
2419    """Read Digitial Elevation model from the following ASCII format (.asc)
2420
2421    Example:
2422
2423    ncols         3121
2424    nrows         1800
2425    xllcorner     722000
2426    yllcorner     5893000
2427    cellsize      25
2428    NODATA_value  -9999
2429    138.3698 137.4194 136.5062 135.5558 ..........
2430
2431    Convert basename_in + '.asc' to NetCDF format (.dem)
2432    mimicking the ASCII format closely.
2433
2434
2435    An accompanying file with same basename_in but extension .prj must exist
2436    and is used to fix the UTM zone, datum, false northings and eastings.
2437
2438    The prj format is assumed to be as
2439
2440    Projection    UTM
2441    Zone          56
2442    Datum         WGS84
2443    Zunits        NO
2444    Units         METERS
2445    Spheroid      WGS84
2446    Xshift        0.0000000000
2447    Yshift        10000000.0000000000
2448    Parameters
2449    """
2450
2451
2452
2453    kwargs = {'basename_out': basename_out, 'verbose': verbose}
2454
2455    if use_cache is True:
2456        from caching import cache
2457        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
[4688]2458                       dependencies = [basename_in + '.asc',
2459                                       basename_in + '.prj'],
[2852]2460                       verbose = verbose)
2461
2462    else:
2463        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
2464
2465    return result
2466
2467
2468
2469
2470
2471
2472def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2473                                  verbose = False):
2474    """Read Digitial Elevation model from the following ASCII format (.asc)
2475
2476    Internal function. See public function convert_dem_from_ascii2netcdf for details.
2477    """
2478
2479    import os
2480    from Scientific.IO.NetCDF import NetCDFFile
2481    from Numeric import Float, array
2482
2483    #root, ext = os.path.splitext(basename_in)
2484    root = basename_in
2485
2486    ###########################################
2487    # Read Meta data
2488    if verbose: print 'Reading METADATA from %s' %root + '.prj'
2489    metadatafile = open(root + '.prj')
2490    metalines = metadatafile.readlines()
2491    metadatafile.close()
2492
2493    L = metalines[0].strip().split()
2494    assert L[0].strip().lower() == 'projection'
2495    projection = L[1].strip()                   #TEXT
2496
2497    L = metalines[1].strip().split()
2498    assert L[0].strip().lower() == 'zone'
2499    zone = int(L[1].strip())
2500
2501    L = metalines[2].strip().split()
2502    assert L[0].strip().lower() == 'datum'
2503    datum = L[1].strip()                        #TEXT
2504
2505    L = metalines[3].strip().split()
2506    assert L[0].strip().lower() == 'zunits'     #IGNORE
2507    zunits = L[1].strip()                       #TEXT
2508
2509    L = metalines[4].strip().split()
2510    assert L[0].strip().lower() == 'units'
2511    units = L[1].strip()                        #TEXT
2512
2513    L = metalines[5].strip().split()
2514    assert L[0].strip().lower() == 'spheroid'   #IGNORE
2515    spheroid = L[1].strip()                     #TEXT
2516
2517    L = metalines[6].strip().split()
2518    assert L[0].strip().lower() == 'xshift'
2519    false_easting = float(L[1].strip())
2520
2521    L = metalines[7].strip().split()
2522    assert L[0].strip().lower() == 'yshift'
2523    false_northing = float(L[1].strip())
2524
2525    #print false_easting, false_northing, zone, datum
2526
2527
2528    ###########################################
2529    #Read DEM data
2530
2531    datafile = open(basename_in + '.asc')
2532
2533    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
2534    lines = datafile.readlines()
2535    datafile.close()
2536
2537    if verbose: print 'Got', len(lines), ' lines'
2538
2539    ncols = int(lines[0].split()[1].strip())
2540    nrows = int(lines[1].split()[1].strip())
[4824]2541
2542    # Do cellsize (line 4) before line 2 and 3
2543    cellsize = float(lines[4].split()[1].strip())   
2544
2545    # Checks suggested by Joaquim Luis
2546    # Our internal representation of xllcorner
2547    # and yllcorner is non-standard.
2548    xref = lines[2].split()
2549    if xref[0].strip() == 'xllcorner':
2550        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
2551    elif xref[0].strip() == 'xllcenter':
2552        xllcorner = float(xref[1].strip())
2553    else:
2554        msg = 'Unknown keyword: %s' %xref[0].strip()
2555        raise Exception, msg
2556
2557    yref = lines[3].split()
2558    if yref[0].strip() == 'yllcorner':
2559        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
2560    elif yref[0].strip() == 'yllcenter':
2561        yllcorner = float(yref[1].strip())
2562    else:
2563        msg = 'Unknown keyword: %s' %yref[0].strip()
2564        raise Exception, msg
2565       
2566
[2852]2567    NODATA_value = int(lines[5].split()[1].strip())
2568
2569    assert len(lines) == nrows + 6
2570
2571
2572    ##########################################
2573
2574
2575    if basename_out == None:
2576        netcdfname = root + '.dem'
2577    else:
2578        netcdfname = basename_out + '.dem'
2579
2580    if verbose: print 'Store to NetCDF file %s' %netcdfname
2581    # NetCDF file definition
2582    fid = NetCDFFile(netcdfname, 'w')
2583
2584    #Create new file
2585    fid.institution = 'Geoscience Australia'
2586    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
2587                      'of spatial point data'
2588
2589    fid.ncols = ncols
2590    fid.nrows = nrows
2591    fid.xllcorner = xllcorner
2592    fid.yllcorner = yllcorner
2593    fid.cellsize = cellsize
2594    fid.NODATA_value = NODATA_value
2595
2596    fid.zone = zone
2597    fid.false_easting = false_easting
2598    fid.false_northing = false_northing
2599    fid.projection = projection
2600    fid.datum = datum
2601    fid.units = units
2602
2603
2604    # dimension definitions
2605    fid.createDimension('number_of_rows', nrows)
2606    fid.createDimension('number_of_columns', ncols)
2607
2608    # variable definitions
2609    fid.createVariable('elevation', Float, ('number_of_rows',
2610                                            'number_of_columns'))
2611
2612    # Get handles to the variables
2613    elevation = fid.variables['elevation']
2614
2615    #Store data
2616    n = len(lines[6:])
2617    for i, line in enumerate(lines[6:]):
2618        fields = line.split()
2619        if verbose and i%((n+10)/10)==0:
2620            print 'Processing row %d of %d' %(i, nrows)
2621
2622        elevation[i, :] = array([float(x) for x in fields])
2623
2624    fid.close()
2625
2626
2627
2628
2629
2630def ferret2sww(basename_in, basename_out = None,
2631               verbose = False,
2632               minlat = None, maxlat = None,
2633               minlon = None, maxlon = None,
2634               mint = None, maxt = None, mean_stage = 0,
2635               origin = None, zscale = 1,
2636               fail_on_NaN = True,
2637               NaN_filler = 0,
2638               elevation = None,
[3694]2639               inverted_bathymetry = True
[2852]2640               ): #FIXME: Bathymetry should be obtained
2641                                  #from MOST somehow.
2642                                  #Alternatively from elsewhere
2643                                  #or, as a last resort,
2644                                  #specified here.
2645                                  #The value of -100 will work
2646                                  #for the Wollongong tsunami
2647                                  #scenario but is very hacky
2648    """Convert MOST and 'Ferret' NetCDF format for wave propagation to
[3560]2649    sww format native to abstract_2d_finite_volumes.
[2852]2650
2651    Specify only basename_in and read files of the form
2652    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
2653    relative height, x-velocity and y-velocity, respectively.
2654
2655    Also convert latitude and longitude to UTM. All coordinates are
2656    assumed to be given in the GDA94 datum.
2657
2658    min's and max's: If omitted - full extend is used.
2659    To include a value min may equal it, while max must exceed it.
2660    Lat and lon are assuemd to be in decimal degrees
2661
2662    origin is a 3-tuple with geo referenced
2663    UTM coordinates (zone, easting, northing)
2664
2665    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
2666    which means that longitude is the fastest
2667    varying dimension (row major order, so to speak)
2668
2669    ferret2sww uses grid points as vertices in a triangular grid
2670    counting vertices from lower left corner upwards, then right
2671    """
2672
2673    import os
2674    from Scientific.IO.NetCDF import NetCDFFile
2675    from Numeric import Float, Int, Int32, searchsorted, zeros, array
2676    from Numeric import allclose, around
2677
2678    precision = Float
2679
2680    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
2681
2682    if minlat != None:
2683        assert -90 < minlat < 90 , msg
2684    if maxlat != None:
2685        assert -90 < maxlat < 90 , msg
[4050]2686        if minlat != None:
2687            assert maxlat > minlat
[2852]2688    if minlon != None:
2689        assert -180 < minlon < 180 , msg
2690    if maxlon != None:
2691        assert -180 < maxlon < 180 , msg
[4050]2692        if minlon != None:
2693            assert maxlon > minlon
2694       
[2852]2695
2696
2697    #Get NetCDF data
2698    if verbose: print 'Reading files %s_*.nc' %basename_in
[3694]2699    #print "basename_in + '_ha.nc'",basename_in + '_ha.nc'
[2852]2700    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
2701    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
2702    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
2703    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
2704
2705    if basename_out is None:
2706        swwname = basename_in + '.sww'
2707    else:
2708        swwname = basename_out + '.sww'
2709
[4418]2710    # Get dimensions of file_h
[2852]2711    for dimension in file_h.dimensions.keys():
2712        if dimension[:3] == 'LON':
2713            dim_h_longitude = dimension
2714        if dimension[:3] == 'LAT':
2715            dim_h_latitude = dimension
2716        if dimension[:4] == 'TIME':
2717            dim_h_time = dimension
2718
2719#    print 'long:', dim_h_longitude
2720#    print 'lats:', dim_h_latitude
2721#    print 'times:', dim_h_time
2722
2723    times = file_h.variables[dim_h_time]
2724    latitudes = file_h.variables[dim_h_latitude]
2725    longitudes = file_h.variables[dim_h_longitude]
[5347]2726
2727    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
2728                                                  longitudes[:],
2729                                                  minlat, maxlat,
2730                                                  minlon, maxlon)
[4418]2731    # get dimensions for file_e
[2852]2732    for dimension in file_e.dimensions.keys():
2733        if dimension[:3] == 'LON':
2734            dim_e_longitude = dimension
2735        if dimension[:3] == 'LAT':
2736            dim_e_latitude = dimension
2737
[4418]2738    # get dimensions for file_u
[2852]2739    for dimension in file_u.dimensions.keys():
2740        if dimension[:3] == 'LON':
2741            dim_u_longitude = dimension
2742        if dimension[:3] == 'LAT':
2743            dim_u_latitude = dimension
2744        if dimension[:4] == 'TIME':
2745            dim_u_time = dimension
2746           
[4418]2747    # get dimensions for file_v
[2852]2748    for dimension in file_v.dimensions.keys():
2749        if dimension[:3] == 'LON':
2750            dim_v_longitude = dimension
2751        if dimension[:3] == 'LAT':
2752            dim_v_latitude = dimension
2753        if dimension[:4] == 'TIME':
2754            dim_v_time = dimension
2755
2756
[4418]2757    # Precision used by most for lat/lon is 4 or 5 decimals
[2852]2758    e_lat = around(file_e.variables[dim_e_latitude][:], 5)
2759    e_lon = around(file_e.variables[dim_e_longitude][:], 5)
2760
[4418]2761    # Check that files are compatible
[2852]2762    assert allclose(latitudes, file_u.variables[dim_u_latitude])
2763    assert allclose(latitudes, file_v.variables[dim_v_latitude])
2764    assert allclose(latitudes, e_lat)
2765
2766    assert allclose(longitudes, file_u.variables[dim_u_longitude])
2767    assert allclose(longitudes, file_v.variables[dim_v_longitude])
2768    assert allclose(longitudes, e_lon)
2769
[4418]2770    if mint is None:
[2852]2771        jmin = 0
[4418]2772        mint = times[0]       
[2852]2773    else:
2774        jmin = searchsorted(times, mint)
[4024]2775       
[4418]2776    if maxt is None:
2777        jmax = len(times)
2778        maxt = times[-1]
[2852]2779    else:
2780        jmax = searchsorted(times, maxt)
2781
[4037]2782    #print "latitudes[:]",latitudes[:]
2783    #print "longitudes[:]",longitudes [:]
[4024]2784    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
2785                                                  longitudes[:],
[4418]2786                                                  minlat, maxlat,
2787                                                  minlon, maxlon)
[2852]2788
2789
2790    times = times[jmin:jmax]
2791    latitudes = latitudes[kmin:kmax]
2792    longitudes = longitudes[lmin:lmax]
2793
[4037]2794    #print "latitudes[:]",latitudes[:]
2795    #print "longitudes[:]",longitudes [:]
[2852]2796
2797    if verbose: print 'cropping'
2798    zname = 'ELEVATION'
2799
2800    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
2801    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
2802    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
2803    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
2804
2805    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
2806    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
2807    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
2808    #        from Numeric import asarray
2809    #        elevations=elevations.tolist()
2810    #        elevations.reverse()
2811    #        elevations=asarray(elevations)
2812    #    else:
2813    #        from Numeric import asarray
2814    #        elevations=elevations.tolist()
2815    #        elevations.reverse()
2816    #        elevations=asarray(elevations)
2817    #        'print hmmm'
2818
2819
2820
2821    #Get missing values
2822    nan_ha = file_h.variables['HA'].missing_value[0]
2823    nan_ua = file_u.variables['UA'].missing_value[0]
2824    nan_va = file_v.variables['VA'].missing_value[0]
2825    if hasattr(file_e.variables[zname],'missing_value'):
2826        nan_e  = file_e.variables[zname].missing_value[0]
2827    else:
2828        nan_e = None
2829
2830    #Cleanup
2831    from Numeric import sometrue
2832
2833    missing = (amplitudes == nan_ha)
2834    if sometrue (missing):
2835        if fail_on_NaN:
2836            msg = 'NetCDFFile %s contains missing values'\
2837                  %(basename_in+'_ha.nc')
2838            raise DataMissingValuesError, msg
2839        else:
2840            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
2841
2842    missing = (uspeed == nan_ua)
2843    if sometrue (missing):
2844        if fail_on_NaN:
2845            msg = 'NetCDFFile %s contains missing values'\
2846                  %(basename_in+'_ua.nc')
2847            raise DataMissingValuesError, msg
2848        else:
2849            uspeed = uspeed*(missing==0) + missing*NaN_filler
2850
2851    missing = (vspeed == nan_va)
2852    if sometrue (missing):
2853        if fail_on_NaN:
2854            msg = 'NetCDFFile %s contains missing values'\
2855                  %(basename_in+'_va.nc')
2856            raise DataMissingValuesError, msg
2857        else:
2858            vspeed = vspeed*(missing==0) + missing*NaN_filler
2859
2860
2861    missing = (elevations == nan_e)
2862    if sometrue (missing):
2863        if fail_on_NaN:
2864            msg = 'NetCDFFile %s contains missing values'\
2865                  %(basename_in+'_e.nc')
2866            raise DataMissingValuesError, msg
2867        else:
2868            elevations = elevations*(missing==0) + missing*NaN_filler
2869
2870    #######
2871
2872
2873
2874    number_of_times = times.shape[0]
2875    number_of_latitudes = latitudes.shape[0]
2876    number_of_longitudes = longitudes.shape[0]
2877
2878    assert amplitudes.shape[0] == number_of_times
2879    assert amplitudes.shape[1] == number_of_latitudes
2880    assert amplitudes.shape[2] == number_of_longitudes
2881
2882    if verbose:
2883        print '------------------------------------------------'
2884        print 'Statistics:'
2885        print '  Extent (lat/lon):'
2886        print '    lat in [%f, %f], len(lat) == %d'\
2887              %(min(latitudes.flat), max(latitudes.flat),
2888                len(latitudes.flat))
2889        print '    lon in [%f, %f], len(lon) == %d'\
2890              %(min(longitudes.flat), max(longitudes.flat),
2891                len(longitudes.flat))
2892        print '    t in [%f, %f], len(t) == %d'\
2893              %(min(times.flat), max(times.flat), len(times.flat))
2894
2895        q = amplitudes.flat
2896        name = 'Amplitudes (ha) [cm]'
2897        print %s in [%f, %f]' %(name, min(q), max(q))
2898
2899        q = uspeed.flat
2900        name = 'Speeds (ua) [cm/s]'
2901        print %s in [%f, %f]' %(name, min(q), max(q))
2902
2903        q = vspeed.flat
2904        name = 'Speeds (va) [cm/s]'
2905        print %s in [%f, %f]' %(name, min(q), max(q))
2906
2907        q = elevations.flat
2908        name = 'Elevations (e) [m]'
2909        print %s in [%f, %f]' %(name, min(q), max(q))
2910
2911
[4704]2912    # print number_of_latitudes, number_of_longitudes
[2852]2913    number_of_points = number_of_latitudes*number_of_longitudes
2914    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
2915
2916
2917    file_h.close()
2918    file_u.close()
2919    file_v.close()
2920    file_e.close()
2921
2922
2923    # NetCDF file definition
2924    outfile = NetCDFFile(swwname, 'w')
2925
[4387]2926    description = 'Converted from Ferret files: %s, %s, %s, %s'\
2927                  %(basename_in + '_ha.nc',
2928                    basename_in + '_ua.nc',
2929                    basename_in + '_va.nc',
2930                    basename_in + '_e.nc')
[4388]2931   
[4704]2932    # Create new file
[4416]2933    starttime = times[0]
[4862]2934   
[4455]2935    sww = Write_sww()
[4704]2936    sww.store_header(outfile, times, number_of_volumes,
2937                     number_of_points, description=description,
[4862]2938                     verbose=verbose,sww_precision=Float)
[2852]2939
[4704]2940    # Store
[3514]2941    from anuga.coordinate_transforms.redfearn import redfearn
[2852]2942    x = zeros(number_of_points, Float)  #Easting
2943    y = zeros(number_of_points, Float)  #Northing
2944
2945
2946    if verbose: print 'Making triangular grid'
[4704]2947
2948    # Check zone boundaries
[2852]2949    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
2950
2951    vertices = {}
2952    i = 0
2953    for k, lat in enumerate(latitudes):       #Y direction
2954        for l, lon in enumerate(longitudes):  #X direction
2955
2956            vertices[l,k] = i
2957
2958            zone, easting, northing = redfearn(lat,lon)
2959
2960            msg = 'Zone boundary crossed at longitude =', lon
2961            #assert zone == refzone, msg
2962            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
2963            x[i] = easting
2964            y[i] = northing
2965            i += 1
2966
2967    #Construct 2 triangles per 'rectangular' element
2968    volumes = []
2969    for l in range(number_of_longitudes-1):    #X direction
2970        for k in range(number_of_latitudes-1): #Y direction
2971            v1 = vertices[l,k+1]
2972            v2 = vertices[l,k]
2973            v3 = vertices[l+1,k+1]
2974            v4 = vertices[l+1,k]
2975
2976            volumes.append([v1,v2,v3]) #Upper element
2977            volumes.append([v4,v3,v2]) #Lower element
2978
2979    volumes = array(volumes)
2980
[4387]2981    if origin is None:
2982        origin = Geo_reference(refzone,min(x),min(y))
2983    geo_ref = write_NetCDF_georeference(origin, outfile)
2984   
[2852]2985    if elevation is not None:
2986        z = elevation
2987    else:
2988        if inverted_bathymetry:
2989            z = -1*elevations
2990        else:
2991            z = elevations
2992    #FIXME: z should be obtained from MOST and passed in here
2993
[4862]2994    #FIXME use the Write_sww instance(sww) to write this info
[2852]2995    from Numeric import resize
2996    z = resize(z,outfile.variables['z'][:].shape)
[4387]2997    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
2998    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
[3954]2999    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
[2852]3000    outfile.variables['elevation'][:] = z
[3954]3001    outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
[2852]3002
3003
3004
3005    #Time stepping
3006    stage = outfile.variables['stage']
3007    xmomentum = outfile.variables['xmomentum']
3008    ymomentum = outfile.variables['ymomentum']
3009
3010    if verbose: print 'Converting quantities'
3011    n = len(times)
3012    for j in range(n):
3013        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
3014        i = 0
3015        for k in range(number_of_latitudes):      #Y direction
3016            for l in range(number_of_longitudes): #X direction
3017                w = zscale*amplitudes[j,k,l]/100 + mean_stage
3018                stage[j,i] = w
3019                h = w - z[i]
3020                xmomentum[j,i] = uspeed[j,k,l]/100*h
3021                ymomentum[j,i] = vspeed[j,k,l]/100*h
3022                i += 1
3023
3024    #outfile.close()
3025
3026    #FIXME: Refactor using code from file_function.statistics
3027    #Something like print swwstats(swwname)
3028    if verbose:
3029        x = outfile.variables['x'][:]
3030        y = outfile.variables['y'][:]
3031        print '------------------------------------------------'
3032        print 'Statistics of output file:'
3033        print '  Name: %s' %swwname
3034        print '  Reference:'
3035        print '    Lower left corner: [%f, %f]'\
[4387]3036              %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
[4416]3037        print ' Start time: %f' %starttime
[4418]3038        print '    Min time: %f' %mint
3039        print '    Max time: %f' %maxt
[2852]3040        print '  Extent:'
3041        print '    x [m] in [%f, %f], len(x) == %d'\
3042              %(min(x.flat), max(x.flat), len(x.flat))
3043        print '    y [m] in [%f, %f], len(y) == %d'\
3044              %(min(y.flat), max(y.flat), len(y.flat))
3045        print '    t [s] in [%f, %f], len(t) == %d'\
3046              %(min(times), max(times), len(times))
3047        print '  Quantities [SI units]:'
3048        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
3049            q = outfile.variables[name][:].flat
3050            print '    %s in [%f, %f]' %(name, min(q), max(q))
3051
3052
3053
3054    outfile.close()
3055
3056
3057
3058
3059
[4303]3060def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
[2852]3061    """Template for converting typical text files with time series to
3062    NetCDF tms file.
3063
3064
3065    The file format is assumed to be either two fields separated by a comma:
3066
3067        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
3068
3069    E.g
3070
3071      31/08/04 00:00:00, 1.328223 0 0
3072      31/08/04 00:15:00, 1.292912 0 0
3073
[4303]3074    or time (seconds), value0 value1 value2 ...
3075   
3076      0.0, 1.328223 0 0
3077      0.1, 1.292912 0 0
3078
[2852]3079    will provide a time dependent function f(t) with three attributes
3080
3081    filename is assumed to be the rootname with extenisons .txt and .sww
3082    """
3083
3084    import time, calendar
3085    from Numeric import array
[3514]3086    from anuga.config import time_format
3087    from anuga.utilities.numerical_tools import ensure_numeric
[2852]3088
3089
3090
3091    fid = open(filename + '.txt')
3092    line = fid.readline()
3093    fid.close()
3094
3095    fields = line.split(',')
3096    msg = 'File %s must have the format date, value0 value1 value2 ...'
3097    assert len(fields) == 2, msg
3098
[4303]3099    if not time_as_seconds:
3100        try:
3101            starttime = calendar.timegm(time.strptime(fields[0], time_format))
3102        except ValueError:
3103            msg = 'First field in file %s must be' %filename
3104            msg += ' date-time with format %s.\n' %time_format
3105            msg += 'I got %s instead.' %fields[0]
3106            raise DataTimeError, msg
3107    else:
3108        try:
3109            starttime = float(fields[0])
3110        except Error:
3111            msg = "Bad time format"
3112            raise DataTimeError, msg
[2852]3113
3114
3115    #Split values
3116    values = []
3117    for value in fields[1].split():
3118        values.append(float(value))
3119
3120    q = ensure_numeric(values)
3121
3122    msg = 'ERROR: File must contain at least one independent value'
3123    assert len(q.shape) == 1, msg
3124
3125
3126
3127    #Read times proper
3128    from Numeric import zeros, Float, alltrue
[3514]3129    from anuga.config import time_format
[2852]3130    import time, calendar
3131
3132    fid = open(filename + '.txt')
3133    lines = fid.readlines()
3134    fid.close()
3135
3136    N = len(lines)
3137    d = len(q)
3138
3139    T = zeros(N, Float)       #Time
3140    Q = zeros((N, d), Float)  #Values
3141
3142    for i, line in enumerate(lines):
3143        fields = line.split(',')
[4303]3144        if not time_as_seconds:
3145            realtime = calendar.timegm(time.strptime(fields[0], time_format))
3146        else:
3147             realtime = float(fields[0])
[2852]3148        T[i] = realtime - starttime
3149
3150        for j, value in enumerate(fields[1].split()):
3151            Q[i, j] = float(value)
3152
3153    msg = 'File %s must list time as a monotonuosly ' %filename
3154    msg += 'increasing sequence'
3155    assert alltrue( T[1:] - T[:-1] > 0 ), msg
3156
3157    #Create NetCDF file
3158    from Scientific.IO.NetCDF import NetCDFFile
3159
3160    fid = NetCDFFile(filename + '.tms', 'w')
3161
3162
3163    fid.institution = 'Geoscience Australia'
3164    fid.description = 'Time series'
3165
3166
3167    #Reference point
3168    #Start time in seconds since the epoch (midnight 1/1/1970)
3169    #FIXME: Use Georef
3170    fid.starttime = starttime
3171
3172    # dimension definitions
3173    #fid.createDimension('number_of_volumes', self.number_of_volumes)
3174    #fid.createDimension('number_of_vertices', 3)
3175
3176
3177    fid.createDimension('number_of_timesteps', len(T))
3178
3179    fid.createVariable('time', Float, ('number_of_timesteps',))
3180
3181    fid.variables['time'][:] = T
3182
3183    for i in range(Q.shape[1]):
3184        try:
3185            name = quantity_names[i]
3186        except:
3187            name = 'Attribute%d'%i
3188
3189        fid.createVariable(name, Float, ('number_of_timesteps',))
3190        fid.variables[name][:] = Q[:,i]
3191
3192    fid.close()
3193
3194
3195def extent_sww(file_name):
3196    """
3197    Read in an sww file.
3198
3199    Input;
3200    file_name - the sww file
3201
3202    Output;
3203    z - Vector of bed elevation
3204    volumes - Array.  Each row has 3 values, representing
3205    the vertices that define the volume
3206    time - Vector of the times where there is stage information
3207    stage - array with respect to time and vertices (x,y)
3208    """
3209
3210
3211    from Scientific.IO.NetCDF import NetCDFFile
3212
3213    #Check contents
3214    #Get NetCDF
3215    fid = NetCDFFile(file_name, 'r')
3216
3217    # Get the variables
3218    x = fid.variables['x'][:]
3219    y = fid.variables['y'][:]
3220    stage = fid.variables['stage'][:]
3221    #print "stage",stage
3222    #print "stage.shap",stage.shape
3223    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
3224    #print "min(stage)",min(stage)
3225
3226    fid.close()
3227
3228    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
3229
3230
[5276]3231def sww2domain(filename, boundary=None, t=None,
3232               fail_if_NaN=True ,NaN_filler=0,
3233               verbose = False, very_verbose = False):
[2852]3234    """
3235    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
3236
3237    Boundary is not recommended if domain.smooth is not selected, as it
3238    uses unique coordinates, but not unique boundaries. This means that
3239    the boundary file will not be compatable with the coordinates, and will
3240    give a different final boundary, or crash.
3241    """
3242    NaN=9.969209968386869e+036
3243    #initialise NaN.
3244
3245    from Scientific.IO.NetCDF import NetCDFFile
3246    from shallow_water import Domain
3247    from Numeric import asarray, transpose, resize
3248
3249    if verbose: print 'Reading from ', filename
3250    fid = NetCDFFile(filename, 'r')    #Open existing file for read
3251    time = fid.variables['time']       #Timesteps
3252    if t is None:
3253        t = time[-1]
3254    time_interp = get_time_interp(time,t)
3255
3256    # Get the variables as Numeric arrays
3257    x = fid.variables['x'][:]             #x-coordinates of vertices
3258    y = fid.variables['y'][:]             #y-coordinates of vertices
3259    elevation = fid.variables['elevation']     #Elevation
3260    stage = fid.variables['stage']     #Water level
3261    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
3262    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
3263
3264    starttime = fid.starttime[0]
3265    volumes = fid.variables['volumes'][:] #Connectivity
[5276]3266    coordinates = transpose(asarray([x.tolist(),y.tolist()]))
3267    #FIXME (Ole): Something like this might be better: concatenate( (x, y), axis=1 )
3268    # or concatenate( (x[:,NewAxis],x[:,NewAxis]), axis=1 )
[2852]3269
3270    conserved_quantities = []
3271    interpolated_quantities = {}
3272    other_quantities = []
3273
3274    # get geo_reference
3275    #sww files don't have to have a geo_ref
3276    try:
3277        geo_reference = Geo_reference(NetCDFObject=fid)
3278    except: #AttributeError, e:
3279        geo_reference = None
3280
3281    if verbose: print '    getting quantities'
3282    for quantity in fid.variables.keys():
3283        dimensions = fid.variables[quantity].dimensions
3284        if 'number_of_timesteps' in dimensions:
3285            conserved_quantities.append(quantity)
3286            interpolated_quantities[quantity]=\
3287                  interpolated_quantity(fid.variables[quantity][:],time_interp)
3288        else: other_quantities.append(quantity)
3289
3290    other_quantities.remove('x')
3291    other_quantities.remove('y')
3292    other_quantities.remove('z')
3293    other_quantities.remove('volumes')
[4455]3294    try:
3295        other_quantities.remove('stage_range')
3296        other_quantities.remove('xmomentum_range')
3297        other_quantities.remove('ymomentum_range')
3298        other_quantities.remove('elevation_range')
3299    except:
3300        pass
3301       
[2852]3302
3303    conserved_quantities.remove('time')
3304
3305    if verbose: print '    building domain'
3306    #    From domain.Domain:
3307    #    domain = Domain(coordinates, volumes,\
3308    #                    conserved_quantities = conserved_quantities,\
3309    #                    other_quantities = other_quantities,zone=zone,\
3310    #                    xllcorner=xllcorner, yllcorner=yllcorner)
3311
3312    #   From shallow_water.Domain:
3313    coordinates=coordinates.tolist()
3314    volumes=volumes.tolist()
3315    #FIXME:should this be in mesh?(peter row)
3316    if fid.smoothing == 'Yes': unique = False
3317    else: unique = True
3318    if unique:
3319        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
3320
3321
3322    try:
3323        domain = Domain(coordinates, volumes, boundary)
3324    except AssertionError, e:
3325        fid.close()
3326        msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
3327        raise DataDomainError, msg
3328
3329    if not boundary is None:
3330        domain.boundary = boundary
3331
3332    domain.geo_reference = geo_reference
3333
3334    domain.starttime=float(starttime)+float(t)
3335    domain.time=0.0
3336
3337    for quantity in other_quantities:
3338        try:
3339            NaN = fid.variables[quantity].missing_value
3340        except:
3341            pass #quantity has no missing_value number
3342        X = fid.variables[quantity][:]
3343        if very_verbose:
3344            print '       ',quantity
3345            print '        NaN =',NaN
3346            print '        max(X)'
3347            print '       ',max(X)
3348            print '        max(X)==NaN'
3349            print '       ',max(X)==NaN
3350            print ''
3351        if (max(X)==NaN) or (min(X)==NaN):
3352            if fail_if_NaN:
3353                msg = 'quantity "%s" contains no_data entry'%quantity
3354                raise DataMissingValuesError, msg
3355            else:
3356                data = (X<>NaN)
3357                X = (X*data)+(data==0)*NaN_filler
3358        if unique:
3359            X = resize(X,(len(X)/3,3))
3360        domain.set_quantity(quantity,X)
3361    #
3362    for quantity in conserved_quantities:
3363        try:
3364            NaN = fid.variables[quantity].missing_value
3365        except:
3366            pass #quantity has no missing_value number
3367        X = interpolated_quantities[quantity]
3368        if very_verbose:
3369            print '       ',quantity
3370            print '        NaN =',NaN
3371            print '        max(X)'
3372            print '       ',max(X)
3373            print '        max(X)==NaN'
3374            print '       ',max(X)==NaN
3375            print ''
3376        if (max(X)==NaN) or (min(X)==NaN):
3377            if fail_if_NaN:
3378                msg = 'quantity "%s" contains no_data entry'%quantity
3379                raise DataMissingValuesError, msg
3380            else:
3381                data = (X<>NaN)
3382                X = (X*data)+(data==0)*NaN_filler
3383        if unique:
3384            X = resize(X,(X.shape[0]/3,3))
3385        domain.set_quantity(quantity,X)
3386
3387    fid.close()
3388    return domain
3389
[5276]3390
[2852]3391def interpolated_quantity(saved_quantity,time_interp):
3392
3393    #given an index and ratio, interpolate quantity with respect to time.
3394    index,ratio = time_interp
3395    Q = saved_quantity
3396    if ratio > 0:
3397        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
3398    else:
3399        q = Q[index]
3400    #Return vector of interpolated values
3401    return q
3402
[5276]3403
[2852]3404def get_time_interp(time,t=None):
3405    #Finds the ratio and index for time interpolation.
[3560]3406    #It is borrowed from previous abstract_2d_finite_volumes code.
[2852]3407    if t is None:
3408        t=time[-1]
3409        index = -1
3410        ratio = 0.
3411    else:
3412        T = time
3413        tau = t
3414        index=0
3415        msg = 'Time interval derived from file %s [%s:%s]'\
3416            %('FIXMEfilename', T[0], T[-1])
3417        msg += ' does not match model time: %s' %tau
3418        if tau < time[0]: raise DataTimeError, msg
3419        if tau > time[-1]: raise DataTimeError, msg
3420        while tau > time[index]: index += 1
3421        while tau < time[index]: index -= 1
3422        if tau == time[index]:
3423            #Protect against case where tau == time[-1] (last time)
3424            # - also works in general when tau == time[i]
3425            ratio = 0
3426        else:
3427            #t is now between index and index+1
3428            ratio = (tau - time[index])/(time[index+1] - time[index])
3429    return (index,ratio)
3430
3431
3432def weed(coordinates,volumes,boundary = None):
[4455]3433    if type(coordinates)==ArrayType:
[2852]3434        coordinates = coordinates.tolist()
[4455]3435    if type(volumes)==ArrayType:
[2852]3436        volumes = volumes.tolist()
3437
3438    unique = False
3439    point_dict = {}
3440    same_point = {}
3441    for i in range(len(coordinates)):
3442        point = tuple(coordinates[i])
3443        if point_dict.has_key(point):
3444            unique = True
3445            same_point[i]=point
3446            #to change all point i references to point j
3447        else:
3448            point_dict[point]=i
3449            same_point[i]=point
3450
3451    coordinates = []
3452    i = 0
3453    for point in point_dict.keys():
3454        point = tuple(point)
3455        coordinates.append(list(point))
3456        point_dict[point]=i
3457        i+=1
3458
3459
3460    for volume in volumes:
3461        for i in range(len(volume)):
3462            index = volume[i]
3463            if index>-1:
3464                volume[i]=point_dict[same_point[index]]
3465
3466    new_boundary = {}
3467    if not boundary is None:
3468        for segment in boundary.keys():
3469            point0 = point_dict[same_point[segment[0]]]
3470            point1 = point_dict[same_point[segment[1]]]
3471            label = boundary[segment]
3472            #FIXME should the bounday attributes be concaterated
3473            #('exterior, pond') or replaced ('pond')(peter row)
3474
3475            if new_boundary.has_key((point0,point1)):
3476                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
3477                                              #+','+label
3478
3479            elif new_boundary.has_key((point1,point0)):
3480                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
3481                                              #+','+label
3482            else: new_boundary[(point0,point1)]=label
3483
3484        boundary = new_boundary
3485
3486    return coordinates,volumes,boundary
3487
3488
3489def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
3490                 verbose=False):
3491    """Read Digitial Elevation model from the following NetCDF format (.dem)
3492
3493    Example:
3494
3495    ncols         3121
3496    nrows         1800
3497    xllcorner     722000
3498    yllcorner     5893000
3499    cellsize      25
3500    NODATA_value  -9999
3501    138.3698 137.4194 136.5062 135.5558 ..........
3502
3503    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
3504    """
3505
3506    import os
3507    from Scientific.IO.NetCDF import NetCDFFile
3508    from Numeric import Float, zeros, sum, reshape, equal
3509
3510    root = basename_in
3511    inname = root + '.dem'
3512
3513    #Open existing netcdf file to read
3514    infile = NetCDFFile(inname, 'r')
3515    if verbose: print 'Reading DEM from %s' %inname
3516
3517    #Read metadata
3518    ncols = infile.ncols[0]
3519    nrows = infile.nrows[0]
3520    xllcorner = infile.xllcorner[0]
3521    yllcorner = infile.yllcorner[0]
3522    cellsize = infile.cellsize[0]
3523    NODATA_value = infile.NODATA_value[0]
3524    zone = infile.zone[0]
3525    false_easting = infile.false_easting[0]
3526    false_northing = infile.false_northing[0]
3527    projection = infile.projection
3528    datum = infile.datum
3529    units = infile.units
3530
3531    dem_elevation = infile.variables['elevation']
3532
3533    #Get output file name
3534    if basename_out == None:
3535        outname = root + '_' + repr(cellsize_new) + '.dem'
3536    else:
3537        outname = basename_out + '.dem'
3538
3539    if verbose: print 'Write decimated NetCDF file to %s' %outname
3540
3541    #Determine some dimensions for decimated grid
3542    (nrows_stencil, ncols_stencil) = stencil.shape
3543    x_offset = ncols_stencil / 2
3544    y_offset = nrows_stencil / 2
3545    cellsize_ratio = int(cellsize_new / cellsize)
3546    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
3547    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
3548
3549    #Open netcdf file for output
3550    outfile = NetCDFFile(outname, 'w')
3551
3552    #Create new file
3553    outfile.institution = 'Geoscience Australia'
3554    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
3555                           'of spatial point data'
3556    #Georeferencing
3557    outfile.zone = zone
3558    outfile.projection = projection
3559    outfile.datum = datum
3560    outfile.units = units
3561
3562    outfile.cellsize = cellsize_new
3563    outfile.NODATA_value = NODATA_value
3564    outfile.false_easting = false_easting
3565    outfile.false_northing = false_northing
3566
3567    outfile.xllcorner = xllcorner + (x_offset * cellsize)
3568    outfile.yllcorner = yllcorner + (y_offset * cellsize)
3569    outfile.ncols = ncols_new
3570    outfile.nrows = nrows_new
3571
3572    # dimension definition
3573    outfile.createDimension('number_of_points', nrows_new*ncols_new)
3574
3575    # variable definition
3576    outfile.createVariable('elevation', Float, ('number_of_points',))
3577
3578    # Get handle to the variable
3579    elevation = outfile.variables['elevation']
3580
3581    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
3582
3583    #Store data
3584    global_index = 0
3585    for i in range(nrows_new):
3586        if verbose: print 'Processing row %d of %d' %(i, nrows_new)
3587        lower_index = global_index
3588        telev =  zeros(ncols_new, Float)
3589        local_index = 0
3590        trow = i * cellsize_ratio
3591
3592        for j in range(ncols_new):
3593            tcol = j * cellsize_ratio
3594            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
3595
3596            #if dem contains 1 or more NODATA_values set value in
3597            #decimated dem to NODATA_value, else compute decimated
3598            #value using stencil
3599            if sum(sum(equal(tmp, NODATA_value))) > 0:
3600                telev[local_index] = NODATA_value
3601            else:
3602                telev[local_index] = sum(sum(tmp * stencil))
3603
3604            global_index += 1
3605            local_index += 1
3606
3607        upper_index = global_index
3608
3609        elevation[lower_index:upper_index] = telev
3610
3611    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
3612
3613    infile.close()
3614    outfile.close()
3615
3616
3617
3618
[4868]3619def tsh2sww(filename, verbose=False): 
[2852]3620    """
3621    to check if a tsh/msh file 'looks' good.
3622    """
3623
3624
3625    if verbose == True:print 'Creating domain from', filename
3626    domain = pmesh_to_domain_instance(filename, Domain)
3627    if verbose == True:print "Number of triangles = ", len(domain)
3628
3629    domain.smooth = True
3630    domain.format = 'sww'   #Native netcdf visualisation format
3631    file_path, filename = path.split(filename)
3632    filename, ext = path.splitext(filename)
[3846]3633    domain.set_name(filename)   
[2852]3634    domain.reduction = mean
3635    if verbose == True:print "file_path",file_path
3636    if file_path == "":file_path = "."
3637    domain.set_datadir(file_path)
3638
3639    if verbose == True:
3640        print "Output written to " + domain.get_datadir() + sep + \
[3846]3641              domain.get_name() + "." + domain.format
[2852]3642    sww = get_dataobject(domain)
3643    sww.store_connectivity()
[4868]3644    sww.store_timestep()
[2852]3645
3646
3647def asc_csiro2sww(bath_dir,
3648                  elevation_dir,
3649                  ucur_dir,
3650                  vcur_dir,
3651                  sww_file,
3652                  minlat = None, maxlat = None,
3653                  minlon = None, maxlon = None,
3654                  zscale=1,
3655                  mean_stage = 0,
3656                  fail_on_NaN = True,
3657                  elevation_NaN_filler = 0,
3658                  bath_prefix='ba',
3659                  elevation_prefix='el',
3660                  verbose=False):
3661    """
3662    Produce an sww boundary file, from esri ascii data from CSIRO.
3663
3664    Also convert latitude and longitude to UTM. All coordinates are
3665    assumed to be given in the GDA94 datum.
3666
3667    assume:
3668    All files are in esri ascii format
3669
3670    4 types of information
3671    bathymetry
3672    elevation
3673    u velocity
3674    v velocity
3675
3676    Assumptions
3677    The metadata of all the files is the same
3678    Each type is in a seperate directory
3679    One bath file with extention .000
3680    The time period is less than 24hrs and uniform.
3681    """
3682    from Scientific.IO.NetCDF import NetCDFFile
3683
[3514]3684    from anuga.coordinate_transforms.redfearn import redfearn
[2852]3685
3686    precision = Float # So if we want to change the precision its done here
3687
3688    # go in to the bath dir and load the only file,
3689    bath_files = os.listdir(bath_dir)
3690
3691    bath_file = bath_files[0]
3692    bath_dir_file =  bath_dir + os.sep + bath_file
3693    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
3694
3695    #Use the date.time of the bath file as a basis for
3696    #the start time for other files
3697    base_start = bath_file[-12:]
3698
3699    #go into the elevation dir and load the 000 file
3700    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
3701                         + base_start
3702
3703    elevation_files = os.listdir(elevation_dir)
3704    ucur_files = os.listdir(ucur_dir)
3705    vcur_files = os.listdir(vcur_dir)
[4031]3706    elevation_files.sort()
[2852]3707    # the first elevation file should be the
3708    # file with the same base name as the bath data
3709    assert elevation_files[0] == 'el' + base_start
3710
3711    number_of_latitudes = bath_grid.shape[0]
3712    number_of_longitudes = bath_grid.shape[1]
3713    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3714
3715    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
3716                  for x in range(number_of_longitudes)]
3717    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
3718                 for y in range(number_of_latitudes)]
3719
3720     # reverse order of lat, so the fist lat represents the first grid row
3721    latitudes.reverse()
3722
[4027]3723    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:],
[2852]3724                                                 minlat=minlat, maxlat=maxlat,
3725                                                 minlon=minlon, maxlon=maxlon)
3726
3727
3728    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
3729    latitudes = latitudes[kmin:kmax]
3730    longitudes = longitudes[lmin:lmax]
3731    number_of_latitudes = len(latitudes)
3732    number_of_longitudes = len(longitudes)
3733    number_of_times = len(os.listdir(elevation_dir))
3734    number_of_points = number_of_latitudes*number_of_longitudes
3735    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3736
3737    #Work out the times
3738    if len(elevation_files) > 1:
3739        # Assume: The time period is less than 24hrs.
3740        time_period = (int(elevation_files[1][-3:]) - \
3741                      int(elevation_files[0][-3:]))*60*60
3742        times = [x*time_period for x in range(len(elevation_files))]
3743    else:
3744        times = [0.0]
3745
3746
3747    if verbose:
3748        print '------------------------------------------------'
3749        print 'Statistics:'
3750        print '  Extent (lat/lon):'
3751        print '    lat in [%f, %f], len(lat) == %d'\
3752              %(min(latitudes), max(latitudes),
3753                len(latitudes))
3754        print '    lon in [%f, %f], len(lon) == %d'\
3755              %(min(longitudes), max(longitudes),
3756                len(longitudes))
3757        print '    t in [%f, %f], len(t) == %d'\
3758              %(min(times), max(times), len(times))
3759
3760    ######### WRITE THE SWW FILE #############
3761    # NetCDF file definition