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

Last change on this file since 3750 was 3750, checked in by duncan, 18 years ago

more error handling

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