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

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

comments and minor changes

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