Changeset 5534


Ignore:
Timestamp:
Jul 18, 2008, 3:14:35 PM (11 years ago)
Author:
ole
Message:

More cleanup prior to implementation of permutation array

File:
1 edited

Legend:

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

    r5485 r5534  
    48414841NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'   
    48424842
    4843 def read_mux2_py(filenames,weights,verbose=False):
     4843def read_mux2_py(filenames,
     4844                 weights,
     4845                 permutation=None,
     4846                 verbose=False):
     4847    """Access the mux files specified in the filenames list. Combine the
     4848       data found therin as a weighted linear sum as specifed by the weights.
     4849       If permutation is None extract timeseries data for all gauges within the
     4850       files.
     4851         
     4852       Input:
     4853           filenames:   List of filenames specifiying the file containing the
     4854                        timeseries data (mux2 format) for each source
     4855           weights:     Weights associated with each source
     4856           permutation: Specifies the gauge numbers that for which data is to be
     4857                        extracted
     4858    """
    48444859
    48454860    from Numeric import ones,Float,compress,zeros,arange
     
    48564871        verbose=0
    48574872       
     4873    # Call underlying C implementation urs2sts_ext.c   
    48584874    data=read_mux2(numSrc,filenames,weights,file_params,verbose)
    48594875
     
    49104926
    49114927
    4912 def urs2sts(basename_in, basename_out = None, weights=None,
    4913             verbose = False, origin = None,
    4914             mean_stage=0.0, zscale=1.0,
    4915             ordering_filename = None,
    4916             minlat = None, maxlat = None,
    4917             minlon = None, maxlon = None):
     4928def urs2sts(basename_in, basename_out=None,
     4929            weights=None,
     4930            verbose=False,
     4931            origin=None,
     4932            mean_stage=0.0,
     4933            zscale=1.0,
     4934            ordering_filename=None):
    49184935    """Convert URS mux2 format for wave propagation to sts format
    49194936
     
    49244941    UTM coordinates (zone, easting, northing)
    49254942   
    4926     input:
     4943    inputs:
     4944       
    49274945    basename_in: list of source file prefixes
    49284946   
    4929     These are combined with the extensions:
    4930     WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage
    4931     EAST_VELOCITY_MUX2_LABEL =  '-e-mux2' xmomentum
    4932     NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2' and ymomentum
    4933    
    4934     to create a 2D list of mux2 file. The rows are associated with each quantity and must have the above extensions
    4935     the columns are the list of file prefixes.
    4936    
    4937     ordering: a .txt file name specifying which mux2 gauge points are sto be stored. This is indicated by the index of
    4938     the gauge in the ordering file.
    4939    
    4940     ordering file format:
    4941     header
    4942     index,longitude,latitude
    4943    
    4944     header='index,longitude,latitude\n'
    4945     If ordering is None all points are taken in the order they appear in the mux2 file subject to rectangular clipping box (see below).
    4946    
    4947    
    4948     minlat, minlon, maxlat, maxlon define a rectangular clipping box and can be used to extract gauge information in a rectangular box. This is useful for extracting
    4949     information at gauge points not on the boundary.
    4950     if ordering_filename is specified clipping box cannot be used
     4947        These are combined with the extensions:
     4948        WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage
     4949        EAST_VELOCITY_MUX2_LABEL =  '-e-mux2' xmomentum
     4950        NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2' and ymomentum
     4951   
     4952        to create a 2D list of mux2 file. The rows are associated with each
     4953        quantity and must have the above extensions
     4954        the columns are the list of file prefixes.
     4955   
     4956    ordering: a .txt file name specifying which mux2 gauge points are
     4957              to be stored. This is indicated by the index of the gauge
     4958              in the ordering file.
     4959   
     4960              ordering file format:
     4961              1st line:    'index,longitude,latitude\n'
     4962              other lines: index,longitude,latitude
     4963             
     4964
     4965              If ordering is None all points are taken in the order they
     4966              appear in the mux2 file.
    49514967   
    49524968   
    49534969    output:
    4954     baename_out: name of sts file in which mux2 data is stored.
     4970      basename_out: name of sts file in which mux2 data is stored.
    49554971   
    49564972   
     
    49684984
    49694985    # Check that basename is a list of strings
    4970     if not reduce(__and__, map(lambda z:isinstance(z,StringType),basename_in)):
     4986    if not reduce(__and__, map(lambda z:isinstance(z,StringType), basename_in)):
    49714987        msg= 'basename_in must be a string or list of strings'
    49724988        raise Exception, msg
     
    49774993    # A weight must be specified for each source
    49784994    if weights is None:
     4995        # Default is equal weighting
    49794996        weights=ones(numSrc,Float)/numSrc
    49804997    else:
     
    49845001        assert len(weights)== numSrc, msg
    49855002
    4986     precision = Float
    4987 
    4988     if ordering_filename is not None:
    4989         msg = 'If ordering_filename is specified,'
    4990         msg += ' rectangular clipping box cannot be specified as well. '
    4991         assert minlat is None and maxlat is None and\
    4992             minlon is None and maxlon is None, msg
    4993    
    4994    
    4995     msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
    4996     if minlat != None:
    4997         assert -90 < minlat < 90 , msg
    4998     if maxlat != None:
    4999         assert -90 < maxlat < 90 , msg
    5000         if minlat != None:
    5001             assert maxlat > minlat
    5002     if minlon != None:
    5003         assert -180 < minlon < 180 , msg
    5004     if maxlon != None:
    5005         assert -180 < maxlon < 180 , msg
    5006         if minlon != None:
    5007             assert maxlon > minlon
    5008 
    5009 
    50105003    # Check output filename   
    50115004    if basename_out is None:
    50125005        msg = 'STS filename must be specified'
    50135006        raise Exception, msg
    5014    
    50155007    stsname = basename_out + '.sts'
    50165008
     5009    # Create input filenames from basenames and check their existence
    50175010    files_in=[[],[],[]]
    50185011    for files in basename_in:
     
    50205013        files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL)
    50215014        files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL)
    5022        
    5023     #files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
    5024     #            basename_in + EAST_VELOCITY_MUX2_LABEL,
    5025     #            basename_in + NORTH_VELOCITY_MUX2_LABEL]
    5026    
    5027     quantities = ['HA','UA','VA']
    5028 
    5029     # For each source file check that there exists three files ending with:
    5030     # WAVEHEIGHT_MUX2_LABEL,
    5031     # EAST_VELOCITY_MUX2_LABEL, and
    5032     # NORTH_VELOCITY_MUX2_LABEL
     5015   
     5016    quantities = ['HA','UA','VA'] # Quantity names used in the MUX2 format
    50335017    for i in range(len(quantities)):
    50345018        for file_in in files_in[i]:
     
    50375021                raise IOError, msg
    50385022
    5039     #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
    5040     #for now set x_momentum and y_moentum quantities to zero
     5023    # Read MUX2 files
    50415024    if (verbose): print 'reading mux2 file'
    50425025    mux={}
    5043     for i,quantity in enumerate(quantities):
    5044         # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
    5045         times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
     5026    for i, quantity in enumerate(quantities):
     5027   
     5028        # For each quantity read the associated list of source mux2 file with
     5029        # extention associated with that quantity
     5030        times,\
     5031        latitudes_urs,\
     5032        longitudes_urs,\
     5033        elevation,\
     5034        mux[quantity],\
     5035        starttime = read_mux2_py(files_in[i], weights, verbose=verbose)
     5036   
     5037        # Check that all quantities have consistent time and space information     
    50465038        if quantity!=quantities[0]:
    5047             msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
    5048             assert allclose(times,times_old),msg
    5049             assert allclose(latitudes_urs,latitudes_urs_old),msg
    5050             assert allclose(longitudes_urs,longitudes_urs_old),msg
    5051             assert allclose(elevation,elevation_old),msg
    5052             assert allclose(starttime,starttime_old)
     5039            msg='%s, %s and %s have inconsistent gauge data'\
     5040                %(files_in[0], files_in[1], files_in[2])
     5041            assert allclose(times, times_old), msg
     5042            assert allclose(latitudes_urs, latitudes_urs_old), msg
     5043            assert allclose(longitudes_urs, longitudes_urs_old), msg
     5044            assert allclose(elevation, elevation_old), msg
     5045            assert allclose(starttime, starttime_old), msg
    50535046        times_old=times
    50545047        latitudes_urs_old=latitudes_urs
     
    51045097        assert allclose(latitudes, ref_latitudes), msg
    51055098       
    5106     elif (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
    5107         #FIXME - this branch is probably not working   
    5108         if verbose: print 'Cliping urs data'
    5109         latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)
    5110         longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)
    5111        
    5112         permutation = range(latitudes.shape[0])       
    5113        
    5114         msg = 'Clipping is not done yet'
    5115         raise Exception, msg
    5116        
    51175099    else:
    51185100        latitudes=latitudes_urs
     
    51205102        permutation = range(latitudes.shape[0])               
    51215103       
    5122 
     5104       
     5105       
     5106    # Store timeseries in NetCDF STS file   
    51235107    msg='File is empty and or clipped region not in file region'
    51245108    assert len(latitudes>0),msg
Note: See TracChangeset for help on using the changeset viewer.