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

Last change on this file since 4416 was 4416, checked in by duncan, 17 years ago

small bug fix

File size: 151.8 KB
Line 
1"""datamanager.py - input output for AnuGA
2
3
4This module takes care of reading and writing datafiles such as topograhies,
5model output, etc
6
7
8Formats used within AnuGA:
9
10.sww: Netcdf format for storing model output f(t,x,y)
11.tms: Netcdf format for storing time series f(t)
12
13.xya: ASCII format for storing arbitrary points and associated attributes
14.pts: NetCDF format for storing arbitrary points and associated attributes
15
16.asc: ASCII format of regular DEMs as output from ArcView
17.prj: Associated ArcView file giving more meta data for asc format
18.ers: ERMapper header format of regular DEMs for ArcView
19
20.dem: NetCDF representation of regular DEM data
21
22.tsh: ASCII format for storing meshes and associated boundary and region info
23.msh: NetCDF format for storing meshes and associated boundary and region info
24
25.nc: Native ferret NetCDF format
26.geo: Houdinis ascii geometry format (?)
27
28
29A typical dataflow can be described as follows
30
31Manually created files:
32ASC, PRJ:     Digital elevation models (gridded)
33TSH:          Triangular meshes (e.g. created from anuga.pmesh)
34NC            Model outputs for use as boundary conditions (e.g from MOST)
35
36
37AUTOMATICALLY CREATED FILES:
38
39ASC, PRJ  ->  DEM  ->  PTS: Conversion of DEM's to native pts file
40
41NC -> SWW: Conversion of MOST bundary files to boundary sww
42
43PTS + TSH -> TSH with elevation: Least squares fit
44
45TSH -> SWW:  Conversion of TSH to sww viewable using Swollen
46
47TSH + Boundary SWW -> SWW: Simluation using abstract_2d_finite_volumes
48
49"""
50
51import exceptions
52class TitleValueError(exceptions.Exception): pass
53class DataMissingValuesError(exceptions.Exception): pass
54class DataFileNotOpenError(exceptions.Exception): pass
55class DataTimeError(exceptions.Exception): pass
56class DataDomainError(exceptions.Exception): pass
57
58
59
60import csv
61import os
62from struct import unpack
63import array as p_array
64#import time, os
65from os import sep, path
66
67from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
68     searchsorted, zeros, allclose, around, reshape, transpose, sort
69from Scientific.IO.NetCDF import NetCDFFile
70
71
72from anuga.coordinate_transforms.redfearn import redfearn, \
73     convert_from_latlon_to_utm
74from anuga.coordinate_transforms.geo_reference import Geo_reference, \
75     write_NetCDF_georeference
76from anuga.geospatial_data.geospatial_data import Geospatial_data,\
77     ensure_absolute
78from anuga.config import minimum_storable_height as default_minimum_storable_height
79from anuga.utilities.numerical_tools import ensure_numeric,  mean
80from anuga.caching.caching import myhash
81from anuga.utilities.anuga_exceptions import ANUGAError
82from anuga.shallow_water import Domain
83from anuga.abstract_2d_finite_volumes.pmesh2domain import \
84     pmesh_to_domain_instance
85   
86def make_filename(s):
87    """Transform argument string into a suitable filename
88    """
89
90    s = s.strip()
91    s = s.replace(' ', '_')
92    s = s.replace('(', '')
93    s = s.replace(')', '')
94    s = s.replace('__', '_')
95
96    return s
97
98
99def check_dir(path, verbose=None):
100    """Check that specified path exists.
101    If path does not exist it will be created if possible
102
103    USAGE:
104       checkdir(path, verbose):
105
106    ARGUMENTS:
107        path -- Directory
108        verbose -- Flag verbose output (default: None)
109
110    RETURN VALUE:
111        Verified path including trailing separator
112
113    """
114
115    import os, sys
116    import os.path
117
118    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
119        unix = 0
120    else:
121        unix = 1
122
123
124    if path[-1] != os.sep:
125        path = path + os.sep  # Add separator for directories
126
127    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
128    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
129        try:
130            exitcode=os.mkdir(path)
131
132            # Change access rights if possible
133            #
134            if unix:
135                exitcode=os.system('chmod 775 '+path)
136            else:
137                pass  # FIXME: What about acces rights under Windows?
138
139            if verbose: print 'MESSAGE: Directory', path, 'created.'
140
141        except:
142            print 'WARNING: Directory', path, 'could not be created.'
143            if unix:
144                path = '/tmp/'
145            else:
146                path = 'C:'
147
148            print 'Using directory %s instead' %path
149
150    return(path)
151
152
153
154def del_dir(path):
155    """Recursively delete directory path and all its contents
156    """
157
158    import os
159
160    if os.path.isdir(path):
161        for file in os.listdir(path):
162            X = os.path.join(path, file)
163
164
165            if os.path.isdir(X) and not os.path.islink(X):
166                del_dir(X)
167            else:
168                try:
169                    os.remove(X)
170                except:
171                    print "Could not remove file %s" %X
172
173        os.rmdir(path)
174
175
176
177def create_filename(datadir, filename, format, size=None, time=None):
178
179    import os
180    #from anuga.config import data_dir
181
182    FN = check_dir(datadir) + filename
183
184    if size is not None:
185        FN += '_size%d' %size
186
187    if time is not None:
188        FN += '_time%.2f' %time
189
190    FN += '.' + format
191    return FN
192
193
194def get_files(datadir, filename, format, size):
195    """Get all file (names) with gven name, size and format
196    """
197
198    import glob
199
200    import os
201    #from anuga.config import data_dir
202
203    dir = check_dir(datadir)
204
205    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
206    return glob.glob(pattern)
207
208
209
210#Generic class for storing output to e.g. visualisation or checkpointing
211class Data_format:
212    """Generic interface to data formats
213    """
214
215
216    def __init__(self, domain, extension, mode = 'w'):
217        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
218                                        '''   'w' (write)'''+\
219                                        '''   'r' (read)''' +\
220                                        '''   'a' (append)'''
221
222        #Create filename
223        self.filename = create_filename(domain.get_datadir(),
224                                        domain.get_name(), extension)
225
226        #print 'F', self.filename
227        self.timestep = 0
228        self.domain = domain
229       
230
231
232        # Exclude ghosts in case this is a parallel domain
233        self.number_of_nodes = domain.number_of_full_nodes       
234        self.number_of_volumes = domain.number_of_full_triangles
235        #self.number_of_volumes = len(domain)       
236
237
238
239
240        #FIXME: Should we have a general set_precision function?
241
242
243
244#Class for storing output to e.g. visualisation
245class Data_format_sww(Data_format):
246    """Interface to native NetCDF format (.sww) for storing model output
247
248    There are two kinds of data
249
250    1: Constant data: Vertex coordinates and field values. Stored once
251    2: Variable data: Conserved quantities. Stored once per timestep.
252
253    All data is assumed to reside at vertex locations.
254    """
255
256
257    def __init__(self, domain, mode = 'w',\
258                 max_size = 2000000000,
259                 recursion = False):
260        from Scientific.IO.NetCDF import NetCDFFile
261        from Numeric import Int, Float, Float32
262
263        self.precision = Float32 #Use single precision for quantities
264        if hasattr(domain, 'max_size'):
265            self.max_size = domain.max_size #file size max is 2Gig
266        else:
267            self.max_size = max_size
268        self.recursion = recursion
269        self.mode = mode
270
271        Data_format.__init__(self, domain, 'sww', mode)
272
273        if hasattr(domain, 'minimum_storable_height'):
274            self.minimum_storable_height =  domain.minimum_storable_height
275        else:
276            self.minimum_storable_height = default_minimum_storable_height
277
278        # NetCDF file definition
279        fid = NetCDFFile(self.filename, mode)
280
281        if mode == 'w':
282
283            #Create new file
284            fid.institution = 'Geoscience Australia'
285            fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
286
287            if domain.smooth:
288                fid.smoothing = 'Yes'
289            else:
290                fid.smoothing = 'No'
291
292            fid.order = domain.default_order
293
294            if hasattr(domain, 'texture'):
295                fid.texture = domain.texture
296            #else:
297            #    fid.texture = 'None'
298
299            #Reference point
300            #Start time in seconds since the epoch (midnight 1/1/1970)
301            #FIXME: Use Georef
302            fid.starttime = domain.starttime
303
304            # dimension definitions
305            fid.createDimension('number_of_volumes', self.number_of_volumes)
306            fid.createDimension('number_of_vertices', 3)
307
308            if domain.smooth is True:
309                #fid.createDimension('number_of_points', len(domain.vertexlist))
310                fid.createDimension('number_of_points', self.domain.number_of_full_nodes)
311
312                # FIXME(Ole): This will cause sww files for paralle domains to
313                # have ghost nodes stored (but not used by triangles).
314                # To clean this up, we have to change get_vertex_values and friends in
315                # quantity.py (but I can't be bothered right now)
316            else:
317                fid.createDimension('number_of_points', 3*self.number_of_volumes)
318
319            fid.createDimension('number_of_timesteps', None) #extensible
320
321            # variable definitions
322            fid.createVariable('x', self.precision, ('number_of_points',))
323            fid.createVariable('y', self.precision, ('number_of_points',))
324            fid.createVariable('elevation', self.precision, ('number_of_points',))
325            if domain.geo_reference is not None:
326                domain.geo_reference.write_NetCDF(fid)
327
328            #FIXME: Backwards compatibility
329            fid.createVariable('z', self.precision, ('number_of_points',))
330            #################################
331
332            fid.createVariable('volumes', Int, ('number_of_volumes',
333                                                'number_of_vertices'))
334
335            fid.createVariable('time', Float,  # Always use full precision lest two timesteps
336                                               # close to each other may appear as the same step
337                               ('number_of_timesteps',))
338
339            fid.createVariable('stage', self.precision,
340                               ('number_of_timesteps',
341                                'number_of_points'))
342
343            fid.createVariable('xmomentum', self.precision,
344                               ('number_of_timesteps',
345                                'number_of_points'))
346
347            fid.createVariable('ymomentum', self.precision,
348                               ('number_of_timesteps',
349                                'number_of_points'))
350
351        #Close
352        fid.close()
353
354
355    def store_connectivity(self):
356        """Specialisation of store_connectivity for net CDF format
357
358        Writes x,y,z coordinates of triangles constituting
359        the bed elevation.
360        """
361
362        from Scientific.IO.NetCDF import NetCDFFile
363
364        from Numeric import concatenate, Int
365
366        domain = self.domain
367
368        #Get NetCDF
369        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
370
371        # Get the variables
372        x = fid.variables['x']
373        y = fid.variables['y']
374        z = fid.variables['elevation']
375
376        volumes = fid.variables['volumes']
377
378        # Get X, Y and bed elevation Z
379        Q = domain.quantities['elevation']
380        X,Y,Z,V = Q.get_vertex_values(xy=True,
381                                      precision=self.precision)
382
383        volumes[:] = V.astype(volumes.typecode())
384        x[:] = X
385        y[:] = Y
386        z[:] = Z
387
388        #FIXME: Backwards compatibility
389        z = fid.variables['z']
390        z[:] = Z
391        ################################
392
393
394
395        #Close
396        fid.close()
397
398    def store_timestep(self, names):
399        """Store time and named quantities to file
400        """
401        from Scientific.IO.NetCDF import NetCDFFile
402        import types
403        from time import sleep
404        from os import stat
405
406        from Numeric import choose
407
408        #Get NetCDF
409        retries = 0
410        file_open = False
411        while not file_open and retries < 10:
412            try:
413                fid = NetCDFFile(self.filename, 'a') #Open existing file
414            except IOError:
415                #This could happen if someone was reading the file.
416                #In that case, wait a while and try again
417                msg = 'Warning (store_timestep): File %s could not be opened'\
418                      %self.filename
419                msg += ' - trying step %s again' %self.domain.time
420                print msg
421                retries += 1
422                sleep(1)
423            else:
424                file_open = True
425
426        if not file_open:
427            msg = 'File %s could not be opened for append' %self.filename
428            raise DataFileNotOpenError, msg
429
430
431
432        #Check to see if the file is already too big:
433        time = fid.variables['time']
434        i = len(time)+1
435        file_size = stat(self.filename)[6]
436        file_size_increase =  file_size/i
437        if file_size + file_size_increase > self.max_size*(2**self.recursion):
438            #in order to get the file name and start time correct,
439            #I change the domian.filename and domain.starttime.
440            #This is the only way to do this without changing
441            #other modules (I think).
442
443            #write a filename addon that won't break swollens reader
444            #(10.sww is bad)
445            filename_ext = '_time_%s'%self.domain.time
446            filename_ext = filename_ext.replace('.', '_')
447            #remember the old filename, then give domain a
448            #name with the extension
449            old_domain_filename = self.domain.get_name()
450            if not self.recursion:
451                self.domain.set_name(old_domain_filename+filename_ext)
452
453
454            #change the domain starttime to the current time
455            old_domain_starttime = self.domain.starttime
456            self.domain.starttime = self.domain.time
457
458            #build a new data_structure.
459            next_data_structure=\
460                Data_format_sww(self.domain, mode=self.mode,\
461                                max_size = self.max_size,\
462                                recursion = self.recursion+1)
463            if not self.recursion:
464                print '    file_size = %s'%file_size
465                print '    saving file to %s'%next_data_structure.filename
466            #set up the new data_structure
467            self.domain.writer = next_data_structure
468
469            #FIXME - could be cleaner to use domain.store_timestep etc.
470            next_data_structure.store_connectivity()
471            next_data_structure.store_timestep(names)
472            fid.sync()
473            fid.close()
474
475            #restore the old starttime and filename
476            self.domain.starttime = old_domain_starttime
477            self.domain.set_name(old_domain_filename)           
478        else:
479            self.recursion = False
480            domain = self.domain
481
482            # Get the variables
483            time = fid.variables['time']
484            stage = fid.variables['stage']
485            xmomentum = fid.variables['xmomentum']
486            ymomentum = fid.variables['ymomentum']
487            i = len(time)
488
489            #Store time
490            time[i] = self.domain.time
491
492
493            if type(names) not in [types.ListType, types.TupleType]:
494                names = [names]
495
496            for name in names:
497                # Get quantity
498                Q = domain.quantities[name]
499                A,V = Q.get_vertex_values(xy = False,
500                                          precision = self.precision)
501
502
503                #FIXME: Make this general (see below)
504                if name == 'stage':
505                    z = fid.variables['elevation']
506                    #print z[:]
507                    #print A-z[:]
508
509
510                    A = choose(A-z[:] >= self.minimum_storable_height,
511                               (z[:], A))
512                    stage[i,:] = A.astype(self.precision)
513                elif name == 'xmomentum':
514                    xmomentum[i,:] = A.astype(self.precision)
515                elif name == 'ymomentum':
516                    ymomentum[i,:] = A.astype(self.precision)
517
518                #As in....
519            #eval( name + '[i,:] = A.astype(self.precision)' )
520            #FIXME: But we need a UNIT test for that before refactoring
521
522
523
524            #Flush and close
525            fid.sync()
526            fid.close()
527
528
529
530#Class for handling checkpoints data
531class Data_format_cpt(Data_format):
532    """Interface to native NetCDF format (.cpt)
533    """
534
535
536    def __init__(self, domain, mode = 'w'):
537        from Scientific.IO.NetCDF import NetCDFFile
538        from Numeric import Int, Float, Float
539
540        self.precision = Float #Use full precision
541
542        Data_format.__init__(self, domain, 'sww', mode)
543
544
545        # NetCDF file definition
546        fid = NetCDFFile(self.filename, mode)
547
548        if mode == 'w':
549            #Create new file
550            fid.institution = 'Geoscience Australia'
551            fid.description = 'Checkpoint data'
552            #fid.smooth = domain.smooth
553            fid.order = domain.default_order
554
555            # dimension definitions
556            fid.createDimension('number_of_volumes', self.number_of_volumes)
557            fid.createDimension('number_of_vertices', 3)
558
559            #Store info at all vertices (no smoothing)
560            fid.createDimension('number_of_points', 3*self.number_of_volumes)
561            fid.createDimension('number_of_timesteps', None) #extensible
562
563            # variable definitions
564
565            #Mesh
566            fid.createVariable('x', self.precision, ('number_of_points',))
567            fid.createVariable('y', self.precision, ('number_of_points',))
568
569
570            fid.createVariable('volumes', Int, ('number_of_volumes',
571                                                'number_of_vertices'))
572
573            fid.createVariable('time', self.precision,
574                               ('number_of_timesteps',))
575
576            #Allocate space for all quantities
577            for name in domain.quantities.keys():
578                fid.createVariable(name, self.precision,
579                                   ('number_of_timesteps',
580                                    'number_of_points'))
581
582        #Close
583        fid.close()
584
585
586    def store_checkpoint(self):
587        """
588        Write x,y coordinates of triangles.
589        Write connectivity (
590        constituting
591        the bed elevation.
592        """
593
594        from Scientific.IO.NetCDF import NetCDFFile
595
596        from Numeric import concatenate
597
598        domain = self.domain
599
600        #Get NetCDF
601        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
602
603        # Get the variables
604        x = fid.variables['x']
605        y = fid.variables['y']
606
607        volumes = fid.variables['volumes']
608
609        # Get X, Y and bed elevation Z
610        Q = domain.quantities['elevation']
611        X,Y,Z,V = Q.get_vertex_values(xy=True,
612                      precision = self.precision)
613
614
615
616        x[:] = X.astype(self.precision)
617        y[:] = Y.astype(self.precision)
618        z[:] = Z.astype(self.precision)
619
620        volumes[:] = V
621
622        #Close
623        fid.close()
624
625
626    def store_timestep(self, name):
627        """Store time and named quantity to file
628        """
629        from Scientific.IO.NetCDF import NetCDFFile
630        from time import sleep
631
632        #Get NetCDF
633        retries = 0
634        file_open = False
635        while not file_open and retries < 10:
636            try:
637                fid = NetCDFFile(self.filename, 'a') #Open existing file
638            except IOError:
639                #This could happen if someone was reading the file.
640                #In that case, wait a while and try again
641                msg = 'Warning (store_timestep): File %s could not be opened'\
642                  %self.filename
643                msg += ' - trying again'
644                print msg
645                retries += 1
646                sleep(1)
647            else:
648                file_open = True
649
650        if not file_open:
651            msg = 'File %s could not be opened for append' %self.filename
652            raise DataFileNotOPenError, msg
653
654
655        domain = self.domain
656
657        # Get the variables
658        time = fid.variables['time']
659        stage = fid.variables['stage']
660        i = len(time)
661
662        #Store stage
663        time[i] = self.domain.time
664
665        # Get quantity
666        Q = domain.quantities[name]
667        A,V = Q.get_vertex_values(xy=False,
668                  precision = self.precision)
669
670        stage[i,:] = A.astype(self.precision)
671
672        #Flush and close
673        fid.sync()
674        fid.close()
675
676
677#### NED is national exposure database (name changed to NEXIS)
678   
679LAT_TITLE = 'LATITUDE'
680LONG_TITLE = 'LONGITUDE'
681X_TITLE = 'x'
682Y_TITLE = 'y'
683class Exposure_csv:
684    def __init__(self,file_name, latitude_title=LAT_TITLE,
685                 longitude_title=LONG_TITLE, is_x_y_locations=None,
686                 x_title=X_TITLE, y_title=Y_TITLE,
687                 refine_polygon=None, title_check_list=None):
688        """
689        This class is for handling the exposure csv file.
690        It reads the file in and converts the lats and longs to a geospatial
691        data object.
692        Use the methods to read and write columns.
693
694        The format of the csv files it reads is;
695           The first row is a title row.
696           comma's are the delimiters
697           each column is a 'set' of data
698
699        Feel free to use/expand it to read other csv files.
700           
701           
702        It is not for adding and deleting rows
703       
704        Can geospatial handle string attributes? It's not made for them.
705        Currently it can't load and save string att's.
706
707        So just use geospatial to hold the x, y and georef? Bad, since
708        different att's are in diferent structures.  Not so bad, the info
709        to write if the .csv file is saved is in attribute_dic
710
711        The location info is in the geospatial attribute.
712       
713       
714        """
715        self._file_name = file_name
716        self._geospatial = None #
717
718        # self._attribute_dic is a dictionary.
719        #The keys are the column titles.
720        #The values are lists of column data
721       
722        # self._title_index_dic is a dictionary.
723        #The keys are the column titles.
724        #The values are the index positions of file columns.
725        self._attribute_dic, self._title_index_dic = \
726        self._load_exposure_csv(self._file_name,
727                                title_check_list=title_check_list)
728        try:
729            #Have code here that handles caps or lower
730            lats = self._attribute_dic[latitude_title]
731            longs = self._attribute_dic[longitude_title]
732           
733        except KeyError:
734            # maybe a warning..
735            #Let's see if this works..
736            if False != is_x_y_locations:
737                is_x_y_locations = True
738            pass
739        else:
740            self._geospatial = Geospatial_data(latitudes = lats,
741                 longitudes = longs)
742
743        if is_x_y_locations is True:
744            if self._geospatial is not None:
745                pass #fixme throw an error
746            try:
747                xs = self._attribute_dic[x_title]
748                ys = self._attribute_dic[y_title]
749                points = [[float(i),float(j)] for i,j in map(None,xs,ys)]
750            except KeyError:
751                # maybe a warning..
752                msg = "Could not find location information."
753                raise TitleValueError, msg
754            else:
755                self._geospatial = Geospatial_data(data_points=points)
756           
757        # create a list of points that are in the refining_polygon
758        # described by a list of indexes representing the points
759
760    def __cmp__(self, other):
761        #print "self._attribute_dic",self._attribute_dic
762        #print "other._attribute_dic",other._attribute_dic
763        #print "self._title_index_dic", self._title_index_dic
764        #print "other._title_index_dic", other._title_index_dic
765       
766        #check that a is an instance of this class
767        if isinstance(self, type(other)):
768            result = cmp(self._attribute_dic, other._attribute_dic)
769            if result <>0:
770                return result
771            # The order of the columns is important. Therefore..
772            result = cmp(self._title_index_dic, other._title_index_dic)
773            if result <>0:
774                return result
775            for self_ls, other_ls in map(None,self._attribute_dic, \
776                    other._attribute_dic):
777                result = cmp(self._attribute_dic[self_ls],
778                             other._attribute_dic[other_ls])
779                if result <>0:
780                    return result
781            return 0
782        else:
783            return 1
784   
785    def _load_exposure_csv(self, file_name, title_check_list=None):
786        """
787        Load in the csv as a dic, title as key and column info as value, .
788        Also, create a dic, title as key and column index as value,
789        to keep track of the column order.
790        """
791        #
792        attribute_dic = {}
793        title_index_dic = {}
794        titles_stripped = [] # list of titles
795        reader = csv.reader(file(file_name))
796
797        # Read in and manipulate the title info
798        titles = reader.next()
799        for i,title in enumerate(titles):
800            titles_stripped.append(title.strip())
801            title_index_dic[title.strip()] = i
802        title_count = len(titles_stripped)       
803        #print "title_index_dic",title_index_dic
804        if title_check_list is not None:
805            for title_check in title_check_list:
806                #msg = "Reading error.  This row is not present ", title_check
807                #assert title_index_dic.has_key(title_check), msg
808                if not title_index_dic.has_key(title_check):
809                    #reader.close()
810                    msg = "Reading error.  This row is not present ", \
811                          title_check                     
812                    raise IOError, msg
813                   
814       
815       
816        #create a dic of colum values, indexed by column title
817        for line in reader:
818            if len(line) <> title_count:
819                raise IOError #FIXME make this nicer
820            for i, value in enumerate(line):
821                attribute_dic.setdefault(titles_stripped[i],[]).append(value)
822           
823        return attribute_dic, title_index_dic
824
825    def get_column(self, column_name, use_refind_polygon=False):
826        """
827        Given a column name return a list of the column values
828
829        Note, the type of the values will be String!
830        do this to change a list of strings to a list of floats
831        time = [float(x) for x in time]
832       
833        Not implemented:
834        if use_refind_polygon is True, only return values in the
835        refined polygon
836        """
837        if not self._attribute_dic.has_key(column_name):
838            msg = 'Therer is  no column called %s!' %column_name
839            raise TitleValueError, msg
840        return self._attribute_dic[column_name]
841
842
843    def get_value(self, value_column_name,
844                  known_column_name,
845                  known_values,
846                  use_refind_polygon=False):
847        """
848        Do linear interpolation on the known_colum, using the known_value,
849        to return a value of the column_value_name.
850        """
851        pass
852
853
854    def get_location(self, use_refind_polygon=False):
855        """
856        Return a geospatial object which describes the
857        locations of the location file.
858
859        Note, if there is not location info, this returns None.
860       
861        Not implemented:
862        if use_refind_polygon is True, only return values in the
863        refined polygon
864        """
865        return self._geospatial
866
867    def set_column(self, column_name, column_values, overwrite=False):
868        """
869        Add a column to the 'end' (with the right most column being the end)
870        of the csv file.
871
872        Set overwrite to True if you want to overwrite a column.
873       
874        Note, in column_name white space is removed and case is not checked.
875        Precondition
876        The column_name and column_values cannot have comma's in it.
877        """
878       
879        value_row_count = \
880            len(self._attribute_dic[self._title_index_dic.keys()[0]])
881        if len(column_values) <> value_row_count: 
882            msg = 'The number of column values must equal the number of rows.'
883            raise DataMissingValuesError, msg
884       
885        if self._attribute_dic.has_key(column_name):
886            if not overwrite:
887                msg = 'Column name %s already in use!' %column_name
888                raise TitleValueError, msg
889        else:
890            # New title.  Add it to the title index.
891            self._title_index_dic[column_name] = len(self._title_index_dic)
892        self._attribute_dic[column_name] = column_values
893        #print "self._title_index_dic[column_name]",self._title_index_dic
894
895    def save(self, file_name=None):
896        """
897        Save the exposure csv file
898        """
899        if file_name is None:
900            file_name = self._file_name
901       
902        fd = open(file_name,'wb')
903        writer = csv.writer(fd)
904       
905        #Write the title to a cvs file
906        line = [None]* len(self._title_index_dic)
907        for title in self._title_index_dic.iterkeys():
908            line[self._title_index_dic[title]]= title
909        writer.writerow(line)
910       
911        # Write the values to a cvs file
912        value_row_count = \
913            len(self._attribute_dic[self._title_index_dic.keys()[0]])
914        for row_i in range(value_row_count):
915            line = [None]* len(self._title_index_dic)
916            for title in self._title_index_dic.iterkeys():
917                line[self._title_index_dic[title]]= \
918                     self._attribute_dic[title][row_i]
919            writer.writerow(line)
920
921
922#Auxiliary
923def write_obj(filename,x,y,z):
924    """Store x,y,z vectors into filename (obj format)
925       Vectors are assumed to have dimension (M,3) where
926       M corresponds to the number elements.
927       triangles are assumed to be disconnected
928
929       The three numbers in each vector correspond to three vertices,
930
931       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
932
933    """
934    #print 'Writing obj to %s' % filename
935
936    import os.path
937
938    root, ext = os.path.splitext(filename)
939    if ext == '.obj':
940        FN = filename
941    else:
942        FN = filename + '.obj'
943
944
945    outfile = open(FN, 'wb')
946    outfile.write("# Triangulation as an obj file\n")
947
948    M, N = x.shape
949    assert N==3  #Assuming three vertices per element
950
951    for i in range(M):
952        for j in range(N):
953            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
954
955    for i in range(M):
956        base = i*N
957        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
958
959    outfile.close()
960
961
962#########################################################
963#Conversion routines
964########################################################
965
966def sww2obj(basefilename, size):
967    """Convert netcdf based data output to obj
968    """
969    from Scientific.IO.NetCDF import NetCDFFile
970
971    from Numeric import Float, zeros
972
973    #Get NetCDF
974    FN = create_filename('.', basefilename, 'sww', size)
975    print 'Reading from ', FN
976    fid = NetCDFFile(FN, 'r')  #Open existing file for read
977
978
979    # Get the variables
980    x = fid.variables['x']
981    y = fid.variables['y']
982    z = fid.variables['elevation']
983    time = fid.variables['time']
984    stage = fid.variables['stage']
985
986    M = size  #Number of lines
987    xx = zeros((M,3), Float)
988    yy = zeros((M,3), Float)
989    zz = zeros((M,3), Float)
990
991    for i in range(M):
992        for j in range(3):
993            xx[i,j] = x[i+j*M]
994            yy[i,j] = y[i+j*M]
995            zz[i,j] = z[i+j*M]
996
997    #Write obj for bathymetry
998    FN = create_filename('.', basefilename, 'obj', size)
999    write_obj(FN,xx,yy,zz)
1000
1001
1002    #Now read all the data with variable information, combine with
1003    #x,y info and store as obj
1004
1005    for k in range(len(time)):
1006        t = time[k]
1007        print 'Processing timestep %f' %t
1008
1009        for i in range(M):
1010            for j in range(3):
1011                zz[i,j] = stage[k,i+j*M]
1012
1013
1014        #Write obj for variable data
1015        #FN = create_filename(basefilename, 'obj', size, time=t)
1016        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
1017        write_obj(FN,xx,yy,zz)
1018
1019
1020def dat2obj(basefilename):
1021    """Convert line based data output to obj
1022    FIXME: Obsolete?
1023    """
1024
1025    import glob, os
1026    from anuga.config import data_dir
1027
1028
1029    #Get bathymetry and x,y's
1030    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
1031
1032    from Numeric import zeros, Float
1033
1034    M = len(lines)  #Number of lines
1035    x = zeros((M,3), Float)
1036    y = zeros((M,3), Float)
1037    z = zeros((M,3), Float)
1038
1039    ##i = 0
1040    for i, line in enumerate(lines):
1041        tokens = line.split()
1042        values = map(float,tokens)
1043
1044        for j in range(3):
1045            x[i,j] = values[j*3]
1046            y[i,j] = values[j*3+1]
1047            z[i,j] = values[j*3+2]
1048
1049        ##i += 1
1050
1051
1052    #Write obj for bathymetry
1053    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
1054
1055
1056    #Now read all the data files with variable information, combine with
1057    #x,y info
1058    #and store as obj
1059
1060    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
1061
1062    for filename in files:
1063        print 'Processing %s' % filename
1064
1065        lines = open(data_dir+os.sep+filename,'r').readlines()
1066        assert len(lines) == M
1067        root, ext = os.path.splitext(filename)
1068
1069        #Get time from filename
1070        i0 = filename.find('_time=')
1071        if i0 == -1:
1072            #Skip bathymetry file
1073            continue
1074
1075        i0 += 6  #Position where time starts
1076        i1 = filename.find('.dat')
1077
1078        if i1 > i0:
1079            t = float(filename[i0:i1])
1080        else:
1081            raise DataTimeError, 'Hmmmm'
1082
1083
1084
1085        ##i = 0
1086        for i, line in enumerate(lines):
1087            tokens = line.split()
1088            values = map(float,tokens)
1089
1090            for j in range(3):
1091                z[i,j] = values[j]
1092
1093            ##i += 1
1094
1095        #Write obj for variable data
1096        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
1097
1098
1099def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
1100    """Read netcdf filename1, pick timesteps first:step:last and save to
1101    nettcdf file filename2
1102    """
1103    from Scientific.IO.NetCDF import NetCDFFile
1104
1105    #Get NetCDF
1106    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
1107    outfile = NetCDFFile(filename2, 'w')  #Open new file
1108
1109
1110    #Copy dimensions
1111    for d in infile.dimensions:
1112        outfile.createDimension(d, infile.dimensions[d])
1113
1114    for name in infile.variables:
1115        var = infile.variables[name]
1116        outfile.createVariable(name, var.typecode(), var.dimensions)
1117
1118
1119    #Copy the static variables
1120    for name in infile.variables:
1121        if name == 'time' or name == 'stage':
1122            pass
1123        else:
1124            #Copy
1125            outfile.variables[name][:] = infile.variables[name][:]
1126
1127    #Copy selected timesteps
1128    time = infile.variables['time']
1129    stage = infile.variables['stage']
1130
1131    newtime = outfile.variables['time']
1132    newstage = outfile.variables['stage']
1133
1134    if last is None:
1135        last = len(time)
1136
1137    selection = range(first, last, step)
1138    for i, j in enumerate(selection):
1139        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
1140        newtime[i] = time[j]
1141        newstage[i,:] = stage[j,:]
1142
1143    #Close
1144    infile.close()
1145    outfile.close()
1146
1147
1148#Get data objects
1149def get_dataobject(domain, mode='w'):
1150    """Return instance of class of given format using filename
1151    """
1152
1153    cls = eval('Data_format_%s' %domain.format)
1154    return cls(domain, mode)
1155
1156
1157
1158
1159def dem2pts(basename_in, basename_out=None,
1160            easting_min=None, easting_max=None,
1161            northing_min=None, northing_max=None,
1162            use_cache=False, verbose=False,):
1163    """Read Digitial Elevation model from the following NetCDF format (.dem)
1164
1165    Example:
1166
1167    ncols         3121
1168    nrows         1800
1169    xllcorner     722000
1170    yllcorner     5893000
1171    cellsize      25
1172    NODATA_value  -9999
1173    138.3698 137.4194 136.5062 135.5558 ..........
1174
1175    Convert to NetCDF pts format which is
1176
1177    points:  (Nx2) Float array
1178    elevation: N Float array
1179    """
1180
1181
1182
1183    kwargs = {'basename_out': basename_out,
1184              'easting_min': easting_min,
1185              'easting_max': easting_max,
1186              'northing_min': northing_min,
1187              'northing_max': northing_max,
1188              'verbose': verbose}
1189
1190    if use_cache is True:
1191        from caching import cache
1192        result = cache(_dem2pts, basename_in, kwargs,
1193                       dependencies = [basename_in + '.dem'],
1194                       verbose = verbose)
1195
1196    else:
1197        result = apply(_dem2pts, [basename_in], kwargs)
1198
1199    return result
1200
1201
1202def _dem2pts(basename_in, basename_out=None, verbose=False,
1203            easting_min=None, easting_max=None,
1204            northing_min=None, northing_max=None):
1205    """Read Digitial Elevation model from the following NetCDF format (.dem)
1206
1207    Internal function. See public function dem2pts for details.
1208    """
1209
1210    #FIXME: Can this be written feasibly using write_pts?
1211
1212    import os
1213    from Scientific.IO.NetCDF import NetCDFFile
1214    from Numeric import Float, zeros, reshape, sum
1215
1216    root = basename_in
1217
1218    #Get NetCDF
1219    infile = NetCDFFile(root + '.dem', 'r')  #Open existing netcdf file for read
1220
1221    if verbose: print 'Reading DEM from %s' %(root + '.dem')
1222
1223    ncols = infile.ncols[0]
1224    nrows = infile.nrows[0]
1225    xllcorner = infile.xllcorner[0]  #Easting of lower left corner
1226    yllcorner = infile.yllcorner[0]  #Northing of lower left corner
1227    cellsize = infile.cellsize[0]
1228    NODATA_value = infile.NODATA_value[0]
1229    dem_elevation = infile.variables['elevation']
1230
1231    zone = infile.zone[0]
1232    false_easting = infile.false_easting[0]
1233    false_northing = infile.false_northing[0]
1234
1235    #Text strings
1236    projection = infile.projection
1237    datum = infile.datum
1238    units = infile.units
1239
1240
1241    #Get output file
1242    if basename_out == None:
1243        ptsname = root + '.pts'
1244    else:
1245        ptsname = basename_out + '.pts'
1246
1247    if verbose: print 'Store to NetCDF file %s' %ptsname
1248    # NetCDF file definition
1249    outfile = NetCDFFile(ptsname, 'w')
1250
1251    #Create new file
1252    outfile.institution = 'Geoscience Australia'
1253    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
1254                      'of spatial point data'
1255    #assign default values
1256    if easting_min is None: easting_min = xllcorner
1257    if easting_max is None: easting_max = xllcorner + ncols*cellsize
1258    if northing_min is None: northing_min = yllcorner
1259    if northing_max is None: northing_max = yllcorner + nrows*cellsize
1260
1261    #compute offsets to update georeferencing
1262    easting_offset = xllcorner - easting_min
1263    northing_offset = yllcorner - northing_min
1264
1265    #Georeferencing
1266    outfile.zone = zone
1267    outfile.xllcorner = easting_min #Easting of lower left corner
1268    outfile.yllcorner = northing_min #Northing of lower left corner
1269    outfile.false_easting = false_easting
1270    outfile.false_northing = false_northing
1271
1272    outfile.projection = projection
1273    outfile.datum = datum
1274    outfile.units = units
1275
1276
1277    #Grid info (FIXME: probably not going to be used, but heck)
1278    outfile.ncols = ncols
1279    outfile.nrows = nrows
1280
1281    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
1282    totalnopoints = nrows*ncols
1283
1284    # calculating number of NODATA_values for each row in clipped region
1285    #FIXME: use array operations to do faster
1286    nn = 0
1287    k = 0
1288    i1_0 = 0
1289    j1_0 = 0
1290    thisj = 0
1291    thisi = 0
1292    for i in range(nrows):
1293        y = (nrows-i-1)*cellsize + yllcorner
1294        for j in range(ncols):
1295            x = j*cellsize + xllcorner
1296            if easting_min <= x <= easting_max and \
1297               northing_min <= y <= northing_max:
1298                thisj = j
1299                thisi = i
1300                if dem_elevation_r[i,j] == NODATA_value: nn += 1
1301
1302                if k == 0:
1303                    i1_0 = i
1304                    j1_0 = j
1305                k += 1
1306
1307    index1 = j1_0
1308    index2 = thisj
1309
1310    # dimension definitions
1311    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
1312    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
1313
1314    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
1315    nopoints = clippednopoints-nn
1316
1317    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
1318
1319    if verbose:
1320        print 'There are %d values in the elevation' %totalnopoints
1321        print 'There are %d values in the clipped elevation' %clippednopoints
1322        print 'There are %d NODATA_values in the clipped elevation' %nn
1323
1324    outfile.createDimension('number_of_points', nopoints)
1325    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1326
1327    # variable definitions
1328    outfile.createVariable('points', Float, ('number_of_points',
1329                                             'number_of_dimensions'))
1330    outfile.createVariable('elevation', Float, ('number_of_points',))
1331
1332    # Get handles to the variables
1333    points = outfile.variables['points']
1334    elevation = outfile.variables['elevation']
1335
1336    lenv = index2-index1+1
1337    #Store data
1338    global_index = 0
1339    #for i in range(nrows):
1340    for i in range(i1_0,thisi+1,1):
1341        if verbose and i%((nrows+10)/10)==0:
1342            print 'Processing row %d of %d' %(i, nrows)
1343
1344        lower_index = global_index
1345
1346        v = dem_elevation_r[i,index1:index2+1]
1347        no_NODATA = sum(v == NODATA_value)
1348        if no_NODATA > 0:
1349            newcols = lenv - no_NODATA #ncols_in_bounding_box - no_NODATA
1350        else:
1351            newcols = lenv #ncols_in_bounding_box
1352
1353        telev = zeros(newcols, Float)
1354        tpoints = zeros((newcols, 2), Float)
1355
1356        local_index = 0
1357
1358        y = (nrows-i-1)*cellsize + yllcorner
1359        #for j in range(ncols):
1360        for j in range(j1_0,index2+1,1):
1361
1362            x = j*cellsize + xllcorner
1363            if easting_min <= x <= easting_max and \
1364               northing_min <= y <= northing_max and \
1365               dem_elevation_r[i,j] <> NODATA_value:
1366                tpoints[local_index, :] = [x-easting_min,y-northing_min]
1367                telev[local_index] = dem_elevation_r[i, j]
1368                global_index += 1
1369                local_index += 1
1370
1371        upper_index = global_index
1372
1373        if upper_index == lower_index + newcols:
1374            points[lower_index:upper_index, :] = tpoints
1375            elevation[lower_index:upper_index] = telev
1376
1377    assert global_index == nopoints, 'index not equal to number of points'
1378
1379    infile.close()
1380    outfile.close()
1381
1382
1383
1384def _read_hecras_cross_sections(lines):
1385    """Return block of surface lines for each cross section
1386    Starts with SURFACE LINE,
1387    Ends with END CROSS-SECTION
1388    """
1389
1390    points = []
1391
1392    reading_surface = False
1393    for i, line in enumerate(lines):
1394
1395        if len(line.strip()) == 0:    #Ignore blanks
1396            continue
1397
1398        if lines[i].strip().startswith('SURFACE LINE'):
1399            reading_surface = True
1400            continue
1401
1402        if lines[i].strip().startswith('END') and reading_surface:
1403            yield points
1404            reading_surface = False
1405            points = []
1406
1407        if reading_surface:
1408            fields = line.strip().split(',')
1409            easting = float(fields[0])
1410            northing = float(fields[1])
1411            elevation = float(fields[2])
1412            points.append([easting, northing, elevation])
1413
1414
1415
1416
1417def hecras_cross_sections2pts(basename_in,
1418                              basename_out=None,
1419                              verbose=False):
1420    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
1421
1422    Example:
1423
1424
1425# RAS export file created on Mon 15Aug2005 11:42
1426# by HEC-RAS Version 3.1.1
1427
1428BEGIN HEADER:
1429  UNITS: METRIC
1430  DTM TYPE: TIN
1431  DTM: v:\1\cit\perth_topo\river_tin
1432  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
1433  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
1434  MAP PROJECTION: UTM
1435  PROJECTION ZONE: 50
1436  DATUM: AGD66
1437  VERTICAL DATUM:
1438  NUMBER OF REACHES:  19
1439  NUMBER OF CROSS-SECTIONS:  14206
1440END HEADER:
1441
1442
1443Only the SURFACE LINE data of the following form will be utilised
1444
1445  CROSS-SECTION:
1446    STREAM ID:Southern-Wungong
1447    REACH ID:Southern-Wungong
1448    STATION:19040.*
1449    CUT LINE:
1450      405548.671603161 , 6438142.7594925
1451      405734.536092045 , 6438326.10404912
1452      405745.130459356 , 6438331.48627354
1453      405813.89633823 , 6438368.6272789
1454    SURFACE LINE:
1455     405548.67,   6438142.76,   35.37
1456     405552.24,   6438146.28,   35.41
1457     405554.78,   6438148.78,   35.44
1458     405555.80,   6438149.79,   35.44
1459     405559.37,   6438153.31,   35.45
1460     405560.88,   6438154.81,   35.44
1461     405562.93,   6438156.83,   35.42
1462     405566.50,   6438160.35,   35.38
1463     405566.99,   6438160.83,   35.37
1464     ...
1465   END CROSS-SECTION
1466
1467    Convert to NetCDF pts format which is
1468
1469    points:  (Nx2) Float array
1470    elevation: N Float array
1471    """
1472
1473    #FIXME: Can this be written feasibly using write_pts?
1474
1475    import os
1476    from Scientific.IO.NetCDF import NetCDFFile
1477    from Numeric import Float, zeros, reshape
1478
1479    root = basename_in
1480
1481    #Get ASCII file
1482    infile = open(root + '.sdf', 'r')  #Open SDF file for read
1483
1484    if verbose: print 'Reading DEM from %s' %(root + '.sdf')
1485
1486    lines = infile.readlines()
1487    infile.close()
1488
1489    if verbose: print 'Converting to pts format'
1490
1491    i = 0
1492    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
1493        i += 1
1494
1495    assert lines[i].strip().upper() == 'BEGIN HEADER:'
1496    i += 1
1497
1498    assert lines[i].strip().upper().startswith('UNITS:')
1499    units = lines[i].strip().split()[1]
1500    i += 1
1501
1502    assert lines[i].strip().upper().startswith('DTM TYPE:')
1503    i += 1
1504
1505    assert lines[i].strip().upper().startswith('DTM:')
1506    i += 1
1507
1508    assert lines[i].strip().upper().startswith('STREAM')
1509    i += 1
1510
1511    assert lines[i].strip().upper().startswith('CROSS')
1512    i += 1
1513
1514    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
1515    projection = lines[i].strip().split(':')[1]
1516    i += 1
1517
1518    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
1519    zone = int(lines[i].strip().split(':')[1])
1520    i += 1
1521
1522    assert lines[i].strip().upper().startswith('DATUM:')
1523    datum = lines[i].strip().split(':')[1]
1524    i += 1
1525
1526    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
1527    i += 1
1528
1529    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
1530    i += 1
1531
1532    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
1533    number_of_cross_sections = int(lines[i].strip().split(':')[1])
1534    i += 1
1535
1536
1537    #Now read all points
1538    points = []
1539    elevation = []
1540    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
1541        for k, entry in enumerate(entries):
1542            points.append(entry[:2])
1543            elevation.append(entry[2])
1544
1545
1546    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
1547          %(j+1, number_of_cross_sections)
1548    assert j+1 == number_of_cross_sections, msg
1549
1550    #Get output file
1551    if basename_out == None:
1552        ptsname = root + '.pts'
1553    else:
1554        ptsname = basename_out + '.pts'
1555
1556    #FIXME (DSG-ON): use loadASCII export_points_file
1557    if verbose: print 'Store to NetCDF file %s' %ptsname
1558    # NetCDF file definition
1559    outfile = NetCDFFile(ptsname, 'w')
1560
1561    #Create new file
1562    outfile.institution = 'Geoscience Australia'
1563    outfile.description = 'NetCDF pts format for compact and portable ' +\
1564                          'storage of spatial point data derived from HEC-RAS'
1565
1566    #Georeferencing
1567    outfile.zone = zone
1568    outfile.xllcorner = 0.0
1569    outfile.yllcorner = 0.0
1570    outfile.false_easting = 500000    #FIXME: Use defaults from e.g. config.py
1571    outfile.false_northing = 1000000
1572
1573    outfile.projection = projection
1574    outfile.datum = datum
1575    outfile.units = units
1576
1577
1578    outfile.createDimension('number_of_points', len(points))
1579    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1580
1581    # variable definitions
1582    outfile.createVariable('points', Float, ('number_of_points',
1583                                             'number_of_dimensions'))
1584    outfile.createVariable('elevation', Float, ('number_of_points',))
1585
1586    # Get handles to the variables
1587    outfile.variables['points'][:] = points
1588    outfile.variables['elevation'][:] = elevation
1589
1590    outfile.close()
1591
1592
1593
1594def sww2dem(basename_in, basename_out = None,
1595            quantity = None,
1596            timestep = None,
1597            reduction = None,
1598            cellsize = 10,
1599            NODATA_value = -9999,
1600            easting_min = None,
1601            easting_max = None,
1602            northing_min = None,
1603            northing_max = None,
1604            verbose = False,
1605            origin = None,
1606            datum = 'WGS84',
1607        format = 'ers'):
1608
1609    """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
1610
1611    Example (ASC):
1612
1613    ncols         3121
1614    nrows         1800
1615    xllcorner     722000
1616    yllcorner     5893000
1617    cellsize      25
1618    NODATA_value  -9999
1619    138.3698 137.4194 136.5062 135.5558 ..........
1620
1621    Also write accompanying file with same basename_in but extension .prj
1622    used to fix the UTM zone, datum, false northings and eastings.
1623
1624    The prj format is assumed to be as
1625
1626    Projection    UTM
1627    Zone          56
1628    Datum         WGS84
1629    Zunits        NO
1630    Units         METERS
1631    Spheroid      WGS84
1632    Xshift        0.0000000000
1633    Yshift        10000000.0000000000
1634    Parameters
1635
1636
1637    The parameter quantity must be the name of an existing quantity or
1638    an expression involving existing quantities. The default is
1639    'elevation'.
1640
1641    if timestep (an index) is given, output quantity at that timestep
1642
1643    if reduction is given use that to reduce quantity over all timesteps.
1644
1645    datum
1646
1647    format can be either 'asc' or 'ers'
1648    """
1649
1650    import sys
1651    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
1652    from Numeric import array2string
1653
1654    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
1655    from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
1656
1657    msg = 'Format must be either asc or ers'
1658    assert format.lower() in ['asc', 'ers'], msg
1659
1660
1661    false_easting = 500000
1662    false_northing = 10000000
1663
1664    if quantity is None:
1665        quantity = 'elevation'
1666
1667    if reduction is None:
1668        reduction = max
1669
1670    if basename_out is None:
1671        basename_out = basename_in + '_%s' %quantity
1672
1673    swwfile = basename_in + '.sww'
1674    demfile = basename_out + '.' + format
1675    # Note the use of a .ers extension is optional (write_ermapper_grid will
1676    # deal with either option
1677   
1678    #if verbose: bye= nsuadsfd[0] #uncomment to check catching verbose errors
1679   
1680    #Read sww file
1681    if verbose: print 'Reading from %s' %swwfile
1682    from Scientific.IO.NetCDF import NetCDFFile
1683    fid = NetCDFFile(swwfile)
1684
1685    #Get extent and reference
1686    x = fid.variables['x'][:]
1687    y = fid.variables['y'][:]
1688    volumes = fid.variables['volumes'][:]
1689    if timestep is not None:
1690        times = fid.variables['time'][timestep]
1691    else:
1692        times = fid.variables['time'][:]
1693
1694    number_of_timesteps = fid.dimensions['number_of_timesteps']
1695    number_of_points = fid.dimensions['number_of_points']
1696   
1697    if origin is None:
1698
1699        #Get geo_reference
1700        #sww files don't have to have a geo_ref
1701        try:
1702            geo_reference = Geo_reference(NetCDFObject=fid)
1703        except AttributeError, e:
1704            geo_reference = Geo_reference() #Default georef object
1705
1706        xllcorner = geo_reference.get_xllcorner()
1707        yllcorner = geo_reference.get_yllcorner()
1708        zone = geo_reference.get_zone()
1709    else:
1710        zone = origin[0]
1711        xllcorner = origin[1]
1712        yllcorner = origin[2]
1713
1714
1715
1716    #FIXME: Refactor using code from file_function.statistics
1717    #Something like print swwstats(swwname)
1718    if verbose:
1719        print '------------------------------------------------'
1720        print 'Statistics of SWW file:'
1721        print '  Name: %s' %swwfile
1722        print '  Reference:'
1723        print '    Lower left corner: [%f, %f]'\
1724              %(xllcorner, yllcorner)
1725        if timestep is not None:
1726            print '    Time: %f' %(times)
1727        else:
1728            print '    Start time: %f' %fid.starttime[0]
1729        print '  Extent:'
1730        print '    x [m] in [%f, %f], len(x) == %d'\
1731              %(min(x.flat), max(x.flat), len(x.flat))
1732        print '    y [m] in [%f, %f], len(y) == %d'\
1733              %(min(y.flat), max(y.flat), len(y.flat))
1734        if timestep is not None:
1735            print '    t [s] = %f, len(t) == %d' %(times, 1)
1736        else:
1737            print '    t [s] in [%f, %f], len(t) == %d'\
1738              %(min(times), max(times), len(times))
1739        print '  Quantities [SI units]:'
1740        for name in ['stage', 'xmomentum', 'ymomentum']:
1741            q = fid.variables[name][:].flat
1742            if timestep is not None:
1743                q = q[timestep*len(x):(timestep+1)*len(x)]
1744            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1745        for name in ['elevation']:
1746            q = fid.variables[name][:].flat
1747            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
1748           
1749    # Get quantity and reduce if applicable
1750    if verbose: print 'Processing quantity %s' %quantity
1751
1752    # Turn NetCDF objects into Numeric arrays
1753    quantity_dict = {}
1754    for name in fid.variables.keys():
1755        quantity_dict[name] = fid.variables[name][:] 
1756
1757
1758    # Convert quantity expression to quantities found in sww file   
1759    q = apply_expression_to_dictionary(quantity, quantity_dict)
1760
1761    if len(q.shape) == 2:
1762        #q has a time component and needs to be reduced along
1763        #the temporal dimension
1764        if verbose: print 'Reducing quantity %s' %quantity
1765        q_reduced = zeros( number_of_points, Float )
1766       
1767        if timestep is not None:
1768            for k in range(number_of_points):
1769                q_reduced[k] = q[timestep,k] 
1770        else:
1771            for k in range(number_of_points):
1772                q_reduced[k] = reduction( q[:,k] )
1773
1774        q = q_reduced
1775
1776    #Post condition: Now q has dimension: number_of_points
1777    assert len(q.shape) == 1
1778    assert q.shape[0] == number_of_points
1779
1780
1781    if verbose:
1782        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
1783
1784
1785    #Create grid and update xll/yll corner and x,y
1786
1787    #Relative extent
1788    if easting_min is None:
1789        xmin = min(x)
1790    else:
1791        xmin = easting_min - xllcorner
1792
1793    if easting_max is None:
1794        xmax = max(x)
1795    else:
1796        xmax = easting_max - xllcorner
1797
1798    if northing_min is None:
1799        ymin = min(y)
1800    else:
1801        ymin = northing_min - yllcorner
1802
1803    if northing_max is None:
1804        ymax = max(y)
1805    else:
1806        ymax = northing_max - yllcorner
1807
1808
1809
1810    if verbose: print 'Creating grid'
1811    ncols = int((xmax-xmin)/cellsize)+1
1812    nrows = int((ymax-ymin)/cellsize)+1
1813
1814
1815    #New absolute reference and coordinates
1816    newxllcorner = xmin+xllcorner
1817    newyllcorner = ymin+yllcorner
1818
1819    x = x+xllcorner-newxllcorner
1820    y = y+yllcorner-newyllcorner
1821   
1822    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
1823    assert len(vertex_points.shape) == 2
1824
1825    grid_points = zeros ( (ncols*nrows, 2), Float )
1826
1827
1828    for i in xrange(nrows):
1829        if format.lower() == 'asc':
1830            yg = i*cellsize
1831        else:
1832        #this will flip the order of the y values for ers
1833            yg = (nrows-i)*cellsize
1834
1835        for j in xrange(ncols):
1836            xg = j*cellsize
1837            k = i*ncols + j
1838
1839            grid_points[k,0] = xg
1840            grid_points[k,1] = yg
1841
1842    #Interpolate
1843    from anuga.fit_interpolate.interpolate import Interpolate
1844
1845    interp = Interpolate(vertex_points, volumes, verbose = verbose)
1846
1847    #Interpolate using quantity values
1848    if verbose: print 'Interpolating'
1849    grid_values = interp.interpolate(q, grid_points).flat
1850
1851
1852    if verbose:
1853        print 'Interpolated values are in [%f, %f]' %(min(grid_values),
1854                                                      max(grid_values))
1855
1856    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
1857    P = interp.mesh.get_boundary_polygon()
1858    outside_indices = outside_polygon(grid_points, P, closed=True)
1859
1860    for i in outside_indices:
1861        grid_values[i] = NODATA_value
1862
1863
1864
1865
1866    if format.lower() == 'ers':
1867        # setup ERS header information
1868        grid_values = reshape(grid_values,(nrows, ncols))
1869        header = {}
1870        header['datum'] = '"' + datum + '"'
1871        # FIXME The use of hardwired UTM and zone number needs to be made optional
1872        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
1873        header['projection'] = '"UTM-' + str(zone) + '"'
1874        header['coordinatetype'] = 'EN'
1875        if header['coordinatetype'] == 'LL':
1876            header['longitude'] = str(newxllcorner)
1877            header['latitude'] = str(newyllcorner)
1878        elif header['coordinatetype'] == 'EN':
1879            header['eastings'] = str(newxllcorner)
1880            header['northings'] = str(newyllcorner)
1881        header['nullcellvalue'] = str(NODATA_value)
1882        header['xdimension'] = str(cellsize)
1883        header['ydimension'] = str(cellsize)
1884        header['value'] = '"' + quantity + '"'
1885        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
1886
1887
1888        #Write
1889        if verbose: print 'Writing %s' %demfile
1890        import ermapper_grids
1891        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
1892
1893        fid.close()
1894    else:
1895        #Write to Ascii format
1896
1897        #Write prj file
1898        prjfile = basename_out + '.prj'
1899
1900        if verbose: print 'Writing %s' %prjfile
1901        prjid = open(prjfile, 'w')
1902        prjid.write('Projection    %s\n' %'UTM')
1903        prjid.write('Zone          %d\n' %zone)
1904        prjid.write('Datum         %s\n' %datum)
1905        prjid.write('Zunits        NO\n')
1906        prjid.write('Units         METERS\n')
1907        prjid.write('Spheroid      %s\n' %datum)
1908        prjid.write('Xshift        %d\n' %false_easting)
1909        prjid.write('Yshift        %d\n' %false_northing)
1910        prjid.write('Parameters\n')
1911        prjid.close()
1912
1913
1914
1915        if verbose: print 'Writing %s' %demfile
1916
1917        ascid = open(demfile, 'w')
1918
1919        ascid.write('ncols         %d\n' %ncols)
1920        ascid.write('nrows         %d\n' %nrows)
1921        ascid.write('xllcorner     %d\n' %newxllcorner)
1922        ascid.write('yllcorner     %d\n' %newyllcorner)
1923        ascid.write('cellsize      %f\n' %cellsize)
1924        ascid.write('NODATA_value  %d\n' %NODATA_value)
1925
1926
1927        #Get bounding polygon from mesh
1928        #P = interp.mesh.get_boundary_polygon()
1929        #inside_indices = inside_polygon(grid_points, P)
1930
1931        for i in range(nrows):
1932            if verbose and i%((nrows+10)/10)==0:
1933                print 'Doing row %d of %d' %(i, nrows)
1934
1935            base_index = (nrows-i-1)*ncols
1936
1937            slice = grid_values[base_index:base_index+ncols]
1938            s = array2string(slice, max_line_width=sys.maxint)
1939            ascid.write(s[1:-1] + '\n')
1940
1941
1942            #print
1943            #for j in range(ncols):
1944            #    index = base_index+j#
1945            #    print grid_values[index],
1946            #    ascid.write('%f ' %grid_values[index])
1947            #ascid.write('\n')
1948
1949        #Close
1950        ascid.close()
1951        fid.close()
1952
1953#Backwards compatibility
1954def sww2asc(basename_in, basename_out = None,
1955            quantity = None,
1956            timestep = None,
1957            reduction = None,
1958            cellsize = 10,
1959            verbose = False,
1960            origin = None):
1961    print 'sww2asc will soon be obsoleted - please use sww2dem'
1962    sww2dem(basename_in,
1963            basename_out = basename_out,
1964            quantity = quantity,
1965            timestep = timestep,
1966            reduction = reduction,
1967            cellsize = cellsize,
1968            verbose = verbose,
1969            origin = origin,
1970        datum = 'WGS84',
1971        format = 'asc')
1972
1973def sww2ers(basename_in, basename_out = None,
1974            quantity = None,
1975            timestep = None,
1976            reduction = None,
1977            cellsize = 10,
1978            verbose = False,
1979            origin = None,
1980            datum = 'WGS84'):
1981    print 'sww2ers will soon be obsoleted - please use sww2dem'
1982    sww2dem(basename_in,
1983            basename_out = basename_out,
1984            quantity = quantity,
1985            timestep = timestep,
1986            reduction = reduction,
1987            cellsize = cellsize,
1988            verbose = verbose,
1989            origin = origin,
1990            datum = datum,
1991            format = 'ers')
1992################################# END COMPATIBILITY ##############
1993
1994
1995
1996def sww2pts(basename_in, basename_out=None,
1997            data_points=None,
1998            quantity=None,
1999            timestep=None,
2000            reduction=None,
2001            NODATA_value=-9999,
2002            verbose=False,
2003            origin=None):
2004            #datum = 'WGS84')
2005
2006
2007    """Read SWW file and convert to interpolated values at selected points
2008
2009    The parameter quantity' must be the name of an existing quantity or
2010    an expression involving existing quantities. The default is
2011    'elevation'.
2012
2013    if timestep (an index) is given, output quantity at that timestep
2014
2015    if reduction is given use that to reduce quantity over all timesteps.
2016
2017    data_points (Nx2 array) give locations of points where quantity is to be computed.
2018   
2019    """
2020
2021    import sys
2022    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
2023    from Numeric import array2string
2024
2025    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
2026    from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
2027
2028    from anuga.geospatial_data.geospatial_data import Geospatial_data
2029
2030    if quantity is None:
2031        quantity = 'elevation'
2032
2033    if reduction is None:
2034        reduction = max
2035
2036    if basename_out is None:
2037        basename_out = basename_in + '_%s' %quantity
2038
2039    swwfile = basename_in + '.sww'
2040    ptsfile = basename_out + '.pts'
2041
2042    # Read sww file
2043    if verbose: print 'Reading from %s' %swwfile
2044    from Scientific.IO.NetCDF import NetCDFFile
2045    fid = NetCDFFile(swwfile)
2046
2047    # Get extent and reference
2048    x = fid.variables['x'][:]
2049    y = fid.variables['y'][:]
2050    volumes = fid.variables['volumes'][:]
2051
2052    number_of_timesteps = fid.dimensions['number_of_timesteps']
2053    number_of_points = fid.dimensions['number_of_points']
2054    if origin is None:
2055
2056        # Get geo_reference
2057        # sww files don't have to have a geo_ref
2058        try:
2059            geo_reference = Geo_reference(NetCDFObject=fid)
2060        except AttributeError, e:
2061            geo_reference = Geo_reference() #Default georef object
2062
2063        xllcorner = geo_reference.get_xllcorner()
2064        yllcorner = geo_reference.get_yllcorner()
2065        zone = geo_reference.get_zone()
2066    else:
2067        zone = origin[0]
2068        xllcorner = origin[1]
2069        yllcorner = origin[2]
2070
2071
2072
2073    # FIXME: Refactor using code from file_function.statistics
2074    # Something like print swwstats(swwname)
2075    if verbose:
2076        x = fid.variables['x'][:]
2077        y = fid.variables['y'][:]
2078        times = fid.variables['time'][:]
2079        print '------------------------------------------------'
2080        print 'Statistics of SWW file:'
2081        print '  Name: %s' %swwfile
2082        print '  Reference:'
2083        print '    Lower left corner: [%f, %f]'\
2084              %(xllcorner, yllcorner)
2085        print '    Start time: %f' %fid.starttime[0]
2086        print '  Extent:'
2087        print '    x [m] in [%f, %f], len(x) == %d'\
2088              %(min(x.flat), max(x.flat), len(x.flat))
2089        print '    y [m] in [%f, %f], len(y) == %d'\
2090              %(min(y.flat), max(y.flat), len(y.flat))
2091        print '    t [s] in [%f, %f], len(t) == %d'\
2092              %(min(times), max(times), len(times))
2093        print '  Quantities [SI units]:'
2094        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2095            q = fid.variables[name][:].flat
2096            print '    %s in [%f, %f]' %(name, min(q), max(q))
2097
2098
2099
2100    # Get quantity and reduce if applicable
2101    if verbose: print 'Processing quantity %s' %quantity
2102
2103    # Turn NetCDF objects into Numeric arrays
2104    quantity_dict = {}
2105    for name in fid.variables.keys():
2106        quantity_dict[name] = fid.variables[name][:]
2107
2108
2109
2110    # Convert quantity expression to quantities found in sww file   
2111    q = apply_expression_to_dictionary(quantity, quantity_dict)
2112
2113
2114
2115    if len(q.shape) == 2:
2116        # q has a time component and needs to be reduced along
2117        # the temporal dimension
2118        if verbose: print 'Reducing quantity %s' %quantity
2119        q_reduced = zeros( number_of_points, Float )
2120
2121        for k in range(number_of_points):
2122            q_reduced[k] = reduction( q[:,k] )
2123
2124        q = q_reduced
2125
2126    # Post condition: Now q has dimension: number_of_points
2127    assert len(q.shape) == 1
2128    assert q.shape[0] == number_of_points
2129
2130
2131    if verbose:
2132        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
2133
2134
2135    # Create grid and update xll/yll corner and x,y
2136    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
2137    assert len(vertex_points.shape) == 2
2138
2139    # Interpolate
2140    from anuga.fit_interpolate.interpolate import Interpolate
2141    interp = Interpolate(vertex_points, volumes, verbose = verbose)
2142
2143    # Interpolate using quantity values
2144    if verbose: print 'Interpolating'
2145    interpolated_values = interp.interpolate(q, data_points).flat
2146
2147
2148    if verbose:
2149        print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
2150                                                      max(interpolated_values))
2151
2152    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
2153    P = interp.mesh.get_boundary_polygon()
2154    outside_indices = outside_polygon(data_points, P, closed=True)
2155
2156    for i in outside_indices:
2157        interpolated_values[i] = NODATA_value
2158
2159
2160    # Store results   
2161    G = Geospatial_data(data_points=data_points,
2162                        attributes=interpolated_values)
2163
2164    G.export_points_file(ptsfile, absolute = True)
2165
2166    fid.close()
2167
2168
2169def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2170                                  use_cache = False,
2171                                  verbose = False):
2172    """Read Digitial Elevation model from the following ASCII format (.asc)
2173
2174    Example:
2175
2176    ncols         3121
2177    nrows         1800
2178    xllcorner     722000
2179    yllcorner     5893000
2180    cellsize      25
2181    NODATA_value  -9999
2182    138.3698 137.4194 136.5062 135.5558 ..........
2183
2184    Convert basename_in + '.asc' to NetCDF format (.dem)
2185    mimicking the ASCII format closely.
2186
2187
2188    An accompanying file with same basename_in but extension .prj must exist
2189    and is used to fix the UTM zone, datum, false northings and eastings.
2190
2191    The prj format is assumed to be as
2192
2193    Projection    UTM
2194    Zone          56
2195    Datum         WGS84
2196    Zunits        NO
2197    Units         METERS
2198    Spheroid      WGS84
2199    Xshift        0.0000000000
2200    Yshift        10000000.0000000000
2201    Parameters
2202    """
2203
2204
2205
2206    kwargs = {'basename_out': basename_out, 'verbose': verbose}
2207
2208    if use_cache is True:
2209        from caching import cache
2210        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
2211                       dependencies = [basename_in + '.asc'],
2212                       verbose = verbose)
2213
2214    else:
2215        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
2216
2217    return result
2218
2219
2220
2221
2222
2223
2224def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
2225                                  verbose = False):
2226    """Read Digitial Elevation model from the following ASCII format (.asc)
2227
2228    Internal function. See public function convert_dem_from_ascii2netcdf for details.
2229    """
2230
2231    import os
2232    from Scientific.IO.NetCDF import NetCDFFile
2233    from Numeric import Float, array
2234
2235    #root, ext = os.path.splitext(basename_in)
2236    root = basename_in
2237
2238    ###########################################
2239    # Read Meta data
2240    if verbose: print 'Reading METADATA from %s' %root + '.prj'
2241    metadatafile = open(root + '.prj')
2242    metalines = metadatafile.readlines()
2243    metadatafile.close()
2244
2245    L = metalines[0].strip().split()
2246    assert L[0].strip().lower() == 'projection'
2247    projection = L[1].strip()                   #TEXT
2248
2249    L = metalines[1].strip().split()
2250    assert L[0].strip().lower() == 'zone'
2251    zone = int(L[1].strip())
2252
2253    L = metalines[2].strip().split()
2254    assert L[0].strip().lower() == 'datum'
2255    datum = L[1].strip()                        #TEXT
2256
2257    L = metalines[3].strip().split()
2258    assert L[0].strip().lower() == 'zunits'     #IGNORE
2259    zunits = L[1].strip()                       #TEXT
2260
2261    L = metalines[4].strip().split()
2262    assert L[0].strip().lower() == 'units'
2263    units = L[1].strip()                        #TEXT
2264
2265    L = metalines[5].strip().split()
2266    assert L[0].strip().lower() == 'spheroid'   #IGNORE
2267    spheroid = L[1].strip()                     #TEXT
2268
2269    L = metalines[6].strip().split()
2270    assert L[0].strip().lower() == 'xshift'
2271    false_easting = float(L[1].strip())
2272
2273    L = metalines[7].strip().split()
2274    assert L[0].strip().lower() == 'yshift'
2275    false_northing = float(L[1].strip())
2276
2277    #print false_easting, false_northing, zone, datum
2278
2279
2280    ###########################################
2281    #Read DEM data
2282
2283    datafile = open(basename_in + '.asc')
2284
2285    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
2286    lines = datafile.readlines()
2287    datafile.close()
2288
2289    if verbose: print 'Got', len(lines), ' lines'
2290
2291    ncols = int(lines[0].split()[1].strip())
2292    nrows = int(lines[1].split()[1].strip())
2293    xllcorner = float(lines[2].split()[1].strip())
2294    yllcorner = float(lines[3].split()[1].strip())
2295    cellsize = float(lines[4].split()[1].strip())
2296    NODATA_value = int(lines[5].split()[1].strip())
2297
2298    assert len(lines) == nrows + 6
2299
2300
2301    ##########################################
2302
2303
2304    if basename_out == None:
2305        netcdfname = root + '.dem'
2306    else:
2307        netcdfname = basename_out + '.dem'
2308
2309    if verbose: print 'Store to NetCDF file %s' %netcdfname
2310    # NetCDF file definition
2311    fid = NetCDFFile(netcdfname, 'w')
2312
2313    #Create new file
2314    fid.institution = 'Geoscience Australia'
2315    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
2316                      'of spatial point data'
2317
2318    fid.ncols = ncols
2319    fid.nrows = nrows
2320    fid.xllcorner = xllcorner
2321    fid.yllcorner = yllcorner
2322    fid.cellsize = cellsize
2323    fid.NODATA_value = NODATA_value
2324
2325    fid.zone = zone
2326    fid.false_easting = false_easting
2327    fid.false_northing = false_northing
2328    fid.projection = projection
2329    fid.datum = datum
2330    fid.units = units
2331
2332
2333    # dimension definitions
2334    fid.createDimension('number_of_rows', nrows)
2335    fid.createDimension('number_of_columns', ncols)
2336
2337    # variable definitions
2338    fid.createVariable('elevation', Float, ('number_of_rows',
2339                                            'number_of_columns'))
2340
2341    # Get handles to the variables
2342    elevation = fid.variables['elevation']
2343
2344    #Store data
2345    n = len(lines[6:])
2346    for i, line in enumerate(lines[6:]):
2347        fields = line.split()
2348        if verbose and i%((n+10)/10)==0:
2349            print 'Processing row %d of %d' %(i, nrows)
2350
2351        elevation[i, :] = array([float(x) for x in fields])
2352
2353    fid.close()
2354
2355
2356
2357
2358
2359def ferret2sww(basename_in, basename_out = None,
2360               verbose = False,
2361               minlat = None, maxlat = None,
2362               minlon = None, maxlon = None,
2363               mint = None, maxt = None, mean_stage = 0,
2364               origin = None, zscale = 1,
2365               fail_on_NaN = True,
2366               NaN_filler = 0,
2367               elevation = None,
2368               inverted_bathymetry = True
2369               ): #FIXME: Bathymetry should be obtained
2370                                  #from MOST somehow.
2371                                  #Alternatively from elsewhere
2372                                  #or, as a last resort,
2373                                  #specified here.
2374                                  #The value of -100 will work
2375                                  #for the Wollongong tsunami
2376                                  #scenario but is very hacky
2377    """Convert MOST and 'Ferret' NetCDF format for wave propagation to
2378    sww format native to abstract_2d_finite_volumes.
2379
2380    Specify only basename_in and read files of the form
2381    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
2382    relative height, x-velocity and y-velocity, respectively.
2383
2384    Also convert latitude and longitude to UTM. All coordinates are
2385    assumed to be given in the GDA94 datum.
2386
2387    min's and max's: If omitted - full extend is used.
2388    To include a value min may equal it, while max must exceed it.
2389    Lat and lon are assuemd to be in decimal degrees
2390
2391    origin is a 3-tuple with geo referenced
2392    UTM coordinates (zone, easting, northing)
2393
2394    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
2395    which means that longitude is the fastest
2396    varying dimension (row major order, so to speak)
2397
2398    ferret2sww uses grid points as vertices in a triangular grid
2399    counting vertices from lower left corner upwards, then right
2400    """
2401
2402    import os
2403    from Scientific.IO.NetCDF import NetCDFFile
2404    from Numeric import Float, Int, Int32, searchsorted, zeros, array
2405    from Numeric import allclose, around
2406
2407    precision = Float
2408
2409    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
2410
2411    if minlat != None:
2412        assert -90 < minlat < 90 , msg
2413    if maxlat != None:
2414        assert -90 < maxlat < 90 , msg
2415        if minlat != None:
2416            assert maxlat > minlat
2417    if minlon != None:
2418        assert -180 < minlon < 180 , msg
2419    if maxlon != None:
2420        assert -180 < maxlon < 180 , msg
2421        if minlon != None:
2422            assert maxlon > minlon
2423       
2424
2425
2426    #Get NetCDF data
2427    if verbose: print 'Reading files %s_*.nc' %basename_in
2428    #print "basename_in + '_ha.nc'",basename_in + '_ha.nc'
2429    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
2430    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
2431    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
2432    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
2433
2434    if basename_out is None:
2435        swwname = basename_in + '.sww'
2436    else:
2437        swwname = basename_out + '.sww'
2438
2439#    Get dimensions of file_h
2440    for dimension in file_h.dimensions.keys():
2441        if dimension[:3] == 'LON':
2442            dim_h_longitude = dimension
2443        if dimension[:3] == 'LAT':
2444            dim_h_latitude = dimension
2445        if dimension[:4] == 'TIME':
2446            dim_h_time = dimension
2447
2448#    print 'long:', dim_h_longitude
2449#    print 'lats:', dim_h_latitude
2450#    print 'times:', dim_h_time
2451
2452    times = file_h.variables[dim_h_time]
2453    latitudes = file_h.variables[dim_h_latitude]
2454    longitudes = file_h.variables[dim_h_longitude]
2455   
2456# get dimensions for file_e
2457    for dimension in file_e.dimensions.keys():
2458        if dimension[:3] == 'LON':
2459            dim_e_longitude = dimension
2460        if dimension[:3] == 'LAT':
2461            dim_e_latitude = dimension
2462
2463# get dimensions for file_u
2464    for dimension in file_u.dimensions.keys():
2465        if dimension[:3] == 'LON':
2466            dim_u_longitude = dimension
2467        if dimension[:3] == 'LAT':
2468            dim_u_latitude = dimension
2469        if dimension[:4] == 'TIME':
2470            dim_u_time = dimension
2471           
2472# get dimensions for file_v
2473    for dimension in file_v.dimensions.keys():
2474        if dimension[:3] == 'LON':
2475            dim_v_longitude = dimension
2476        if dimension[:3] == 'LAT':
2477            dim_v_latitude = dimension
2478        if dimension[:4] == 'TIME':
2479            dim_v_time = dimension
2480
2481
2482    #Precision used by most for lat/lon is 4 or 5 decimals
2483    e_lat = around(file_e.variables[dim_e_latitude][:], 5)
2484    e_lon = around(file_e.variables[dim_e_longitude][:], 5)
2485
2486    #Check that files are compatible
2487    assert allclose(latitudes, file_u.variables[dim_u_latitude])
2488    assert allclose(latitudes, file_v.variables[dim_v_latitude])
2489    assert allclose(latitudes, e_lat)
2490
2491    assert allclose(longitudes, file_u.variables[dim_u_longitude])
2492    assert allclose(longitudes, file_v.variables[dim_v_longitude])
2493    assert allclose(longitudes, e_lon)
2494
2495    if mint == None:
2496        jmin = 0
2497    else:
2498        jmin = searchsorted(times, mint)
2499       
2500    if maxt == None:
2501        jmax=len(times)
2502    else:
2503        jmax = searchsorted(times, maxt)
2504
2505    #print "latitudes[:]",latitudes[:]
2506    #print "longitudes[:]",longitudes [:]
2507    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
2508                                                  longitudes[:],
2509                                                 minlat, maxlat,
2510                                                 minlon, maxlon)
2511
2512
2513#    print' j', jmin, jmax
2514    times = times[jmin:jmax]
2515    latitudes = latitudes[kmin:kmax]
2516    longitudes = longitudes[lmin:lmax]
2517
2518    #print "latitudes[:]",latitudes[:]
2519    #print "longitudes[:]",longitudes [:]
2520
2521    if verbose: print 'cropping'
2522    zname = 'ELEVATION'
2523
2524    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
2525    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
2526    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
2527    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
2528
2529    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
2530    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
2531    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
2532    #        from Numeric import asarray
2533    #        elevations=elevations.tolist()
2534    #        elevations.reverse()
2535    #        elevations=asarray(elevations)
2536    #    else:
2537    #        from Numeric import asarray
2538    #        elevations=elevations.tolist()
2539    #        elevations.reverse()
2540    #        elevations=asarray(elevations)
2541    #        'print hmmm'
2542
2543
2544
2545    #Get missing values
2546    nan_ha = file_h.variables['HA'].missing_value[0]
2547    nan_ua = file_u.variables['UA'].missing_value[0]
2548    nan_va = file_v.variables['VA'].missing_value[0]
2549    if hasattr(file_e.variables[zname],'missing_value'):
2550        nan_e  = file_e.variables[zname].missing_value[0]
2551    else:
2552        nan_e = None
2553
2554    #Cleanup
2555    from Numeric import sometrue
2556
2557    missing = (amplitudes == nan_ha)
2558    if sometrue (missing):
2559        if fail_on_NaN:
2560            msg = 'NetCDFFile %s contains missing values'\
2561                  %(basename_in+'_ha.nc')
2562            raise DataMissingValuesError, msg
2563        else:
2564            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
2565
2566    missing = (uspeed == nan_ua)
2567    if sometrue (missing):
2568        if fail_on_NaN:
2569            msg = 'NetCDFFile %s contains missing values'\
2570                  %(basename_in+'_ua.nc')
2571            raise DataMissingValuesError, msg
2572        else:
2573            uspeed = uspeed*(missing==0) + missing*NaN_filler
2574
2575    missing = (vspeed == nan_va)
2576    if sometrue (missing):
2577        if fail_on_NaN:
2578            msg = 'NetCDFFile %s contains missing values'\
2579                  %(basename_in+'_va.nc')
2580            raise DataMissingValuesError, msg
2581        else:
2582            vspeed = vspeed*(missing==0) + missing*NaN_filler
2583
2584
2585    missing = (elevations == nan_e)
2586    if sometrue (missing):
2587        if fail_on_NaN:
2588            msg = 'NetCDFFile %s contains missing values'\
2589                  %(basename_in+'_e.nc')
2590            raise DataMissingValuesError, msg
2591        else:
2592            elevations = elevations*(missing==0) + missing*NaN_filler
2593
2594    #######
2595
2596
2597
2598    number_of_times = times.shape[0]
2599    number_of_latitudes = latitudes.shape[0]
2600    number_of_longitudes = longitudes.shape[0]
2601
2602    assert amplitudes.shape[0] == number_of_times
2603    assert amplitudes.shape[1] == number_of_latitudes
2604    assert amplitudes.shape[2] == number_of_longitudes
2605
2606    if verbose:
2607        print '------------------------------------------------'
2608        print 'Statistics:'
2609        print '  Extent (lat/lon):'
2610        print '    lat in [%f, %f], len(lat) == %d'\
2611              %(min(latitudes.flat), max(latitudes.flat),
2612                len(latitudes.flat))
2613        print '    lon in [%f, %f], len(lon) == %d'\
2614              %(min(longitudes.flat), max(longitudes.flat),
2615                len(longitudes.flat))
2616        print '    t in [%f, %f], len(t) == %d'\
2617              %(min(times.flat), max(times.flat), len(times.flat))
2618
2619        q = amplitudes.flat
2620        name = 'Amplitudes (ha) [cm]'
2621        print %s in [%f, %f]' %(name, min(q), max(q))
2622
2623        q = uspeed.flat
2624        name = 'Speeds (ua) [cm/s]'
2625        print %s in [%f, %f]' %(name, min(q), max(q))
2626
2627        q = vspeed.flat
2628        name = 'Speeds (va) [cm/s]'
2629        print %s in [%f, %f]' %(name, min(q), max(q))
2630
2631        q = elevations.flat
2632        name = 'Elevations (e) [m]'
2633        print %s in [%f, %f]' %(name, min(q), max(q))
2634
2635
2636    #print number_of_latitudes, number_of_longitudes
2637    number_of_points = number_of_latitudes*number_of_longitudes
2638    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
2639
2640
2641    file_h.close()
2642    file_u.close()
2643    file_v.close()
2644    file_e.close()
2645
2646
2647    # NetCDF file definition
2648    outfile = NetCDFFile(swwname, 'w')
2649
2650    description = 'Converted from Ferret files: %s, %s, %s, %s'\
2651                  %(basename_in + '_ha.nc',
2652                    basename_in + '_ua.nc',
2653                    basename_in + '_va.nc',
2654                    basename_in + '_e.nc')
2655   
2656    #Create new file
2657    starttime = times[0]
2658    write_sww_header(outfile, times, number_of_volumes,
2659                     number_of_points, description=description) 
2660
2661    #Store
2662    from anuga.coordinate_transforms.redfearn import redfearn
2663    x = zeros(number_of_points, Float)  #Easting
2664    y = zeros(number_of_points, Float)  #Northing
2665
2666
2667    if verbose: print 'Making triangular grid'
2668    #Check zone boundaries
2669    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
2670
2671    vertices = {}
2672    i = 0
2673    for k, lat in enumerate(latitudes):       #Y direction
2674        for l, lon in enumerate(longitudes):  #X direction
2675
2676            vertices[l,k] = i
2677
2678            zone, easting, northing = redfearn(lat,lon)
2679
2680            msg = 'Zone boundary crossed at longitude =', lon
2681            #assert zone == refzone, msg
2682            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
2683            x[i] = easting
2684            y[i] = northing
2685            i += 1
2686
2687    #Construct 2 triangles per 'rectangular' element
2688    volumes = []
2689    for l in range(number_of_longitudes-1):    #X direction
2690        for k in range(number_of_latitudes-1): #Y direction
2691            v1 = vertices[l,k+1]
2692            v2 = vertices[l,k]
2693            v3 = vertices[l+1,k+1]
2694            v4 = vertices[l+1,k]
2695
2696            volumes.append([v1,v2,v3]) #Upper element
2697            volumes.append([v4,v3,v2]) #Lower element
2698
2699    volumes = array(volumes)
2700
2701    if origin is None:
2702        origin = Geo_reference(refzone,min(x),min(y))
2703    geo_ref = write_NetCDF_georeference(origin, outfile)
2704   
2705    if elevation is not None:
2706        z = elevation
2707    else:
2708        if inverted_bathymetry:
2709            z = -1*elevations
2710        else:
2711            z = elevations
2712    #FIXME: z should be obtained from MOST and passed in here
2713
2714    from Numeric import resize
2715    z = resize(z,outfile.variables['z'][:].shape)
2716    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
2717    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
2718    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
2719    outfile.variables['elevation'][:] = z
2720    outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
2721
2722
2723
2724    #Time stepping
2725    stage = outfile.variables['stage']
2726    xmomentum = outfile.variables['xmomentum']
2727    ymomentum = outfile.variables['ymomentum']
2728
2729    if verbose: print 'Converting quantities'
2730    n = len(times)
2731    for j in range(n):
2732        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
2733        i = 0
2734        for k in range(number_of_latitudes):      #Y direction
2735            for l in range(number_of_longitudes): #X direction
2736                w = zscale*amplitudes[j,k,l]/100 + mean_stage
2737                stage[j,i] = w
2738                h = w - z[i]
2739                xmomentum[j,i] = uspeed[j,k,l]/100*h
2740                ymomentum[j,i] = vspeed[j,k,l]/100*h
2741                i += 1
2742
2743    #outfile.close()
2744
2745    #FIXME: Refactor using code from file_function.statistics
2746    #Something like print swwstats(swwname)
2747    if verbose:
2748        x = outfile.variables['x'][:]
2749        y = outfile.variables['y'][:]
2750        times = outfile.variables['time'][:]
2751        print '------------------------------------------------'
2752        print 'Statistics of output file:'
2753        print '  Name: %s' %swwname
2754        print '  Reference:'
2755        print '    Lower left corner: [%f, %f]'\
2756              %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
2757        print ' Start time: %f' %starttime
2758        #print '    Min time: %f' %mint
2759        #print '    Max time: %f' %maxt
2760        print '  Extent:'
2761        print '    x [m] in [%f, %f], len(x) == %d'\
2762              %(min(x.flat), max(x.flat), len(x.flat))
2763        print '    y [m] in [%f, %f], len(y) == %d'\
2764              %(min(y.flat), max(y.flat), len(y.flat))
2765        print '    t [s] in [%f, %f], len(t) == %d'\
2766              %(min(times), max(times), len(times))
2767        print '  Quantities [SI units]:'
2768        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
2769            q = outfile.variables[name][:].flat
2770            print '    %s in [%f, %f]' %(name, min(q), max(q))
2771
2772
2773
2774    outfile.close()
2775
2776
2777
2778
2779
2780def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
2781    """Template for converting typical text files with time series to
2782    NetCDF tms file.
2783
2784
2785    The file format is assumed to be either two fields separated by a comma:
2786
2787        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
2788
2789    E.g
2790
2791      31/08/04 00:00:00, 1.328223 0 0
2792      31/08/04 00:15:00, 1.292912 0 0
2793
2794    or time (seconds), value0 value1 value2 ...
2795   
2796      0.0, 1.328223 0 0
2797      0.1, 1.292912 0 0
2798
2799    will provide a time dependent function f(t) with three attributes
2800
2801    filename is assumed to be the rootname with extenisons .txt and .sww
2802    """
2803
2804    import time, calendar
2805    from Numeric import array
2806    from anuga.config import time_format
2807    from anuga.utilities.numerical_tools import ensure_numeric
2808
2809
2810
2811    fid = open(filename + '.txt')
2812    line = fid.readline()
2813    fid.close()
2814
2815    fields = line.split(',')
2816    msg = 'File %s must have the format date, value0 value1 value2 ...'
2817    assert len(fields) == 2, msg
2818
2819    if not time_as_seconds:
2820        try:
2821            starttime = calendar.timegm(time.strptime(fields[0], time_format))
2822        except ValueError:
2823            msg = 'First field in file %s must be' %filename
2824            msg += ' date-time with format %s.\n' %time_format
2825            msg += 'I got %s instead.' %fields[0]
2826            raise DataTimeError, msg
2827    else:
2828        try:
2829            starttime = float(fields[0])
2830        except Error:
2831            msg = "Bad time format"
2832            raise DataTimeError, msg
2833
2834
2835    #Split values
2836    values = []
2837    for value in fields[1].split():
2838        values.append(float(value))
2839
2840    q = ensure_numeric(values)
2841
2842    msg = 'ERROR: File must contain at least one independent value'
2843    assert len(q.shape) == 1, msg
2844
2845
2846
2847    #Read times proper
2848    from Numeric import zeros, Float, alltrue
2849    from anuga.config import time_format
2850    import time, calendar
2851
2852    fid = open(filename + '.txt')
2853    lines = fid.readlines()
2854    fid.close()
2855
2856    N = len(lines)
2857    d = len(q)
2858
2859    T = zeros(N, Float)       #Time
2860    Q = zeros((N, d), Float)  #Values
2861
2862    for i, line in enumerate(lines):
2863        fields = line.split(',')
2864        if not time_as_seconds:
2865            realtime = calendar.timegm(time.strptime(fields[0], time_format))
2866        else:
2867             realtime = float(fields[0])
2868        T[i] = realtime - starttime
2869
2870        for j, value in enumerate(fields[1].split()):
2871            Q[i, j] = float(value)
2872
2873    msg = 'File %s must list time as a monotonuosly ' %filename
2874    msg += 'increasing sequence'
2875    assert alltrue( T[1:] - T[:-1] > 0 ), msg
2876
2877    #Create NetCDF file
2878    from Scientific.IO.NetCDF import NetCDFFile
2879
2880    fid = NetCDFFile(filename + '.tms', 'w')
2881
2882
2883    fid.institution = 'Geoscience Australia'
2884    fid.description = 'Time series'
2885
2886
2887    #Reference point
2888    #Start time in seconds since the epoch (midnight 1/1/1970)
2889    #FIXME: Use Georef
2890    fid.starttime = starttime
2891
2892    # dimension definitions
2893    #fid.createDimension('number_of_volumes', self.number_of_volumes)
2894    #fid.createDimension('number_of_vertices', 3)
2895
2896
2897    fid.createDimension('number_of_timesteps', len(T))
2898
2899    fid.createVariable('time', Float, ('number_of_timesteps',))
2900
2901    fid.variables['time'][:] = T
2902
2903    for i in range(Q.shape[1]):
2904        try:
2905            name = quantity_names[i]
2906        except:
2907            name = 'Attribute%d'%i
2908
2909        fid.createVariable(name, Float, ('number_of_timesteps',))
2910        fid.variables[name][:] = Q[:,i]
2911
2912    fid.close()
2913
2914
2915def extent_sww(file_name):
2916    """
2917    Read in an sww file.
2918
2919    Input;
2920    file_name - the sww file
2921
2922    Output;
2923    z - Vector of bed elevation
2924    volumes - Array.  Each row has 3 values, representing
2925    the vertices that define the volume
2926    time - Vector of the times where there is stage information
2927    stage - array with respect to time and vertices (x,y)
2928    """
2929
2930
2931    from Scientific.IO.NetCDF import NetCDFFile
2932
2933    #Check contents
2934    #Get NetCDF
2935    fid = NetCDFFile(file_name, 'r')
2936
2937    # Get the variables
2938    x = fid.variables['x'][:]
2939    y = fid.variables['y'][:]
2940    stage = fid.variables['stage'][:]
2941    #print "stage",stage
2942    #print "stage.shap",stage.shape
2943    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
2944    #print "min(stage)",min(stage)
2945
2946    fid.close()
2947
2948    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
2949
2950
2951def sww2domain(filename,boundary=None,t=None,\
2952               fail_if_NaN=True,NaN_filler=0\
2953               ,verbose = False,very_verbose = False):
2954    """
2955    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
2956
2957    Boundary is not recommended if domain.smooth is not selected, as it
2958    uses unique coordinates, but not unique boundaries. This means that
2959    the boundary file will not be compatable with the coordinates, and will
2960    give a different final boundary, or crash.
2961    """
2962    NaN=9.969209968386869e+036
2963    #initialise NaN.
2964
2965    from Scientific.IO.NetCDF import NetCDFFile
2966    from shallow_water import Domain
2967    from Numeric import asarray, transpose, resize
2968
2969    if verbose: print 'Reading from ', filename
2970    fid = NetCDFFile(filename, 'r')    #Open existing file for read
2971    time = fid.variables['time']       #Timesteps
2972    if t is None:
2973        t = time[-1]
2974    time_interp = get_time_interp(time,t)
2975
2976    # Get the variables as Numeric arrays
2977    x = fid.variables['x'][:]             #x-coordinates of vertices
2978    y = fid.variables['y'][:]             #y-coordinates of vertices
2979    elevation = fid.variables['elevation']     #Elevation
2980    stage = fid.variables['stage']     #Water level
2981    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
2982    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
2983
2984    starttime = fid.starttime[0]
2985    volumes = fid.variables['volumes'][:] #Connectivity
2986    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
2987
2988    conserved_quantities = []
2989    interpolated_quantities = {}
2990    other_quantities = []
2991
2992    # get geo_reference
2993    #sww files don't have to have a geo_ref
2994    try:
2995        geo_reference = Geo_reference(NetCDFObject=fid)
2996    except: #AttributeError, e:
2997        geo_reference = None
2998
2999    if verbose: print '    getting quantities'
3000    for quantity in fid.variables.keys():
3001        dimensions = fid.variables[quantity].dimensions
3002        if 'number_of_timesteps' in dimensions:
3003            conserved_quantities.append(quantity)
3004            interpolated_quantities[quantity]=\
3005                  interpolated_quantity(fid.variables[quantity][:],time_interp)
3006        else: other_quantities.append(quantity)
3007
3008    other_quantities.remove('x')
3009    other_quantities.remove('y')
3010    other_quantities.remove('z')
3011    other_quantities.remove('volumes')
3012
3013    conserved_quantities.remove('time')
3014
3015    if verbose: print '    building domain'
3016    #    From domain.Domain:
3017    #    domain = Domain(coordinates, volumes,\
3018    #                    conserved_quantities = conserved_quantities,\
3019    #                    other_quantities = other_quantities,zone=zone,\
3020    #                    xllcorner=xllcorner, yllcorner=yllcorner)
3021
3022    #   From shallow_water.Domain:
3023    coordinates=coordinates.tolist()
3024    volumes=volumes.tolist()
3025    #FIXME:should this be in mesh?(peter row)
3026    if fid.smoothing == 'Yes': unique = False
3027    else: unique = True
3028    if unique:
3029        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
3030
3031
3032    try:
3033        domain = Domain(coordinates, volumes, boundary)
3034    except AssertionError, e:
3035        fid.close()
3036        msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
3037        raise DataDomainError, msg
3038
3039    if not boundary is None:
3040        domain.boundary = boundary
3041
3042    domain.geo_reference = geo_reference
3043
3044    domain.starttime=float(starttime)+float(t)
3045    domain.time=0.0
3046
3047    for quantity in other_quantities:
3048        try:
3049            NaN = fid.variables[quantity].missing_value
3050        except:
3051            pass #quantity has no missing_value number
3052        X = fid.variables[quantity][:]
3053        if very_verbose:
3054            print '       ',quantity
3055            print '        NaN =',NaN
3056            print '        max(X)'
3057            print '       ',max(X)
3058            print '        max(X)==NaN'
3059            print '       ',max(X)==NaN
3060            print ''
3061        if (max(X)==NaN) or (min(X)==NaN):
3062            if fail_if_NaN:
3063                msg = 'quantity "%s" contains no_data entry'%quantity
3064                raise DataMissingValuesError, msg
3065            else:
3066                data = (X<>NaN)
3067                X = (X*data)+(data==0)*NaN_filler
3068        if unique:
3069            X = resize(X,(len(X)/3,3))
3070        domain.set_quantity(quantity,X)
3071    #
3072    for quantity in conserved_quantities:
3073        try:
3074            NaN = fid.variables[quantity].missing_value
3075        except:
3076            pass #quantity has no missing_value number
3077        X = interpolated_quantities[quantity]
3078        if very_verbose:
3079            print '       ',quantity
3080            print '        NaN =',NaN
3081            print '        max(X)'
3082            print '       ',max(X)
3083            print '        max(X)==NaN'
3084            print '       ',max(X)==NaN
3085            print ''
3086        if (max(X)==NaN) or (min(X)==NaN):
3087            if fail_if_NaN:
3088                msg = 'quantity "%s" contains no_data entry'%quantity
3089                raise DataMissingValuesError, msg
3090            else:
3091                data = (X<>NaN)
3092                X = (X*data)+(data==0)*NaN_filler
3093        if unique:
3094            X = resize(X,(X.shape[0]/3,3))
3095        domain.set_quantity(quantity,X)
3096
3097    fid.close()
3098    return domain
3099
3100def interpolated_quantity(saved_quantity,time_interp):
3101
3102    #given an index and ratio, interpolate quantity with respect to time.
3103    index,ratio = time_interp
3104    Q = saved_quantity
3105    if ratio > 0:
3106        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
3107    else:
3108        q = Q[index]
3109    #Return vector of interpolated values
3110    return q
3111
3112def get_time_interp(time,t=None):
3113    #Finds the ratio and index for time interpolation.
3114    #It is borrowed from previous abstract_2d_finite_volumes code.
3115    if t is None:
3116        t=time[-1]
3117        index = -1
3118        ratio = 0.
3119    else:
3120        T = time
3121        tau = t
3122        index=0
3123        msg = 'Time interval derived from file %s [%s:%s]'\
3124            %('FIXMEfilename', T[0], T[-1])
3125        msg += ' does not match model time: %s' %tau
3126        if tau < time[0]: raise DataTimeError, msg
3127        if tau > time[-1]: raise DataTimeError, msg
3128        while tau > time[index]: index += 1
3129        while tau < time[index]: index -= 1
3130        if tau == time[index]:
3131            #Protect against case where tau == time[-1] (last time)
3132            # - also works in general when tau == time[i]
3133            ratio = 0
3134        else:
3135            #t is now between index and index+1
3136            ratio = (tau - time[index])/(time[index+1] - time[index])
3137    return (index,ratio)
3138
3139
3140def weed(coordinates,volumes,boundary = None):
3141    if type(coordinates)=='array':
3142        coordinates = coordinates.tolist()
3143    if type(volumes)=='array':
3144        volumes = volumes.tolist()
3145
3146    unique = False
3147    point_dict = {}
3148    same_point = {}
3149    for i in range(len(coordinates)):
3150        point = tuple(coordinates[i])
3151        if point_dict.has_key(point):
3152            unique = True
3153            same_point[i]=point
3154            #to change all point i references to point j
3155        else:
3156            point_dict[point]=i
3157            same_point[i]=point
3158
3159    coordinates = []
3160    i = 0
3161    for point in point_dict.keys():
3162        point = tuple(point)
3163        coordinates.append(list(point))
3164        point_dict[point]=i
3165        i+=1
3166
3167
3168    for volume in volumes:
3169        for i in range(len(volume)):
3170            index = volume[i]
3171            if index>-1:
3172                volume[i]=point_dict[same_point[index]]
3173
3174    new_boundary = {}
3175    if not boundary is None:
3176        for segment in boundary.keys():
3177            point0 = point_dict[same_point[segment[0]]]
3178            point1 = point_dict[same_point[segment[1]]]
3179            label = boundary[segment]
3180            #FIXME should the bounday attributes be concaterated
3181            #('exterior, pond') or replaced ('pond')(peter row)
3182
3183            if new_boundary.has_key((point0,point1)):
3184                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
3185                                              #+','+label
3186
3187            elif new_boundary.has_key((point1,point0)):
3188                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
3189                                              #+','+label
3190            else: new_boundary[(point0,point1)]=label
3191
3192        boundary = new_boundary
3193
3194    return coordinates,volumes,boundary
3195
3196
3197def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
3198                 verbose=False):
3199    """Read Digitial Elevation model from the following NetCDF format (.dem)
3200
3201    Example:
3202
3203    ncols         3121
3204    nrows         1800
3205    xllcorner     722000
3206    yllcorner     5893000
3207    cellsize      25
3208    NODATA_value  -9999
3209    138.3698 137.4194 136.5062 135.5558 ..........
3210
3211    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
3212    """
3213
3214    import os
3215    from Scientific.IO.NetCDF import NetCDFFile
3216    from Numeric import Float, zeros, sum, reshape, equal
3217
3218    root = basename_in
3219    inname = root + '.dem'
3220
3221    #Open existing netcdf file to read
3222    infile = NetCDFFile(inname, 'r')
3223    if verbose: print 'Reading DEM from %s' %inname
3224
3225    #Read metadata
3226    ncols = infile.ncols[0]
3227    nrows = infile.nrows[0]
3228    xllcorner = infile.xllcorner[0]
3229    yllcorner = infile.yllcorner[0]
3230    cellsize = infile.cellsize[0]
3231    NODATA_value = infile.NODATA_value[0]
3232    zone = infile.zone[0]
3233    false_easting = infile.false_easting[0]
3234    false_northing = infile.false_northing[0]
3235    projection = infile.projection
3236    datum = infile.datum
3237    units = infile.units
3238
3239    dem_elevation = infile.variables['elevation']
3240
3241    #Get output file name
3242    if basename_out == None:
3243        outname = root + '_' + repr(cellsize_new) + '.dem'
3244    else:
3245        outname = basename_out + '.dem'
3246
3247    if verbose: print 'Write decimated NetCDF file to %s' %outname
3248
3249    #Determine some dimensions for decimated grid
3250    (nrows_stencil, ncols_stencil) = stencil.shape
3251    x_offset = ncols_stencil / 2
3252    y_offset = nrows_stencil / 2
3253    cellsize_ratio = int(cellsize_new / cellsize)
3254    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
3255    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
3256
3257    #Open netcdf file for output
3258    outfile = NetCDFFile(outname, 'w')
3259
3260    #Create new file
3261    outfile.institution = 'Geoscience Australia'
3262    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
3263                           'of spatial point data'
3264    #Georeferencing
3265    outfile.zone = zone
3266    outfile.projection = projection
3267    outfile.datum = datum
3268    outfile.units = units
3269
3270    outfile.cellsize = cellsize_new
3271    outfile.NODATA_value = NODATA_value
3272    outfile.false_easting = false_easting
3273    outfile.false_northing = false_northing
3274
3275    outfile.xllcorner = xllcorner + (x_offset * cellsize)
3276    outfile.yllcorner = yllcorner + (y_offset * cellsize)
3277    outfile.ncols = ncols_new
3278    outfile.nrows = nrows_new
3279
3280    # dimension definition
3281    outfile.createDimension('number_of_points', nrows_new*ncols_new)
3282
3283    # variable definition
3284    outfile.createVariable('elevation', Float, ('number_of_points',))
3285
3286    # Get handle to the variable
3287    elevation = outfile.variables['elevation']
3288
3289    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
3290
3291    #Store data
3292    global_index = 0
3293    for i in range(nrows_new):
3294        if verbose: print 'Processing row %d of %d' %(i, nrows_new)
3295        lower_index = global_index
3296        telev =  zeros(ncols_new, Float)
3297        local_index = 0
3298        trow = i * cellsize_ratio
3299
3300        for j in range(ncols_new):
3301            tcol = j * cellsize_ratio
3302            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
3303
3304            #if dem contains 1 or more NODATA_values set value in
3305            #decimated dem to NODATA_value, else compute decimated
3306            #value using stencil
3307            if sum(sum(equal(tmp, NODATA_value))) > 0:
3308                telev[local_index] = NODATA_value
3309            else:
3310                telev[local_index] = sum(sum(tmp * stencil))
3311
3312            global_index += 1
3313            local_index += 1
3314
3315        upper_index = global_index
3316
3317        elevation[lower_index:upper_index] = telev
3318
3319    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
3320
3321    infile.close()
3322    outfile.close()
3323
3324
3325
3326
3327def tsh2sww(filename, verbose=False): #test_tsh2sww
3328    """
3329    to check if a tsh/msh file 'looks' good.
3330    """
3331
3332
3333    if verbose == True:print 'Creating domain from', filename
3334    domain = pmesh_to_domain_instance(filename, Domain)
3335    if verbose == True:print "Number of triangles = ", len(domain)
3336
3337    domain.smooth = True
3338    domain.format = 'sww'   #Native netcdf visualisation format
3339    file_path, filename = path.split(filename)
3340    filename, ext = path.splitext(filename)
3341    domain.set_name(filename)   
3342    domain.reduction = mean
3343    if verbose == True:print "file_path",file_path
3344    if file_path == "":file_path = "."
3345    domain.set_datadir(file_path)
3346
3347    if verbose == True:
3348        print "Output written to " + domain.get_datadir() + sep + \
3349              domain.get_name() + "." + domain.format
3350    sww = get_dataobject(domain)
3351    sww.store_connectivity()
3352    sww.store_timestep('stage')
3353
3354
3355def asc_csiro2sww(bath_dir,
3356                  elevation_dir,
3357                  ucur_dir,
3358                  vcur_dir,
3359                  sww_file,
3360                  minlat = None, maxlat = None,
3361                  minlon = None, maxlon = None,
3362                  zscale=1,
3363                  mean_stage = 0,
3364                  fail_on_NaN = True,
3365                  elevation_NaN_filler = 0,
3366                  bath_prefix='ba',
3367                  elevation_prefix='el',
3368                  verbose=False):
3369    """
3370    Produce an sww boundary file, from esri ascii data from CSIRO.
3371
3372    Also convert latitude and longitude to UTM. All coordinates are
3373    assumed to be given in the GDA94 datum.
3374
3375    assume:
3376    All files are in esri ascii format
3377
3378    4 types of information
3379    bathymetry
3380    elevation
3381    u velocity
3382    v velocity
3383
3384    Assumptions
3385    The metadata of all the files is the same
3386    Each type is in a seperate directory
3387    One bath file with extention .000
3388    The time period is less than 24hrs and uniform.
3389    """
3390    from Scientific.IO.NetCDF import NetCDFFile
3391
3392    from anuga.coordinate_transforms.redfearn import redfearn
3393
3394    precision = Float # So if we want to change the precision its done here
3395
3396    # go in to the bath dir and load the only file,
3397    bath_files = os.listdir(bath_dir)
3398
3399    bath_file = bath_files[0]
3400    bath_dir_file =  bath_dir + os.sep + bath_file
3401    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
3402
3403    #Use the date.time of the bath file as a basis for
3404    #the start time for other files
3405    base_start = bath_file[-12:]
3406
3407    #go into the elevation dir and load the 000 file
3408    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
3409                         + base_start
3410
3411    elevation_files = os.listdir(elevation_dir)
3412    ucur_files = os.listdir(ucur_dir)
3413    vcur_files = os.listdir(vcur_dir)
3414    elevation_files.sort()
3415    # the first elevation file should be the
3416    # file with the same base name as the bath data
3417    assert elevation_files[0] == 'el' + base_start
3418
3419    number_of_latitudes = bath_grid.shape[0]
3420    number_of_longitudes = bath_grid.shape[1]
3421    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3422
3423    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
3424                  for x in range(number_of_longitudes)]
3425    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
3426                 for y in range(number_of_latitudes)]
3427
3428     # reverse order of lat, so the fist lat represents the first grid row
3429    latitudes.reverse()
3430
3431    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:],
3432                                                 minlat=minlat, maxlat=maxlat,
3433                                                 minlon=minlon, maxlon=maxlon)
3434
3435
3436    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
3437    latitudes = latitudes[kmin:kmax]
3438    longitudes = longitudes[lmin:lmax]
3439    number_of_latitudes = len(latitudes)
3440    number_of_longitudes = len(longitudes)
3441    number_of_times = len(os.listdir(elevation_dir))
3442    number_of_points = number_of_latitudes*number_of_longitudes
3443    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
3444
3445    #Work out the times
3446    if len(elevation_files) > 1:
3447        # Assume: The time period is less than 24hrs.
3448        time_period = (int(elevation_files[1][-3:]) - \
3449                      int(elevation_files[0][-3:]))*60*60
3450        times = [x*time_period for x in range(len(elevation_files))]
3451    else:
3452        times = [0.0]
3453
3454
3455    if verbose:
3456        print '------------------------------------------------'
3457        print 'Statistics:'
3458        print '  Extent (lat/lon):'
3459        print '    lat in [%f, %f], len(lat) == %d'\
3460              %(min(latitudes), max(latitudes),
3461                len(latitudes))
3462        print '    lon in [%f, %f], len(lon) == %d'\
3463              %(min(longitudes), max(longitudes),
3464                len(longitudes))
3465        print '    t in [%f, %f], len(t) == %d'\
3466              %(min(times), max(times), len(times))
3467
3468    ######### WRITE THE SWW FILE #############
3469    # NetCDF file definition
3470    outfile = NetCDFFile(sww_file, 'w')
3471
3472    #Create new file
3473    outfile.institution = 'Geoscience Australia'
3474    outfile.description = 'Converted from XXX'
3475
3476
3477    #For sww compatibility
3478    outfile.smoothing = 'Yes'
3479    outfile.order = 1
3480
3481    #Start time in seconds since the epoch (midnight 1/1/1970)
3482    outfile.starttime = starttime = times[0]
3483
3484
3485    # dimension definitions
3486    outfile.createDimension('number_of_volumes', number_of_volumes)
3487
3488    outfile.createDimension('number_of_vertices', 3)
3489    outfile.createDimension('number_of_points', number_of_points)
3490    outfile.createDimension('number_of_timesteps', number_of_times)
3491
3492    # variable definitions
3493    outfile.createVariable('x', precision, ('number_of_points',))
3494    outfile.createVariable('y', precision, ('number_of_points',))
3495    outfile.createVariable('elevation', precision, ('number_of_points',))
3496
3497    #FIXME: Backwards compatibility
3498    outfile.createVariable('z', precision, ('number_of_points',))
3499    #################################
3500
3501    outfile.createVariable('volumes', Int, ('number_of_volumes',
3502                                            'number_of_vertices'))
3503
3504    outfile.createVariable('time', precision,
3505                           ('number_of_timesteps',))
3506
3507    outfile.createVariable('stage', precision,
3508                           ('number_of_timesteps',
3509                            'number_of_points'))
3510
3511    outfile.createVariable('xmomentum', precision,
3512                           ('number_of_timesteps',
3513                            'number_of_points'))
3514
3515    outfile.createVariable('ymomentum', precision,
3516                           ('number_of_timesteps',
3517                            'number_of_points'))
3518
3519    #Store
3520    from anuga.coordinate_transforms.redfearn import redfearn
3521    x = zeros(number_of_points, Float)  #Easting
3522    y = zeros(number_of_points, Float)  #Northing
3523
3524    if verbose: print 'Making triangular grid'
3525    #Get zone of 1st point.
3526    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
3527
3528    vertices = {}
3529    i = 0
3530    for k, lat in enumerate(latitudes):
3531        for l, lon in enumerate(longitudes):
3532
3533            vertices[l,k] = i
3534
3535            zone, easting, northing = redfearn(lat,lon)
3536
3537            msg = 'Zone boundary crossed at longitude =', lon
3538            #assert zone == refzone, msg
3539            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
3540            x[i] = easting
3541            y[i] = northing
3542            i += 1
3543
3544
3545    #Construct 2 triangles per 'rectangular' element
3546    volumes = []
3547    for l in range(number_of_longitudes-1):    #X direction
3548        for k in range(number_of_latitudes-1): #Y direction
3549            v1 = vertices[l,k+1]
3550            v2 = vertices[l,k]
3551            v3 = vertices[l+1,k+1]
3552            v4 = vertices[l+1,k]
3553
3554            #Note, this is different to the ferrit2sww code
3555            #since the order of the lats is reversed.
3556            volumes.append([v1,v3,v2]) #Upper element
3557            volumes.append([v4,v2,v3]) #Lower element
3558
3559    volumes = array(volumes)
3560
3561    geo_ref = Geo_reference(refzone,min(x),min(y))
3562    geo_ref.write_NetCDF(outfile)
3563
3564    # This will put the geo ref in the middle
3565    #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
3566
3567
3568    if verbose:
3569        print '------------------------------------------------'
3570        print 'More Statistics:'
3571        print '  Extent (/lon):'
3572        print '    x in [%f, %f], len(lat) == %d'\
3573              %(min(x), max(x),
3574                len(x))
3575        print '    y in [%f, %f], len(lon) == %d'\
3576              %(min(y), max(y),
3577                len(y))
3578        print 'geo_ref: ',geo_ref
3579
3580    z = resize(bath_grid,outfile.variables['z'][:].shape)
3581    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
3582    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
3583    outfile.variables['z'][:] = z
3584    outfile.variables['elevation'][:] = z  #FIXME HACK
3585    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
3586
3587    stage = outfile.variables['stage']
3588    xmomentum = outfile.variables['xmomentum']
3589    ymomentum = outfile.variables['ymomentum']
3590
3591    outfile.variables['time'][:] = times   #Store time relative
3592
3593    if verbose: print 'Converting quantities'
3594    n = number_of_times
3595    for j in range(number_of_times):
3596        # load in files
3597        elevation_meta, elevation_grid = \
3598           _read_asc(elevation_dir + os.sep + elevation_files[j])
3599
3600        _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
3601        _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
3602
3603        #cut matrix to desired size
3604        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
3605        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
3606        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
3607       
3608        # handle missing values
3609        missing = (elevation_grid == elevation_meta['NODATA_value'])
3610        if sometrue (missing):
3611            if fail_on_NaN:
3612                msg = 'File %s contains missing values'\
3613                      %(elevation_files[j])
3614                raise DataMissingValuesError, msg
3615            else:
3616                elevation_grid = elevation_grid*(missing==0) + \
3617                                 missing*elevation_NaN_filler
3618
3619
3620        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
3621        i = 0
3622        for k in range(number_of_latitudes):      #Y direction
3623            for l in range(number_of_longitudes): #X direction
3624                w = zscale*elevation_grid[k,l] + mean_stage
3625                stage[j,i] = w
3626                h = w - z[i]
3627                xmomentum[j,i] = u_momentum_grid[k,l]*h
3628                ymomentum[j,i] = v_momentum_grid[k,l]*h
3629                i += 1
3630    outfile.close()
3631
3632def _get_min_max_indexes(latitudes_ref,longitudes_ref,
3633                        minlat=None, maxlat=None,
3634                        minlon=None, maxlon=None):
3635    """
3636    return max, min indexes (for slicing) of the lat and long arrays to cover the area
3637    specified with min/max lat/long
3638
3639    Think of the latitudes and longitudes describing a 2d surface.
3640    The area returned is, if possible, just big enough to cover the
3641    inputed max/min area. (This will not be possible if the max/min area
3642    has a section outside of the latitudes/longitudes area.)
3643
3644    asset  longitudes are sorted,
3645    long - from low to high (west to east, eg 148 - 151)
3646    assert latitudes are sorted, ascending or decending
3647    """
3648    latitudes = latitudes_ref[:]
3649    longitudes = longitudes_ref[:]
3650
3651    latitudes = ensure_numeric(latitudes)
3652    longitudes = ensure_numeric(longitudes)
3653   
3654    assert allclose(sort(longitudes), longitudes)
3655   
3656    lat_ascending = True
3657    if not allclose(sort(latitudes), latitudes):
3658        lat_ascending = False
3659        # reverse order of lat, so it's in ascending order         
3660        latitudes = latitudes[::-1]
3661        assert allclose(sort(latitudes), latitudes)
3662    #print "latitudes  in funct", latitudes
3663   
3664    largest_lat_index = len(latitudes)-1
3665    #Cut out a smaller extent.
3666    if minlat == None:
3667        lat_min_index = 0
3668    else:
3669        lat_min_index = searchsorted(latitudes, minlat)-1
3670        if lat_min_index <0:
3671            lat_min_index = 0
3672
3673
3674    if maxlat == None:
3675        lat_max_index = largest_lat_index #len(latitudes)
3676    else:
3677        lat_max_index = searchsorted(latitudes, maxlat)
3678        if lat_max_index > largest_lat_index:
3679            lat_max_index = largest_lat_index
3680
3681    if minlon == None:
3682        lon_min_index = 0
3683    else:
3684        lon_min_index = searchsorted(longitudes, minlon)-1
3685        if lon_min_index <0:
3686            lon_min_index = 0
3687
3688    if maxlon == None:
3689        lon_max_index = len(longitudes)
3690    else:
3691        lon_max_index = searchsorted(longitudes, maxlon)
3692
3693    # Reversing the indexes, if the lat array is decending
3694    if lat_ascending is False:
3695        lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
3696                                       largest_lat_index - lat_min_index
3697    lat_max_index = lat_max_index + 1 # taking into account how slicing works
3698    lon_max_index = lon_max_index + 1 # taking into account how slicing works
3699
3700    return lat_min_index, lat_max_index, lon_min_index, lon_max_index
3701
3702
3703def _read_asc(filename, verbose=False):
3704    """Read esri file from the following ASCII format (.asc)
3705
3706    Example:
3707
3708    ncols         3121
3709    nrows         1800
3710    xllcorner     722000
3711    yllcorner     5893000
3712    cellsize      25
3713    NODATA_value  -9999
3714    138.3698 137.4194 136.5062 135.5558 ..........
3715
3716    """
3717
3718    datafile = open(filename)
3719
3720    if verbose: print 'Reading DEM from %s' %(filename)
3721    lines = datafile.readlines()
3722    datafile.close()
3723
3724    if verbose: print 'Got', len(lines), ' lines'
3725
3726    ncols = int(lines.pop(0).split()[1].strip())
3727    nrows = int(lines.pop(0).split()[1].strip())
3728    xllcorner = float(lines.pop(0).split()[1].strip())
3729    yllcorner = float(lines.pop(0).split()[1].strip())
3730    cellsize = float(lines.pop(0).split()[1].strip())
3731    NODATA_value = float(lines.pop(0).split()[1].strip())
3732
3733    assert len(lines) == nrows
3734
3735    #Store data
3736    grid = []
3737
3738    n = len(lines)
3739    for i, line in enumerate(lines):
3740        cells = line.split()
3741        assert len(cells) == ncols
3742        grid.append(array([float(x) for x in cells]))
3743    grid = array(grid)
3744
3745    return {'xllcorner':xllcorner,
3746            'yllcorner':yllcorner,
3747            'cellsize':cellsize,
3748            'NODATA_value':NODATA_value}, grid
3749
3750
3751
3752    ####  URS 2 SWW  ###
3753
3754lon_name = 'LON'
3755lat_name = 'LAT'
3756time_name = 'TIME'
3757precision = Float # So if we want to change the precision its done here       
3758class Write_nc:
3759    """
3760    Write an nc file.
3761   
3762    Note, this should be checked to meet cdc netcdf conventions for gridded
3763    data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
3764   
3765    """
3766    def __init__(self,
3767                 quantity_name,
3768                 file_name,
3769                 time_step_count,
3770                 time_step,
3771                 lon,
3772                 lat):
3773        """
3774        time_step_count is the number of time steps.
3775        time_step is the time step size
3776       
3777        pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
3778        """
3779        self.quantity_name = quantity_name
3780        quantity_units = {'HA':'CENTIMETERS',
3781                              'UA':'CENTIMETERS/SECOND',
3782                              'VA':'CENTIMETERS/SECOND'}       
3783       
3784        multiplier_dic = {'HA':100.0, # To convert from m to cm
3785                              'UA':100.0,  #  m/s to cm/sec
3786                              'VA':-100.0}  # MUX files have positve x in the
3787        # Southern direction.  This corrects for it, when writing nc files.
3788       
3789        self.quantity_multiplier =  multiplier_dic[self.quantity_name]
3790       
3791        #self.file_name = file_name
3792        self.time_step_count = time_step_count
3793        self.time_step = time_step
3794
3795        # NetCDF file definition
3796        self.outfile = NetCDFFile(file_name, 'w')
3797        outfile = self.outfile       
3798
3799        #Create new file
3800        nc_lon_lat_header(outfile, lon, lat)
3801   
3802        # TIME
3803        outfile.createDimension(time_name, None)
3804        outfile.createVariable(time_name, precision, (time_name,))
3805
3806        #QUANTITY
3807        outfile.createVariable(self.quantity_name, precision,
3808                               (time_name, lat_name, lon_name))
3809        outfile.variables[self.quantity_name].missing_value=-1.e+034
3810        outfile.variables[self.quantity_name].units= \
3811                                 quantity_units[self.quantity_name]
3812        outfile.variables[lon_name][:]= ensure_numeric(lon)
3813        outfile.variables[lat_name][:]= ensure_numeric(lat)
3814
3815        #Assume no one will be wanting to read this, while we are writing
3816        #outfile.close()
3817       
3818    def store_timestep(self, quantity_slice):
3819        """
3820        Write a time slice of quantity info
3821        quantity_slice is the data to be stored at this time step
3822        """
3823       
3824        outfile = self.outfile
3825       
3826        # Get the variables
3827        time = outfile.variables[time_name]
3828        quantity = outfile.variables[self.quantity_name]
3829           
3830        i = len(time)
3831
3832        #Store time
3833        time[i] = i*self.time_step #self.domain.time
3834        quantity[i,:] = quantity_slice* self.quantity_multiplier
3835       
3836    def close(self):
3837        self.outfile.close()
3838
3839def urs2sww(basename_in='o', basename_out=None, verbose=False,
3840            remove_nc_files=True,
3841            minlat=None, maxlat=None,
3842            minlon= None, maxlon=None,
3843            mint=None, maxt=None,
3844            mean_stage=0,
3845            origin = None,
3846            zscale=1,
3847            fail_on_NaN=True,
3848            NaN_filler=0,
3849            elevation=None):
3850    """
3851    Convert URS C binary format for wave propagation to
3852    sww format native to abstract_2d_finite_volumes.
3853
3854    Specify only basename_in and read files of the form
3855    basefilename_velocity-z-mux, basefilename_velocity-e-mux and
3856    basefilename_waveheight-n-mux containing relative height,
3857    x-velocity and y-velocity, respectively.
3858
3859    Also convert latitude and longitude to UTM. All coordinates are
3860    assumed to be given in the GDA94 datum. The latitude and longitude
3861    information is for  a grid.
3862
3863    min's and max's: If omitted - full extend is used.
3864    To include a value min may equal it, while max must exceed it.
3865    Lat and lon are assumed to be in decimal degrees.
3866    NOTE: minlon is the most east boundary.
3867   
3868    origin is a 3-tuple with geo referenced
3869    UTM coordinates (zone, easting, northing)
3870    It will be the origin of the sww file. This shouldn't be used,
3871    since all of anuga should be able to handle an arbitary origin.
3872
3873
3874    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
3875    which means that latitude is the fastest
3876    varying dimension (row major order, so to speak)
3877
3878    In URS C binary the latitudes and longitudes are in assending order.
3879    """
3880    if basename_out == None:
3881        basename_out = basename_in
3882    files_out = urs2nc(basename_in, basename_out)
3883    ferret2sww(basename_out,
3884               minlat=minlat,
3885               maxlat=maxlat,
3886               minlon=minlon,
3887               maxlon=maxlon,
3888               mint=mint,
3889               maxt=maxt,
3890               mean_stage=mean_stage,
3891               origin=origin,
3892               zscale=zscale,
3893               fail_on_NaN=fail_on_NaN,
3894               NaN_filler=NaN_filler,
3895               inverted_bathymetry=True,
3896               verbose=verbose)
3897    #print "files_out",files_out
3898    if remove_nc_files:
3899        for file_out in files_out:
3900            os.remove(file_out)
3901   
3902def urs2nc(basename_in = 'o', basename_out = 'urs'):
3903    """
3904    Convert the 3 urs files to 4 nc files.
3905
3906    The name of the urs file names must be;
3907    [basename_in]_velocity-z-mux
3908    [basename_in]_velocity-e-mux
3909    [basename_in]_waveheight-n-mux
3910   
3911    """
3912   
3913    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
3914                basename_in + EAST_VELOCITY_LABEL,
3915                basename_in + NORTH_VELOCITY_LABEL]
3916    files_out = [basename_out+'_ha.nc',
3917                 basename_out+'_ua.nc',
3918                 basename_out+'_va.nc']
3919    quantities = ['HA','UA','VA']
3920
3921    for file_name in files_in:
3922        if os.access(file_name, os.F_OK) == 0 :
3923            msg = 'File %s does not exist or is not accessible' %file_name
3924            raise IOError, msg
3925       
3926    hashed_elevation = None
3927    for file_in, file_out, quantity in map(None, files_in,
3928                                           files_out,
3929                                           quantities):
3930        lonlatdep, lon, lat, depth = _binary_c2nc(file_in,
3931                                         file_out,
3932                                         quantity)
3933        #print "lonlatdep", lonlatdep
3934        if hashed_elevation == None:
3935            elevation_file = basename_out+'_e.nc'
3936            write_elevation_nc(elevation_file,
3937                                lon,
3938                                lat,
3939                                depth)
3940            hashed_elevation = myhash(lonlatdep)
3941        else:
3942            msg = "The elevation information in the mux files is inconsistent"
3943            assert hashed_elevation == myhash(lonlatdep), msg
3944    files_out.append(elevation_file)
3945    return files_out
3946   
3947def _binary_c2nc(file_in, file_out, quantity):
3948    """
3949    Reads in a quantity urs file and writes a quantity nc file.
3950    additionally, returns the depth and lat, long info,
3951    so it can be written to a file.
3952    """
3953    columns = 3 # long, lat , depth
3954    mux_file = open(file_in, 'rb')
3955   
3956    # Number of points/stations
3957    (points_num,)= unpack('i',mux_file.read(4))
3958
3959    # nt, int - Number of time steps
3960    (time_step_count,)= unpack('i',mux_file.read(4))
3961
3962    #dt, float - time step, seconds
3963    (time_step,) = unpack('f', mux_file.read(4))
3964   
3965    msg = "Bad data in the mux file."
3966    if points_num < 0:
3967        mux_file.close()
3968        raise ANUGAError, msg
3969    if time_step_count < 0:
3970        mux_file.close()
3971        raise ANUGAError, msg
3972    if time_step < 0:
3973        mux_file.close()
3974        raise ANUGAError, msg
3975   
3976    lonlatdep = p_array.array('f')
3977    lonlatdep.read(mux_file, columns * points_num)
3978    lonlatdep = array(lonlatdep, typecode=Float)   
3979    lonlatdep = reshape(lonlatdep, (points_num, columns))
3980   
3981    lon, lat, depth = lon_lat2grid(lonlatdep)
3982    lon_sorted = list(lon)
3983    lon_sorted.sort()
3984
3985    if not lon == lon_sorted:
3986        msg = "Longitudes in mux file are not in ascending order"
3987        raise IOError, msg
3988    lat_sorted = list(lat)
3989    lat_sorted.sort()
3990
3991    if not lat == lat_sorted:
3992        msg = "Latitudes in mux file are not in ascending order"
3993   
3994    nc_file = Write_nc(quantity,
3995                       file_out,
3996                       time_step_count,
3997                       time_step,
3998                       lon,
3999                       lat)
4000
4001    for i in range(time_step_count):
4002        #Read in a time slice  from mux file 
4003        hz_p_array = p_array.array('f')
4004        hz_p_array.read(mux_file, points_num)
4005        hz_p = array(hz_p_array, typecode=Float)
4006        hz_p = reshape(hz_p, (len(lon), len(lat)))
4007        hz_p = transpose(hz_p) #mux has lat varying fastest, nc has long v.f.
4008
4009        #write time slice to nc file
4010        nc_file.store_timestep(hz_p)
4011    mux_file.close()
4012    nc_file.close()
4013
4014    return lonlatdep, lon, lat, depth
4015   
4016
4017def write_elevation_nc(file_out, lon, lat, depth_vector):
4018    """
4019    Write an nc elevation file.
4020    """
4021   
4022    # NetCDF file definition
4023    outfile = NetCDFFile(file_out, 'w')
4024
4025    #Create new file
4026    nc_lon_lat_header(outfile, lon, lat)
4027   
4028    # ELEVATION
4029    zname = 'ELEVATION'
4030    outfile.createVariable(zname, precision, (lat_name, lon_name))
4031    outfile.variables[zname].units='CENTIMETERS'
4032    outfile.variables[zname].missing_value=-1.e+034
4033
4034    outfile.variables[lon_name][:]= ensure_numeric(lon)
4035    outfile.variables[lat_name][:]= ensure_numeric(lat)
4036
4037    depth = reshape(depth_vector, ( len(lat), len(lon)))
4038    outfile.variables[zname][:]= depth
4039   
4040    outfile.close()
4041   
4042def nc_lon_lat_header(outfile, lon, lat):
4043    """
4044    outfile is the netcdf file handle.
4045    lon - a list/array of the longitudes
4046    lat - a list/array of the latitudes
4047    """
4048   
4049    outfile.institution = 'Geoscience Australia'
4050    outfile.description = 'Converted from URS binary C'
4051   
4052    # Longitude
4053    outfile.createDimension(lon_name, len(lon))
4054    outfile.createVariable(lon_name, precision, (lon_name,))
4055    outfile.variables[lon_name].point_spacing='uneven'
4056    outfile.variables[lon_name].units='degrees_east'
4057    outfile.variables[lon_name].assignValue(lon)
4058
4059
4060    # Latitude
4061    outfile.createDimension(lat_name, len(lat))
4062    outfile.createVariable(lat_name, precision, (lat_name,))
4063    outfile.variables[lat_name].point_spacing='uneven'
4064    outfile.variables[lat_name].units='degrees_north'
4065    outfile.variables[lat_name].assignValue(lat)
4066
4067
4068   
4069def lon_lat2grid(long_lat_dep):
4070    """
4071    given a list of points that are assumed to be an a grid,
4072    return the long's and lat's of the grid.
4073    long_lat_dep is an array where each row is a position.
4074    The first column is longitudes.
4075    The second column is latitudes.
4076
4077    The latitude is the fastest varying dimension - in mux files
4078    """
4079    LONG = 0
4080    LAT = 1
4081    QUANTITY = 2
4082
4083    long_lat_dep = ensure_numeric(long_lat_dep, Float)
4084   
4085    num_points = long_lat_dep.shape[0]
4086    this_rows_long = long_lat_dep[0,LONG]
4087
4088    # Count the length of unique latitudes
4089    i = 0
4090    while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
4091        i += 1
4092    # determine the lats and longsfrom the grid
4093    lat = long_lat_dep[:i, LAT]       
4094    long = long_lat_dep[::i, LONG]
4095   
4096    lenlong = len(long)
4097    lenlat = len(lat)
4098    #print 'len lat', lat, len(lat)
4099    #print 'len long', long, len(long)
4100         
4101    msg = 'Input data is not gridded'     
4102    assert num_points % lenlat == 0, msg
4103    assert num_points % lenlong == 0, msg
4104         
4105    # Test that data is gridded       
4106    for i in range(lenlong):
4107        msg = 'Data is not gridded.  It must be for this operation'
4108        first = i*lenlat
4109        last = first + lenlat
4110               
4111        assert allclose(long_lat_dep[first:last,LAT], lat), msg
4112        assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
4113   
4114   
4115#    print 'range long', min(long), max(long)
4116#    print 'range lat', min(lat), max(lat)
4117#    print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0])
4118#    print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1])
4119   
4120   
4121   
4122    msg = 'Out of range latitudes/longitudes'
4123    for l in lat:assert -90 < l < 90 , msg
4124    for l in long:assert -180 < l < 180 , msg
4125
4126    # Changing quantity from lat being the fastest varying dimension to
4127    # long being the fastest varying dimension
4128    # FIXME - make this faster/do this a better way
4129    # use numeric transpose, after reshaping the quantity vector
4130#    quantity = zeros(len(long_lat_dep), Float)
4131    quantity = zeros(num_points, Float)
4132   
4133#    print 'num',num_points
4134    for lat_i, _ in enumerate(lat):
4135        for long_i, _ in enumerate(long):
4136            q_index = lat_i*lenlong+long_i
4137            lld_index = long_i*lenlat+lat_i
4138#            print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index
4139            temp = long_lat_dep[lld_index, QUANTITY]
4140            quantity[q_index] = temp
4141           
4142    return long, lat, quantity
4143
4144    ####  END URS 2 SWW  ###
4145
4146    #### URS UNGRIDDED 2 SWW ###
4147
4148    ### PRODUCING THE POINTS NEEDED FILE ###
4149LL_LAT = -50.0
4150LL_LONG = 80.0
4151GRID_SPACING = 1.0/60.0
4152LAT_AMOUNT = 4800
4153LONG_AMOUNT = 3600
4154def URS_points_needed_to_file(file_name, boundary_polygon, zone,
4155                              ll_lat=LL_LAT, ll_long=LL_LONG,
4156                              grid_spacing=GRID_SPACING, 
4157                              lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4158                              export_csv=False, use_cache=False,
4159                              verbose=False):
4160    """
4161    file_name - name of the urs file produced for David.
4162    boundary_polygon - a list of points that describes a polygon.
4163                      The last point is assumed ot join the first point.
4164                      This is in UTM (lat long would be better though)
4165
4166    ll_lat - lower left latitude, in decimal degrees
4167    ll-long - lower left longitude, in decimal degrees
4168    grid_spacing - in deciamal degrees
4169
4170
4171    Don't add the file extension.  It will be added.
4172    """
4173    geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long,
4174                            grid_spacing, 
4175                      lat_amount, long_amount,use_cache, verbose)
4176    if not file_name[-4:] == ".urs":
4177        file_name += ".urs"
4178    geo.export_points_file(file_name)
4179    if export_csv:
4180        if file_name[-4:] == ".urs":
4181            file_name = file_name[:-4] + ".csv"
4182        geo.export_points_file(file_name)
4183
4184def URS_points_needed(boundary_polygon, zone, ll_lat=LL_LAT,
4185                      ll_long=LL_LONG, grid_spacing=GRID_SPACING, 
4186                      lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT,
4187                      use_cache=False, verbose=False):
4188    args = (boundary_polygon,
4189                      zone)
4190    kwargs = {'ll_lat': ll_lat,
4191              'll_long': ll_long,
4192              'grid_spacing': grid_spacing,
4193              'lat_amount': lat_amount,
4194              'long_amount': long_amount} 
4195    if use_cache is True:
4196        try:
4197            from anuga.caching import cache
4198        except:
4199            msg = 'Caching was requested, but caching module'+\
4200                  'could not be imported'
4201            raise msg
4202
4203
4204        geo = cache(_URS_points_needed,
4205                  args, kwargs,
4206                  verbose=verbose,
4207                  compression=False)
4208    else:
4209        #I was getting 'got multiple values for keyword argument' errors
4210        #geo = apply(_URS_points_needed, args, kwargs)
4211        geo = _URS_points_needed(boundary_polygon,
4212                      zone, ll_lat,
4213                      ll_long, grid_spacing, 
4214                      lat_amount, long_amount)
4215
4216    return geo   
4217def _URS_points_needed(boundary_polygon,
4218                      zone, ll_lat=LL_LAT,
4219                      ll_long=LL_LONG, grid_spacing=GRID_SPACING, 
4220                      lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT):
4221    """
4222
4223    boundary_polygon - a list of points that describes a polygon.
4224                      The last point is assumed ot join the first point.
4225                      This is in UTM (lat long would be better though)
4226
4227    ll_lat - lower left latitude, in decimal degrees
4228    ll-long - lower left longitude, in decimal degrees
4229    grid_spacing - in deciamal degrees
4230
4231    """
4232   
4233    from sets import ImmutableSet
4234   
4235    msg = "grid_spacing can not be zero"
4236    assert not grid_spacing ==0, msg
4237    a = boundary_polygon
4238    # List of segments.  Each segment is two points.
4239    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
4240
4241    # convert the segs to Lat's and longs.
4242   
4243    # Don't assume the zone of the segments is the same as the lower left
4244    # corner of the lat long data!!  They can easily be in different zones
4245   
4246    lat_long_set = ImmutableSet()
4247    for seg in segs:
4248        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, 
4249                      lat_amount, long_amount, zone)
4250        lat_long_set |= ImmutableSet(points_lat_long)
4251    #print "lat_long_set",lat_long_set
4252    geo = Geospatial_data(data_points=list(lat_long_set),
4253                              points_are_lats_longs=True)
4254    return geo
4255   
4256def points_needed(seg, ll_lat, ll_long, grid_spacing, 
4257                  lat_amount, long_amount, zone):
4258    """
4259    return a list of the points, in lats and longs that are needed to
4260    interpolate any point on the segment.
4261    """
4262    from math import sqrt
4263    #print "zone",zone
4264    geo_reference = Geo_reference(zone=zone)
4265    #print "seg",seg
4266    geo = Geospatial_data(seg,geo_reference=geo_reference)
4267    seg_lat_long = geo.get_data_points(as_lat_long=True)
4268    #print "seg_lat_long", seg_lat_long
4269    # 1.415 = 2^0.5, rounded up....
4270    sqrt_2_rounded_up = 1.415
4271    buffer = sqrt_2_rounded_up * grid_spacing
4272   
4273    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
4274    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
4275    min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
4276    min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
4277
4278    first_row = (min_long - ll_long)/grid_spacing
4279    # To round up
4280    first_row_long = int(round(first_row + 0.5))
4281    #print "first_row", first_row_long
4282
4283    last_row = (max_long - ll_long)/grid_spacing # round down
4284    last_row_long = int(round(last_row))
4285    #print "last_row",last_row _long
4286   
4287    first_row = (min_lat - ll_lat)/grid_spacing
4288    # To round up
4289    first_row_lat = int(round(first_row + 0.5))
4290    #print "first_row", first_row_lat
4291
4292    last_row = (max_lat - ll_lat)/grid_spacing # round down
4293    last_row_lat = int(round(last_row))
4294    #print "last_row",last_row_lat
4295
4296    # to work out the max distance -
4297    # 111120 - horizontal distance between 1 deg latitude.
4298    #max_distance = sqrt_2_rounded_up * 111120 * grid_spacing
4299    max_distance = 157147.4112 * grid_spacing
4300    #print "max_distance", max_distance #2619.12 m for 1 minute
4301    points_lat_long = []
4302    # Create a list of the lat long points to include.
4303    for index_lat in range(first_row_lat, last_row_lat + 1):
4304        for index_long in range(first_row_long, last_row_long + 1):
4305            lat = ll_lat + index_lat*grid_spacing
4306            long = ll_long + index_long*grid_spacing
4307
4308            #filter here to keep good points
4309            if keep_point(lat, long, seg, max_distance):
4310                points_lat_long.append((lat, long)) #must be hashable
4311    #print "points_lat_long", points_lat_long
4312
4313    # Now that we have these points, lets throw ones out that are too far away
4314    return points_lat_long
4315
4316def keep_point(lat, long, seg, max_distance):
4317    """
4318    seg is two points, UTM
4319    """
4320    from math import sqrt
4321    _ , x0, y0 = redfearn(lat, long)
4322    x1 = seg[0][0]
4323    y1 = seg[0][1]
4324    x2 = seg[1][0]
4325    y2 = seg[1][1]
4326
4327    x2_1 = x2-x1
4328    y2_1 = y2-y1
4329    d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1))
4330    if d <= max_distance:
4331        return True
4332    else:
4333        return False
4334   
4335    #### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE ####
4336   
4337WAVEHEIGHT_MUX_LABEL = '_waveheight-z-mux'
4338EAST_VELOCITY_LABEL =  '_velocity-e-mux'
4339NORTH_VELOCITY_LABEL =  '_velocity-n-mux' 
4340def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
4341                      mint=None, maxt=None,
4342                      mean_stage=0,
4343                      origin=None,
4344                      hole_points_UTM=None,
4345                      zscale=1):
4346    """   
4347    Convert URS C binary format for wave propagation to
4348    sww format native to abstract_2d_finite_volumes.
4349
4350
4351    Specify only basename_in and read files of the form
4352    basefilename_velocity-z-mux, basefilename_velocity-e-mux and
4353    basefilename_waveheight-n-mux containing relative height,
4354    x-velocity and y-velocity, respectively.
4355
4356    Also convert latitude and longitude to UTM. All coordinates are
4357    assumed to be given in the GDA94 datum. The latitude and longitude
4358    information is assumed ungridded grid.
4359
4360    min's and max's: If omitted - full extend is used.
4361    To include a value min ans max may equal it.
4362    Lat and lon are assumed to be in decimal degrees.
4363   
4364    origin is a 3-tuple with geo referenced
4365    UTM coordinates (zone, easting, northing)
4366    It will be the origin of the sww file. This shouldn't be used,
4367    since all of anuga should be able to handle an arbitary origin.
4368
4369
4370    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
4371    which means that latitude is the fastest
4372    varying dimension (row major order, so to speak)
4373
4374    In URS C binary the latitudes and longitudes are in assending order.
4375
4376    Note, interpolations of the resulting sww file will be different
4377    from results of urs2sww.  This is due to the interpolation
4378    function used, and the different grid structure between urs2sww
4379    and this function.
4380   
4381    Interpolating data that has an underlying gridded source can
4382    easily end up with different values, depending on the underlying
4383    mesh.
4384
4385    consider these 4 points
4386    50  -50
4387
4388    0     0
4389
4390    The grid can be
4391     -
4392    |\|    A
4393     -
4394     or;
4395      -
4396     |/|   B
4397      -
4398      If a point is just below the center of the midpoint, it will have a
4399      +ve value in grid A and a -ve value in grid B.
4400    """ 
4401    from anuga.pmesh.mesh import Mesh, NoTrianglesError
4402
4403    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4404                basename_in + EAST_VELOCITY_LABEL,
4405                basename_in + NORTH_VELOCITY_LABEL]
4406    quantities = ['HA','UA','VA']
4407
4408    # instanciate urs_points of the three mux files.
4409    mux = {}
4410    for quantity, file in map(None, quantities, files_in):
4411        mux[quantity] = Urs_points(file)
4412       
4413    # Could check that the depth is the same. (hashing)
4414
4415    # handle to a mux file to do depth stuff
4416    a_mux = mux[quantities[0]]
4417   
4418    # Convert to utm
4419    lat = a_mux.lonlatdep[:,1]
4420    long = a_mux.lonlatdep[:,0]
4421    points_utm, zone = convert_from_latlon_to_utm( \
4422        latitudes=lat, longitudes=long)
4423    #print "points_utm", points_utm
4424    #print "zone", zone
4425
4426    elevation = a_mux.lonlatdep[:,2] * -1 #
4427   
4428    # grid ( create a mesh from the selected points)
4429    # This mesh has a problem.  Triangles are streched over ungridded areas.
4430    #  If these areas could be described as holes in pmesh, that would be great
4431
4432    # I can't just get the user to selection a point in the middle.
4433    # A boundary is needed around these points.
4434    # But if the zone of points is obvious enough auto-segment should do
4435    # a good boundary.
4436    mesh = Mesh()
4437    mesh.add_vertices(points_utm)
4438    mesh.auto_segment()
4439    if hole_points_UTM is not None:
4440        point = ensure_absolute(hole_points_UTM)
4441        mesh.add_hole(point[0], point[1])
4442    try:
4443        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4444    except NoTrianglesError:
4445        # This is a bit of a hack, going in and changing the
4446        # data structure.
4447        mesh.holes = []
4448        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
4449    mesh_dic = mesh.Mesh2MeshList()
4450
4451    #mesh.export_mesh_file(basename_in + '.tsh')
4452    # These are the times of the mux file
4453    mux_times = []
4454    for i in range(a_mux.time_step_count):
4455        mux_times.append(a_mux.time_step * i) 
4456    mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt)
4457    times = mux_times[mux_times_start_i:mux_times_fin_i]
4458   
4459    if mux_times_start_i == mux_times_fin_i:
4460        # Close the mux files
4461        for quantity, file in map(None, quantities, files_in):
4462            mux[quantity].close()
4463        msg="Due to mint and maxt there's no time info in the boundary SWW."
4464        assert 1==0 , msg #Happy to know a better way of doing this..
4465       
4466    # If this raise is removed there is currently no downstream errors
4467           
4468    points_utm=ensure_numeric(points_utm)
4469    #print "mesh_dic['generatedpointlist']", mesh_dic['generatedpointlist']
4470    #print "points_utm", points_utm
4471    assert ensure_numeric(mesh_dic['generatedpointlist']) == \
4472           ensure_numeric(points_utm)
4473   
4474    volumes = mesh_dic['generatedtrianglelist']
4475   
4476    # write sww intro and grid stuff.   
4477    if basename_out is None:
4478        swwname = basename_in + '.sww'
4479    else:
4480        swwname = basename_out + '.sww'
4481
4482    if verbose: print 'Output to ', swwname
4483    outfile = NetCDFFile(swwname, 'w')
4484    # For a different way of doing this, check out tsh2sww
4485    # work out sww_times and the index range this covers
4486    write_sww_header(outfile, times, len(volumes), len(points_utm))
4487    outfile.mean_stage = mean_stage
4488    outfile.zscale = zscale
4489    write_sww_triangulation(outfile, points_utm, volumes,
4490                            elevation, zone,  origin=origin, verbose=verbose)
4491   
4492    if verbose: print 'Converting quantities'
4493    j = 0
4494    # Read in a time slice from each mux file and write it to the sww file
4495    for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
4496        if j >= mux_times_start_i and j < mux_times_fin_i:
4497            write_sww_time_slice(outfile, HA, UA, VA,  elevation,
4498                                 j - mux_times_start_i,
4499                                 mean_stage=mean_stage, zscale=zscale,
4500                                 verbose=verbose)
4501        j += 1
4502    outfile.close()
4503    #
4504    # Do some conversions while writing the sww file
4505def mux2sww_time(mux_times, mint, maxt):
4506    """
4507    """
4508
4509    if mint == None:
4510        mux_times_start_i = 0
4511    else:
4512        mux_times_start_i = searchsorted(mux_times, mint)
4513       
4514    if maxt == None:
4515        mux_times_fin_i = len(mux_times)
4516    else:
4517        maxt += 0.5 # so if you specify a time where there is
4518                    # data that time will be included
4519        mux_times_fin_i = searchsorted(mux_times, maxt)
4520
4521    return mux_times_start_i, mux_times_fin_i
4522
4523def write_sww_header(outfile, times, number_of_volumes, number_of_points, description='Converted from XXX'):
4524    """
4525    outfile - the name of the file that will be written
4526    times - A list of the time slice times
4527    number_of_volumes - the number of triangles
4528    """
4529    times = ensure_numeric(times)
4530   
4531    number_of_times = len(times)
4532   
4533    outfile.institution = 'Geoscience Australia'
4534    outfile.description = 'Converted from XXX'
4535
4536    #For sww compatibility
4537    outfile.smoothing = 'Yes'
4538    outfile.order = 1
4539
4540    #Start time in seconds since the epoch (midnight 1/1/1970)
4541    if len(times) == 0:
4542        outfile.starttime = starttime = 0
4543    else:
4544        outfile.starttime = starttime = times[0]
4545    times = times - starttime  #Store relative times
4546
4547    # dimension definitions
4548    outfile.createDimension('number_of_volumes', number_of_volumes)
4549    outfile.createDimension('number_of_vertices', 3)
4550    outfile.createDimension('number_of_points', number_of_points)
4551    outfile.createDimension('number_of_timesteps', number_of_times)
4552
4553    # variable definitions
4554    outfile.createVariable('x', precision, ('number_of_points',))
4555    outfile.createVariable('y', precision, ('number_of_points',))
4556    outfile.createVariable('elevation', precision, ('number_of_points',))
4557
4558    #FIXME: Backwards compatibility
4559    outfile.createVariable('z', precision, ('number_of_points',))
4560    #################################
4561
4562    outfile.createVariable('volumes', Int, ('number_of_volumes',
4563                                            'number_of_vertices'))
4564
4565    outfile.createVariable('time', precision,
4566                           ('number_of_timesteps',))
4567
4568    outfile.createVariable('stage', precision,
4569                           ('number_of_timesteps',
4570                            'number_of_points'))
4571
4572    outfile.createVariable('xmomentum', precision,
4573                           ('number_of_timesteps',
4574                            'number_of_points'))
4575
4576    outfile.createVariable('ymomentum', precision,
4577                           ('number_of_timesteps',
4578                            'number_of_points'))
4579   
4580    outfile.variables['time'][:] = times    #Store time relative
4581
4582
4583def write_sww_triangulation(outfile, points_utm, volumes,
4584                            elevation, zone, origin=None, verbose=False):
4585
4586    number_of_points = len(points_utm)   
4587    volumes = array(volumes)
4588
4589    if origin is None:
4590        origin = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
4591    geo_ref = write_NetCDF_georeference(origin, outfile)
4592   
4593#     if origin is not None:
4594#         if isinstance(origin, Geo_reference):
4595#             geo_ref = origin
4596#         else:
4597#             geo_ref = apply(Geo_reference,origin)
4598#     else:
4599#         geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
4600#     #geo_ref = Geo_reference(zone,0,0)
4601#     geo_ref.write_NetCDF(outfile)
4602   
4603
4604    # This will put the geo ref in the middle
4605    #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
4606    x =  points_utm[:,0]
4607    y =  points_utm[:,1]
4608    z = outfile.variables['z'][:]
4609    #FIXME the w.r.t time is lost now..
4610    if verbose:
4611        print '------------------------------------------------'
4612        print 'More Statistics:'
4613        print '  Extent (/lon):'
4614        print '    x in [%f, %f], len(lat) == %d'\
4615              %(min(x), max(x),
4616                len(x))
4617        print '    y in [%f, %f], len(lon) == %d'\
4618              %(min(y), max(y),
4619                len(y))
4620        print '    z in [%f, %f], len(z) == %d'\
4621              %(min(elevation), max(elevation),
4622                len(elevation))
4623        print 'geo_ref: ',geo_ref
4624        print '------------------------------------------------'
4625
4626    #z = resize(bath_grid,outfile.variables['z'][:].shape)
4627    outfile.variables['x'][:] = points_utm[:,0] - geo_ref.get_xllcorner()
4628    outfile.variables['y'][:] = points_utm[:,1] - geo_ref.get_yllcorner()
4629    outfile.variables['z'][:] = elevation
4630    outfile.variables['elevation'][:] = elevation  #FIXME HACK
4631    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
4632
4633
4634def write_sww_time_slices(outfile, has, uas, vas, elevation,
4635                         mean_stage=0, zscale=1,
4636                         verbose=False):   
4637    #Time stepping
4638    stage = outfile.variables['stage']
4639    xmomentum = outfile.variables['xmomentum']
4640    ymomentum = outfile.variables['ymomentum']
4641
4642    n = len(has)
4643    j=0
4644    for ha, ua, va in map(None, has, uas, vas):
4645        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
4646        w = zscale*ha + mean_stage
4647        stage[j] = w
4648        h = w - elevation
4649        xmomentum[j] = ua*h
4650        ymomentum[j] = -1*va*#  -1 since in mux files south is positive.
4651        j += 1
4652
4653def write_sww_time_slice(outfile, ha, ua, va, elevation, slice_index,
4654                         mean_stage=0, zscale=1,
4655                         verbose=False):   
4656    #Time stepping
4657    stage = outfile.variables['stage']
4658    xmomentum = outfile.variables['xmomentum']
4659    ymomentum = outfile.variables['ymomentum']
4660   
4661    w = zscale*ha + mean_stage
4662    stage[slice_index] = w
4663    h = w - elevation
4664    xmomentum[slice_index] = ua*h
4665    ymomentum[slice_index] = -1*va*h # -1 since in mux files south is positive.
4666
4667def urs2txt(basename_in, location_index=None):
4668    """
4669    Not finished or tested
4670    """
4671   
4672    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
4673                basename_in + EAST_VELOCITY_LABEL,
4674                basename_in + NORTH_VELOCITY_LABEL]
4675    quantities = ['HA','UA','VA']
4676
4677    d = ","
4678   
4679    # instanciate urs_points of the three mux files.
4680    mux = {}
4681    for quantity, file in map(None, quantities, files_in):
4682        mux[quantity] = Urs_points(file)
4683       
4684    # Could check that the depth is the same. (hashing)
4685
4686    # handle to a mux file to do depth stuff
4687    a_mux = mux[quantities[0]]
4688   
4689    # Convert to utm
4690    latitudes = a_mux.lonlatdep[:,1]
4691    longitudes = a_mux.lonlatdep[:,0]
4692    points_utm, zone = convert_from_latlon_to_utm( \
4693        latitudes=latitudes, longitudes=longitudes)
4694    #print "points_utm", points_utm
4695    #print "zone", zone
4696    depths = a_mux.lonlatdep[:,2]  #
4697   
4698    fid = open(basename_in + '.txt', 'w')
4699
4700    fid.write("zone: " + str(zone) + "\n")
4701
4702    if location_index is not None:
4703        #Title
4704        li = location_index
4705        fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
4706                  'Northing' + "\n")
4707        fid.write(str(li) +d+ str(latitudes[li])+d+ \
4708              str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \
4709              str(points_utm[li][01]) + "\n")
4710
4711    # the non-time dependent stuff
4712    #Title
4713    fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
4714                  'Northing' +d+ 'depth m' + "\n")
4715    i = 0
4716    for depth, point_utm, lat, long in map(None, depths,
4717                                               points_utm, latitudes,
4718                                               longitudes):
4719       
4720        fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \
4721                  str(point_utm[01]) +d+ str(depth) + "\n")
4722        i +=1
4723    #Time dependent
4724    if location_index is not None:
4725        time_step = a_mux.time_step
4726        i = 0
4727        #Title
4728        fid.write('time' +d+ 'HA depth m'+d+ \
4729                      'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \
4730                      + "\n")
4731        for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
4732            fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \
4733                      str(UA[location_index]) +d+ str(VA[location_index]) \
4734                      + "\n")
4735           
4736            i +=1
4737   
4738class Urs_points:
4739    """
4740    Read the info in URS mux files.
4741
4742    for the quantities heres a correlation between the file names and
4743    what they mean;
4744    z-mux is height above sea level, m
4745    e-mux is velocity is Eastern direction, m/s
4746    n-mux is velocity is Northern direction, m/s   
4747    """
4748    def __init__(self,urs_file):
4749        self.iterated = False
4750        columns = 3 # long, lat , depth
4751        mux_file = open(urs_file, 'rb')
4752       
4753        # Number of points/stations
4754        (self.points_num,)= unpack('i',mux_file.read(4))
4755       
4756        # nt, int - Number of time steps
4757        (self.time_step_count,)= unpack('i',mux_file.read(4))
4758        #print "self.time_step_count", self.time_step_count
4759        #dt, float - time step, seconds
4760        (self.time_step,) = unpack('f', mux_file.read(4))
4761        #print "self.time_step", self.time_step
4762        msg = "Bad data in the urs file."
4763        if self.points_num < 0:
4764            mux_file.close()
4765            raise ANUGAError, msg
4766        if self.time_step_count < 0:
4767            mux_file.close()
4768            raise ANUGAError, msg
4769        if self.time_step < 0:
4770            mux_file.close()
4771            raise ANUGAError, msg
4772
4773        # the depth is in meters, and it is the distance from the ocean
4774        # to the sea bottom.
4775        lonlatdep = p_array.array('f')
4776        lonlatdep.read(mux_file, columns * self.points_num)
4777        lonlatdep = array(lonlatdep, typecode=Float)   
4778        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
4779        #print 'lonlatdep',lonlatdep
4780        self.lonlatdep = lonlatdep
4781       
4782        self.mux_file = mux_file
4783        # check this array
4784
4785    def __iter__(self):
4786        """
4787        iterate over quantity data which is with respect to time.
4788
4789        Note: You can only interate once over an object
4790       
4791        returns quantity infomation for each time slice
4792        """
4793        msg =  "You can only interate once over a urs file."
4794        assert not self.iterated, msg
4795        self.iter_time_step = 0
4796        self.iterated = True
4797        return self
4798   
4799    def next(self):
4800        if self.time_step_count == self.iter_time_step:
4801            self.close()
4802            raise StopIteration
4803        #Read in a time slice  from mux file 
4804        hz_p_array = p_array.array('f')
4805        hz_p_array.read(self.mux_file, self.points_num)
4806        hz_p = array(hz_p_array, typecode=Float)
4807        self.iter_time_step += 1
4808       
4809        return hz_p
4810
4811    def close(self):
4812        self.mux_file.close()
4813
4814    #### END URS UNGRIDDED 2 SWW ###
4815
4816#-------------------------------------------------------------
4817if __name__ == "__main__":   
4818        pass
4819
Note: See TracBrowser for help on using the repository browser.