Changeset 4455


Ignore:
Timestamp:
May 16, 2007, 3:13:52 PM (17 years ago)
Author:
duncan
Message:

Adding changes to sww file format. This revision has been tested as stable.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r4418 r4455  
    5555class DataTimeError(exceptions.Exception): pass
    5656class DataDomainError(exceptions.Exception): pass
     57class NewQuantity(exceptions.Exception): pass
    5758
    5859
     
    6667
    6768from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
    68      searchsorted, zeros, allclose, around, reshape, transpose, sort
     69     searchsorted, zeros, allclose, around, reshape, transpose, sort, \
     70     NewAxis, ArrayType
    6971from Scientific.IO.NetCDF import NetCDFFile
    7072
     
    7375     convert_from_latlon_to_utm
    7476from anuga.coordinate_transforms.geo_reference import Geo_reference, \
    75      write_NetCDF_georeference
     77     write_NetCDF_georeference, ensure_geo_reference
    7678from anuga.geospatial_data.geospatial_data import Geospatial_data,\
    7779     ensure_absolute
     
    8385from anuga.abstract_2d_finite_volumes.pmesh2domain import \
    8486     pmesh_to_domain_instance
     87from anuga.abstract_2d_finite_volumes.util import get_revision_number
    8588   
    8689def make_filename(s):
     
    280283
    281284        if mode == 'w':
    282 
    283             #Create new file
    284             fid.institution = 'Geoscience Australia'
    285             fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
    286 
    287             if domain.smooth:
    288                 fid.smoothing = 'Yes'
    289             else:
    290                 fid.smoothing = 'No'
    291 
    292             fid.order = domain.default_order
    293 
     285            description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
     286            self.writer = Write_sww()
     287            self.writer.header(fid, domain.starttime,
     288                                         self.number_of_volumes,
     289                                         self.domain.number_of_full_nodes,
     290                                         description=description,
     291                                         smoothing=domain.smooth,
     292                                         order=domain.default_order)
    294293            if hasattr(domain, 'texture'):
    295                 fid.texture = domain.texture
    296             #else:
    297             #    fid.texture = 'None'
    298 
    299             #Reference point
    300             #Start time in seconds since the epoch (midnight 1/1/1970)
    301             #FIXME: Use Georef
    302             fid.starttime = domain.starttime
    303 
    304             # dimension definitions
    305             fid.createDimension('number_of_volumes', self.number_of_volumes)
    306             fid.createDimension('number_of_vertices', 3)
    307 
    308             if domain.smooth is True:
    309                 #fid.createDimension('number_of_points', len(domain.vertexlist))
    310                 fid.createDimension('number_of_points', self.domain.number_of_full_nodes)
    311 
    312                 # FIXME(Ole): This will cause sww files for paralle domains to
    313                 # have ghost nodes stored (but not used by triangles).
    314                 # To clean this up, we have to change get_vertex_values and friends in
    315                 # quantity.py (but I can't be bothered right now)
    316             else:
    317                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
    318 
    319             fid.createDimension('number_of_timesteps', None) #extensible
    320 
    321             # variable definitions
    322             fid.createVariable('x', self.precision, ('number_of_points',))
    323             fid.createVariable('y', self.precision, ('number_of_points',))
    324             fid.createVariable('elevation', self.precision, ('number_of_points',))
    325             if domain.geo_reference is not None:
    326                 domain.geo_reference.write_NetCDF(fid)
    327 
    328             #FIXME: Backwards compatibility
    329             fid.createVariable('z', self.precision, ('number_of_points',))
    330             #################################
    331 
    332             fid.createVariable('volumes', Int, ('number_of_volumes',
    333                                                 'number_of_vertices'))
    334 
    335             fid.createVariable('time', Float,  # Always use full precision lest two timesteps
    336                                                # close to each other may appear as the same step
    337                                ('number_of_timesteps',))
    338 
    339             fid.createVariable('stage', self.precision,
    340                                ('number_of_timesteps',
    341                                 'number_of_points'))
    342 
    343             fid.createVariable('xmomentum', self.precision,
    344                                ('number_of_timesteps',
    345                                 'number_of_points'))
    346 
    347             fid.createVariable('ymomentum', self.precision,
    348                                ('number_of_timesteps',
    349                                 'number_of_points'))
     294                fid.texture = domain.texture               
     295    #        if domain.geo_reference is not None:
     296    #            domain.geo_reference.write_NetCDF(fid)
     297           
     298#             fid.institution = 'Geoscience Australia'
     299#             fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
     300
     301#             if domain.smooth:
     302#                 fid.smoothing = 'Yes'
     303#             else:
     304#                 fid.smoothing = 'No'
     305
     306#             fid.order = domain.default_order
     307
     308#             #Reference point
     309#             #Start time in seconds since the epoch (midnight 1/1/1970)
     310#             #FIXME: Use Georef
     311#             fid.starttime = domain.starttime
     312
     313#             # dimension definitions
     314#             fid.createDimension('number_of_volumes', self.number_of_volumes)
     315#             fid.createDimension('number_of_vertices', 3)
     316
     317#             if domain.smooth is True:
     318#                 #fid.createDimension('number_of_points', len(domain.vertexlist))
     319#                 fid.createDimension('number_of_points', self.domain.number_of_full_nodes)
     320
     321#                 # FIXME(Ole): This will cause sww files for paralle domains to
     322#                 # have ghost nodes stored (but not used by triangles).
     323#                 # To clean this up, we have to change get_vertex_values and friends in
     324#                 # quantity.py (but I can't be bothered right now)
     325#             else:
     326#                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
     327
     328#             fid.createDimension('number_of_timesteps', None) #extensible
     329
     330#             # variable definitions
     331#             fid.createVariable('x', self.precision, ('number_of_points',))
     332#             fid.createVariable('y', self.precision, ('number_of_points',))
     333#             fid.createVariable('elevation', self.precision, ('number_of_points',))
     334#             if domain.geo_reference is not None:
     335#                 domain.geo_reference.write_NetCDF(fid)
     336
     337#             #FIXME: Backwards compatibility
     338#             fid.createVariable('z', self.precision, ('number_of_points',))
     339#             #################################
     340
     341#             fid.createVariable('volumes', Int, ('number_of_volumes',
     342#                                                 'number_of_vertices'))
     343
     344#             fid.createVariable('time', Float,  # Always use full precision lest two timesteps
     345#                                              # close to each other may appear as the same step
     346#                                ('number_of_timesteps',))
     347
     348#             fid.createVariable('stage', self.precision,
     349#                                ('number_of_timesteps',
     350#                                 'number_of_points'))
     351
     352#             fid.createVariable('xmomentum', self.precision,
     353#                                ('number_of_timesteps',
     354#                                 'number_of_points'))
     355
     356#             fid.createVariable('ymomentum', self.precision,
     357#                                ('number_of_timesteps',
     358#                                 'number_of_points'))
    350359
    351360        #Close
     
    381390                                      precision=self.precision)
    382391
    383         volumes[:] = V.astype(volumes.typecode())
    384         x[:] = X
    385         y[:] = Y
    386         z[:] = Z
    387 
    388         #FIXME: Backwards compatibility
    389         z = fid.variables['z']
    390         z[:] = Z
     392        #
     393        points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
     394        self.writer.triangulation(fid, points,
     395                                            V.astype(volumes.typecode()),
     396                                            Z,
     397                                            points_georeference= \
     398                                            domain.geo_reference)
     399        #
     400                               
     401#         volumes[:] = V.astype(volumes.typecode())
     402#         x[:] = X
     403#         y[:] = Y
     404#         z[:] = Z
     405
     406#         #FIXME: Backwards compatibility
     407#         z = fid.variables['z']
     408#         z[:] = Z
    391409        ################################
    392410
     
    405423
    406424        from Numeric import choose
    407 
    408         #Get NetCDF
     425       
     426        #Get NetCDF       
    409427        retries = 0
    410428        file_open = False
     
    486504            ymomentum = fid.variables['ymomentum']
    487505            i = len(time)
    488 
    489             #Store time
    490             time[i] = self.domain.time
    491 
    492 
    493506            if type(names) not in [types.ListType, types.TupleType]:
    494507                names = [names]
    495508
    496             for name in names:
    497                 # Get quantity
    498                 Q = domain.quantities[name]
    499                 A,V = Q.get_vertex_values(xy = False,
     509            if 'stage' in names and 'xmomentum' in names and \
     510                   'ymomentum' in names:
     511
     512                # Get stage
     513                Q = domain.quantities['stage']
     514                A,_ = Q.get_vertex_values(xy = False,
     515                                          precision = self.precision)               
     516                z = fid.variables['elevation']
     517                stage = choose(A-z[:] >= self.minimum_storable_height,
     518                           (z[:], A))
     519               
     520                # Get xmomentum
     521                Q = domain.quantities['xmomentum']
     522                xmomentum, _ = Q.get_vertex_values(xy = False,
    500523                                          precision = self.precision)
    501 
    502 
    503                 #FIXME: Make this general (see below)
    504                 if name == 'stage':
    505                     z = fid.variables['elevation']
    506                     #print z[:]
    507                     #print A-z[:]
    508 
    509 
    510                     A = choose(A-z[:] >= self.minimum_storable_height,
    511                                (z[:], A))
    512                     stage[i,:] = A.astype(self.precision)
    513                 elif name == 'xmomentum':
    514                     xmomentum[i,:] = A.astype(self.precision)
    515                 elif name == 'ymomentum':
    516                     ymomentum[i,:] = A.astype(self.precision)
    517 
    518                 #As in....
    519             #eval( name + '[i,:] = A.astype(self.precision)' )
    520             #FIXME: But we need a UNIT test for that before refactoring
     524               
     525                # Get ymomentum
     526                Q = domain.quantities['ymomentum']
     527                ymomentum, _ = Q.get_vertex_values(xy = False,
     528                                          precision = self.precision)
     529                self.writer.quantities(fid,
     530                                                 time=self.domain.time,
     531                                                 precision=self.precision,
     532                                                 stage=stage,
     533                                                 xmomentum=xmomentum,
     534                                                 ymomentum=ymomentum)
     535            else:
     536                # This is producing a sww that is not standard.
     537                #Store time
     538                time[i] = self.domain.time
     539               
     540                for name in names:
     541                    # Get quantity
     542                    Q = domain.quantities[name]
     543                    A,V = Q.get_vertex_values(xy = False,
     544                                          precision = self.precision)
     545
     546                    #FIXME: Make this general (see below)
     547                    if name == 'stage':
     548                        z = fid.variables['elevation']
     549                        A = choose(A-z[:] >= self.minimum_storable_height,
     550                                   (z[:], A))
     551                        stage[i,:] = A.astype(self.precision)
     552                    elif name == 'xmomentum':
     553                        xmomentum[i,:] = A.astype(self.precision)
     554                    elif name == 'ymomentum':
     555                        ymomentum[i,:] = A.astype(self.precision)
     556
     557                   #As in....
     558                   #eval( name + '[i,:] = A.astype(self.precision)' )
     559                   #FIXME: But we need a UNIT test for that before
     560                   # refactoring
    521561
    522562
     
    14711511    """
    14721512
    1473     #FIXME: Can this be written feasibly using write_pts?
    1474 
    14751513    import os
    14761514    from Scientific.IO.NetCDF import NetCDFFile
     
    15541592        ptsname = basename_out + '.pts'
    15551593
    1556     #FIXME (DSG-ON): use loadASCII export_points_file
    1557     if verbose: print 'Store to NetCDF file %s' %ptsname
    1558     # NetCDF file definition
    1559     outfile = NetCDFFile(ptsname, 'w')
    1560 
    1561     #Create new file
    1562     outfile.institution = 'Geoscience Australia'
    1563     outfile.description = 'NetCDF pts format for compact and portable ' +\
    1564                           'storage of spatial point data derived from HEC-RAS'
    1565 
    1566     #Georeferencing
    1567     outfile.zone = zone
    1568     outfile.xllcorner = 0.0
    1569     outfile.yllcorner = 0.0
    1570     outfile.false_easting = 500000    #FIXME: Use defaults from e.g. config.py
    1571     outfile.false_northing = 1000000
    1572 
    1573     outfile.projection = projection
    1574     outfile.datum = datum
    1575     outfile.units = units
    1576 
    1577 
    1578     outfile.createDimension('number_of_points', len(points))
    1579     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
    1580 
    1581     # variable definitions
    1582     outfile.createVariable('points', Float, ('number_of_points',
    1583                                              'number_of_dimensions'))
    1584     outfile.createVariable('elevation', Float, ('number_of_points',))
    1585 
    1586     # Get handles to the variables
    1587     outfile.variables['points'][:] = points
    1588     outfile.variables['elevation'][:] = elevation
    1589 
    1590     outfile.close()
    1591 
     1594    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
     1595    geo = Geospatial_data(points, {"elevation":elevation},
     1596                          verbose=verbose, geo_reference=geo_ref)
     1597    geo.export_points_file(ptsname)
    15921598
    15931599
     
    26572663    #Create new file
    26582664    starttime = times[0]
    2659     write_sww_header(outfile, times, number_of_volumes,
    2660                      number_of_points, description=description)
     2665    sww = Write_sww()
     2666    sww.header(outfile, times, number_of_volumes,
     2667                         number_of_points, description=description,
     2668                         verbose=verbose)
    26612669
    26622670    #Store
     
    27492757        x = outfile.variables['x'][:]
    27502758        y = outfile.variables['y'][:]
    2751         times = outfile.variables['time'][:]
    27522759        print '------------------------------------------------'
    27532760        print 'Statistics of output file:'
     
    30113018    other_quantities.remove('z')
    30123019    other_quantities.remove('volumes')
     3020    try:
     3021        other_quantities.remove('stage_range')
     3022        other_quantities.remove('xmomentum_range')
     3023        other_quantities.remove('ymomentum_range')
     3024        other_quantities.remove('elevation_range')
     3025    except:
     3026        pass
     3027       
    30133028
    30143029    conserved_quantities.remove('time')
     
    31403155
    31413156def weed(coordinates,volumes,boundary = None):
    3142     if type(coordinates)=='array':
     3157    if type(coordinates)==ArrayType:
    31433158        coordinates = coordinates.tolist()
    3144     if type(volumes)=='array':
     3159    if type(volumes)==ArrayType:
    31453160        volumes = volumes.tolist()
    31463161
     
    43674382    It will be the origin of the sww file. This shouldn't be used,
    43684383    since all of anuga should be able to handle an arbitary origin.
     4384    The mux point info is NOT relative to this origin.
    43694385
    43704386
     
    44374453    mesh = Mesh()
    44384454    mesh.add_vertices(points_utm)
    4439     mesh.auto_segment()
     4455    mesh.auto_segment(smooth_indents=True, expand_pinch=True)
     4456    # To try and avoid alpha shape 'hugging' too much
     4457    mesh.auto_segment( mesh.shape.get_alpha()*1.1 )
    44404458    if hole_points_UTM is not None:
    44414459        point = ensure_absolute(hole_points_UTM)
     
    44854503    # For a different way of doing this, check out tsh2sww
    44864504    # work out sww_times and the index range this covers
    4487     write_sww_header(outfile, times, len(volumes), len(points_utm))
     4505    sww = Write_sww()
     4506    sww.header(outfile, times, len(volumes), len(points_utm),
     4507                         verbose=verbose)
    44884508    outfile.mean_stage = mean_stage
    44894509    outfile.zscale = zscale
    4490     write_sww_triangulation(outfile, points_utm, volumes,
    4491                             elevation, zone,  origin=origin, verbose=verbose)
     4510
     4511    sww.triangulation(outfile, points_utm, volumes,
     4512                            elevation, zone,  new_origin=origin,
     4513                            verbose=verbose)
    44924514   
    44934515    if verbose: print 'Converting quantities'
    44944516    j = 0
    44954517    # Read in a time slice from each mux file and write it to the sww file
    4496     for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
     4518    for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
    44974519        if j >= mux_times_start_i and j < mux_times_fin_i:
    4498             write_sww_time_slice(outfile, HA, UA, VA,  elevation,
    4499                                  j - mux_times_start_i,
    4500                                  mean_stage=mean_stage, zscale=zscale,
    4501                                  verbose=verbose)
     4520            stage = zscale*ha + mean_stage
     4521            h = stage - elevation
     4522            xmomentum = ua*h
     4523            ymomentum = -1*va*h # -1 since in mux files south is positive.
     4524            sww.quantities(outfile,
     4525                                     slice_index=j - mux_times_start_i,
     4526                                     verbose=verbose,
     4527                                     stage=stage,
     4528                                     xmomentum=xmomentum,
     4529                                     ymomentum=ymomentum)
    45024530        j += 1
     4531    if verbose: sww.verbose_quantities(outfile)
    45034532    outfile.close()
    45044533    #
    45054534    # Do some conversions while writing the sww file
     4535
     4536   
    45064537def mux2sww_time(mux_times, mint, maxt):
    45074538    """
     
    45224553    return mux_times_start_i, mux_times_fin_i
    45234554
    4524 def write_sww_header(outfile, times, number_of_volumes, number_of_points, description='Converted from XXX'):
    4525     """
    4526     outfile - the name of the file that will be written
    4527     times - A list of the time slice times
    4528     number_of_volumes - the number of triangles
    4529     """
    4530     times = ensure_numeric(times)
    4531    
    4532     number_of_times = len(times)
    4533    
    4534     outfile.institution = 'Geoscience Australia'
    4535     outfile.description = 'Converted from XXX'
    4536 
    4537     #For sww compatibility
    4538     outfile.smoothing = 'Yes'
    4539     outfile.order = 1
    4540 
    4541     #Start time in seconds since the epoch (midnight 1/1/1970)
    4542     if len(times) == 0:
    4543         outfile.starttime = starttime = 0
    4544     else:
    4545         outfile.starttime = starttime = times[0]
    4546     times = times - starttime  #Store relative times
    4547 
    4548     # dimension definitions
    4549     outfile.createDimension('number_of_volumes', number_of_volumes)
    4550     outfile.createDimension('number_of_vertices', 3)
    4551     outfile.createDimension('number_of_points', number_of_points)
    4552     outfile.createDimension('number_of_timesteps', number_of_times)
    4553 
    4554     # variable definitions
    4555     outfile.createVariable('x', precision, ('number_of_points',))
    4556     outfile.createVariable('y', precision, ('number_of_points',))
    4557     outfile.createVariable('elevation', precision, ('number_of_points',))
    4558 
    4559     #FIXME: Backwards compatibility
    4560     outfile.createVariable('z', precision, ('number_of_points',))
    4561     #################################
    4562 
    4563     outfile.createVariable('volumes', Int, ('number_of_volumes',
    4564                                             'number_of_vertices'))
    4565 
    4566     outfile.createVariable('time', precision,
    4567                            ('number_of_timesteps',))
    4568 
    4569     outfile.createVariable('stage', precision,
    4570                            ('number_of_timesteps',
    4571                             'number_of_points'))
    4572 
    4573     outfile.createVariable('xmomentum', precision,
    4574                            ('number_of_timesteps',
    4575                             'number_of_points'))
    4576 
    4577     outfile.createVariable('ymomentum', precision,
    4578                            ('number_of_timesteps',
    4579                             'number_of_points'))
    4580    
    4581     outfile.variables['time'][:] = times    #Store time relative
    4582 
    4583 
    4584 def write_sww_triangulation(outfile, points_utm, volumes,
    4585                             elevation, zone, origin=None, verbose=False):
    4586 
    4587     number_of_points = len(points_utm)   
    4588     volumes = array(volumes)
    4589 
    4590     if origin is None:
    4591         origin = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
    4592     geo_ref = write_NetCDF_georeference(origin, outfile)
    4593    
    4594 #     if origin is not None:
    4595 #         if isinstance(origin, Geo_reference):
    4596 #             geo_ref = origin
    4597 #         else:
    4598 #             geo_ref = apply(Geo_reference,origin)
    4599 #     else:
    4600 #         geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
    4601 #     #geo_ref = Geo_reference(zone,0,0)
    4602 #     geo_ref.write_NetCDF(outfile)
    4603    
    4604 
    4605     # This will put the geo ref in the middle
    4606     #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
    4607     x =  points_utm[:,0]
    4608     y =  points_utm[:,1]
    4609     z = outfile.variables['z'][:]
    4610     #FIXME the w.r.t time is lost now..
    4611     if verbose:
     4555
     4556class Write_sww:
     4557    from anuga.shallow_water.shallow_water_domain import Domain
     4558    sww_quantities = Domain.conserved_quantities
     4559    RANGE = '_range'
     4560    def __init__(self):
     4561        pass
     4562   
     4563    def header(self,outfile, times, number_of_volumes,
     4564                         number_of_points,
     4565                         description='Converted from XXX',
     4566                         smoothing=True,
     4567                         order=1, verbose=False):
     4568        """
     4569        outfile - the name of the file that will be written
     4570        times - A list of the time slice times OR a start time
     4571        Note, if a list is given the info will be made relative.
     4572        number_of_volumes - the number of triangles
     4573        """
     4574   
     4575        outfile.institution = 'Geoscience Australia'
     4576        outfile.description = description
     4577
     4578        #For sww compatibility
     4579        if smoothing is True:
     4580            # Smoothing to be depreciated
     4581            outfile.smoothing = 'Yes'
     4582            outfile.vertices_are_stored_uniquely = 'False'
     4583        else:
     4584            # Smoothing to be depreciated
     4585            outfile.smoothing = 'No'
     4586            outfile.vertices_are_stored_uniquely = 'True'
     4587        outfile.order = order
     4588
     4589        try:
     4590            revision_number = get_revision_number()
     4591        except:
     4592            revision_number = None
     4593        # writing a string so None can be written   
     4594        outfile.revision_number = str(revision_number)
     4595        #times - A list or array of the time slice times OR a start time
     4596        #times = ensure_numeric(times)
     4597        #Start time in seconds since the epoch (midnight 1/1/1970)
     4598
     4599        # this is being used to seperate one number from a list.
     4600        # what it is actually doing is sorting lists from numeric arrays.
     4601        if type(times) is list or type(times) is ArrayType: 
     4602            number_of_times = len(times)
     4603            times = ensure_numeric(times) 
     4604            if number_of_times == 0:
     4605                starttime = 0
     4606            else:
     4607                starttime = times[0]
     4608                times = times - starttime  #Store relative times
     4609        else:
     4610            number_of_times = 0
     4611            starttime = times
     4612        #print "times",times
     4613        #print "starttime", starttime
     4614        outfile.starttime = starttime
     4615        # dimension definitions
     4616        outfile.createDimension('number_of_volumes', number_of_volumes)
     4617        outfile.createDimension('number_of_vertices', 3)
     4618        outfile.createDimension('numbers_in_range', 2)
     4619   
     4620        if smoothing is True:
     4621            outfile.createDimension('number_of_points', number_of_points)
     4622       
     4623            # FIXME(Ole): This will cause sww files for paralle domains to
     4624            # have ghost nodes stored (but not used by triangles).
     4625            # To clean this up, we have to change get_vertex_values and
     4626            # friends in quantity.py (but I can't be bothered right now)
     4627        else:
     4628            outfile.createDimension('number_of_points', 3*number_of_volumes)
     4629        outfile.createDimension('number_of_timesteps', number_of_times)
     4630
     4631        # variable definitions
     4632        outfile.createVariable('x', precision, ('number_of_points',))
     4633        outfile.createVariable('y', precision, ('number_of_points',))
     4634        outfile.createVariable('elevation', precision, ('number_of_points',))
     4635        q = 'elevation'
     4636        outfile.createVariable(q+Write_sww.RANGE, precision,
     4637                               ('numbers_in_range',))
     4638        # The values are initally filled with large (10e+36) numbers.
     4639        # I'm relying on this feature.  Maybe I shouldn't?
     4640        outfile.variables[q+Write_sww.RANGE][1] = \
     4641                          -1*outfile.variables[q+Write_sww.RANGE][1]
     4642
     4643        #FIXME: Backwards compatibility
     4644        outfile.createVariable('z', precision, ('number_of_points',))
     4645        #################################
     4646
     4647        outfile.createVariable('volumes', Int, ('number_of_volumes',
     4648                                                'number_of_vertices'))
     4649       
     4650        outfile.createVariable('time', precision,
     4651                               ('number_of_timesteps',))
     4652       
     4653        for q in Write_sww.sww_quantities:
     4654            outfile.createVariable(q, precision,
     4655                                   ('number_of_timesteps',
     4656                                    'number_of_points')) 
     4657            outfile.createVariable(q+Write_sww.RANGE, precision,
     4658                                   ('numbers_in_range',))
     4659            # Initialising the max value to a very small number.
     4660            # It assumes that netcdf initialises it to a very large number
     4661            outfile.variables[q+Write_sww.RANGE][1] = \
     4662                -outfile.variables[q+Write_sww.RANGE][1]
     4663        outfile.variables['time'][:] = times    #Store time relative
     4664           
     4665        if verbose:
     4666            print '------------------------------------------------'
     4667            print 'Statistics:'
     4668            print '    t in [%f, %f], len(t) == %d'\
     4669                  %(min(times.flat), max(times.flat), len(times.flat))
     4670
     4671       
     4672    def triangulation(self,outfile, points_utm, volumes,
     4673                                elevation, zone=None, new_origin=None,
     4674                                points_georeference=None, verbose=False):
     4675        """
     4676       
     4677        new_origin - qa georeference that the points can be set to. (Maybe
     4678        do this before calling this function.)
     4679       
     4680        points_utm - currently a list or array of the points in UTM.
     4681        points_georeference - the georeference of the points_utm
     4682       
     4683        How about passing new_origin and current_origin.
     4684        If you get both, do a convertion from the old to the new.
     4685       
     4686        If you only get new_origin, the points are absolute, convert to relative
     4687       
     4688        if you only get the current_origin the points are relative, store
     4689        as relative.
     4690       
     4691        if you get no georefs create a new georef based on the minimums of
     4692        points_utm.  (Another option would be to default to absolute)
     4693       
     4694        Yes, and this is done in another part of the code.
     4695        Probably geospatial.
     4696       
     4697        If you don't supply either geo_refs, then supply a zone. If not
     4698        the default zone will be used.
     4699       
     4700       
     4701        precon
     4702       
     4703        header has been called.
     4704        """
     4705       
     4706        number_of_points = len(points_utm)   
     4707        volumes = array(volumes)
     4708
     4709        # given the two geo_refs and the points, do the stuff
     4710        # described in the method header
     4711        # if this is needed else where, pull out as a function
     4712        points_georeference = ensure_geo_reference(points_georeference)
     4713        new_origin = ensure_geo_reference(new_origin)
     4714        if new_origin is None and points_georeference is not None:
     4715            points = points_utm
     4716            geo_ref = points_georeference
     4717        else:
     4718            if new_origin is None:
     4719                new_origin = Geo_reference(zone,min(points_utm[:,0]),
     4720                                           min(points_utm[:,1]))
     4721            points = new_origin.change_points_geo_ref(points_utm,
     4722                                                      points_georeference)
     4723            geo_ref = new_origin
     4724
     4725        # At this stage I need a georef and points
     4726        # the points are relative to the georef
     4727        geo_ref.write_NetCDF(outfile)
     4728   
     4729        # This will put the geo ref in the middle
     4730        #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
     4731        x =  points[:,0]
     4732        y =  points[:,1]
     4733        z = outfile.variables['z'][:]
     4734   
     4735        if verbose:
     4736            print '------------------------------------------------'
     4737            print 'More Statistics:'
     4738            print '  Extent (/lon):'
     4739            print '    x in [%f, %f], len(lat) == %d'\
     4740                  %(min(x), max(x),
     4741                    len(x))
     4742            print '    y in [%f, %f], len(lon) == %d'\
     4743                  %(min(y), max(y),
     4744                    len(y))
     4745            print '    z in [%f, %f], len(z) == %d'\
     4746                  %(min(elevation), max(elevation),
     4747                    len(elevation))
     4748            print 'geo_ref: ',geo_ref
     4749            print '------------------------------------------------'
     4750           
     4751        #z = resize(bath_grid,outfile.variables['z'][:].shape)
     4752        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
     4753        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
     4754        outfile.variables['z'][:] = elevation
     4755        outfile.variables['elevation'][:] = elevation  #FIXME HACK
     4756        outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
     4757
     4758        q = 'elevation'
     4759        # This updates the _range values
     4760        outfile.variables[q+Write_sww.RANGE][0] = min(elevation)
     4761        outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
     4762
     4763    def quantities(self, outfile, precision=Float,
     4764                             slice_index=None, time=None,
     4765                             verbose=False, **quant):
     4766        """
     4767        Write the quantity info.
     4768
     4769        **quant is extra keyword arguments passed in. These must be
     4770          the sww quantities, currently; stage, xmomentum, ymomentum.
     4771       
     4772        if the time array is already been built, use the slice_index
     4773        to specify the index.
     4774       
     4775        Otherwise, use time to increase the time dimension
     4776
     4777        Maybe make this general, but the viewer assumes these quantities,
     4778        so maybe we don't want it general - unless the viewer is general
     4779       
     4780        precon
     4781        triangulation and
     4782        header have been called.
     4783        """
     4784
     4785        if time is not None:
     4786            file_time = outfile.variables['time']
     4787            slice_index = len(file_time)
     4788            file_time[slice_index] = time   
     4789
     4790        # write the conserved quantities from Doamin.
     4791        # Typically stage,  xmomentum, ymomentum
     4792        # other quantities will be ignored, silently.
     4793        for q in Write_sww.sww_quantities:
     4794            if not quant.has_key(q):
     4795                msg = 'SWW file can not write quantity %s' %q
     4796                raise NewQuantity, msg
     4797            else:
     4798                q_values = quant[q]
     4799                outfile.variables[q][slice_index] = q_values.astype(precision)
     4800
     4801                # This updates the _range values
     4802                q_range = outfile.variables[q+Write_sww.RANGE][:]
     4803                q_values_min = min(q_values)
     4804                if q_values_min < q_range[0]:
     4805                    outfile.variables[q+Write_sww.RANGE][0] = q_values_min
     4806                q_values_max = max(q_values)
     4807                if q_values_max > q_range[1]:
     4808                    outfile.variables[q+Write_sww.RANGE][1] = q_values_max
     4809
     4810    def verbose_quantities(self, outfile):
    46124811        print '------------------------------------------------'
    46134812        print 'More Statistics:'
    4614         print '  Extent (/lon):'
    4615         print '    x in [%f, %f], len(lat) == %d'\
    4616               %(min(x), max(x),
    4617                 len(x))
    4618         print '    y in [%f, %f], len(lon) == %d'\
    4619               %(min(y), max(y),
    4620                 len(y))
    4621         print '    z in [%f, %f], len(z) == %d'\
    4622               %(min(elevation), max(elevation),
    4623                 len(elevation))
    4624         print 'geo_ref: ',geo_ref
     4813        for q in Write_sww.sww_quantities:
     4814            print '  %s in [%f, %f]' %(q
     4815                                       ,outfile.variables[q+Write_sww.RANGE][0]
     4816                                       ,outfile.variables[q+Write_sww.RANGE][1]
     4817                                       )
    46254818        print '------------------------------------------------'
    4626 
    4627     #z = resize(bath_grid,outfile.variables['z'][:].shape)
    4628     outfile.variables['x'][:] = points_utm[:,0] - geo_ref.get_xllcorner()
    4629     outfile.variables['y'][:] = points_utm[:,1] - geo_ref.get_yllcorner()
    4630     outfile.variables['z'][:] = elevation
    4631     outfile.variables['elevation'][:] = elevation  #FIXME HACK
    4632     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
    4633 
    4634 
    4635 def write_sww_time_slices(outfile, has, uas, vas, elevation,
     4819       
     4820def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation,
    46364821                         mean_stage=0, zscale=1,
    46374822                         verbose=False):   
     
    46514836        ymomentum[j] = -1*va*h  #  -1 since in mux files south is positive.
    46524837        j += 1
    4653 
    4654 def write_sww_time_slice(outfile, ha, ua, va, elevation, slice_index,
    4655                          mean_stage=0, zscale=1,
    4656                          verbose=False):   
    4657     #Time stepping
    4658     stage = outfile.variables['stage']
    4659     xmomentum = outfile.variables['xmomentum']
    4660     ymomentum = outfile.variables['ymomentum']
    4661    
    4662     w = zscale*ha + mean_stage
    4663     stage[slice_index] = w
    4664     h = w - elevation
    4665     xmomentum[slice_index] = ua*h
    4666     ymomentum[slice_index] = -1*va*h # -1 since in mux files south is positive.
    4667 
     4838   
    46684839def urs2txt(basename_in, location_index=None):
    46694840    """
Note: See TracChangeset for help on using the changeset viewer.