Ignore:
Timestamp:
Jul 7, 2008, 3:34:26 PM (14 years ago)
Author:
ole
Message:

Allow urs_ext to store only a subset of gauges in mux file according to a list of ordered points

File:
1 edited

Legend:

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

    r5470 r5473  
    48384838NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'   
    48394839
    4840 def read_mux2_py(filenames,weights,verbose=False):
    4841 
     4840def read_mux2_py(filenames,weights,permutation=None,verbose=False):
     4841    """
     4842    Access the mux files specified in the filenames list. Combine the
     4843    data found therin as a weighted linear sum as specifed by the weights.
     4844    If permutation is None extract timeseries data for all gauges within the
     4845    files.
     4846   
     4847    Input:
     4848    filenames:   List of filenames specifiying the file containing the
     4849                 timeseries data (mux2 format) for each source
     4850    weights:     Weights associated with each source
     4851    permutation: Specifies the gauge numbers that for which data is to be
     4852                 extracted
     4853    """
    48424854    from Numeric import ones,Float,compress,zeros,arange
    48434855    from urs_ext import read_mux2
     
    48534865        verbose=0
    48544866       
    4855     data=read_mux2(numSrc,filenames,weights,file_params,verbose)
     4867    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    48564868
    48574869    msg='File parameter values were not read in correctly from c file'
     
    50345046                raise IOError, msg
    50355047
    5036     #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
    5037     #for now set x_momentum and y_moentum quantities to zero
     5048    if ordering_filename is not None:
     5049        # Read ordering file
     5050       
     5051        try:
     5052            fid=open(ordering_filename, 'r')
     5053            file_header=fid.readline().split(',')
     5054            ordering_lines=fid.readlines()
     5055            fid.close()
     5056        except:
     5057            msg = 'Cannot open %s'%ordering_filename
     5058            raise Exception, msg
     5059
     5060        reference_header = 'index, longitude, latitude\n'
     5061        reference_header_split = reference_header.split(',')
     5062        for i in range(3):
     5063            if not file_header[i] == reference_header_split[i]:
     5064                msg = 'File must contain header: '+reference_header+'\n'
     5065                raise Exception, msg
     5066
     5067        if len(ordering_lines)<2:
     5068            msg = 'File must contain at least two points'
     5069            raise Exception, msg
     5070
     5071   
     5072       
     5073        permutation = [int(line.split(',')[0]) for line in ordering_lines]
     5074    else:
     5075        permutation = None 
     5076
    50385077    if (verbose): print 'reading mux2 file'
    50395078    mux={}
    50405079    for i,quantity in enumerate(quantities):
    50415080        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
    5042         times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
     5081        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
    50435082        if quantity!=quantities[0]:
    50445083            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
     
    50535092        elevation_old=elevation
    50545093        starttime_old=starttime
    5055 
    50565094       
    50575095    if ordering_filename is not None:
    5058         # Read ordering file
    5059        
    5060         try:
    5061             fid=open(ordering_filename, 'r')
    5062             file_header=fid.readline().split(',')
    5063             ordering_lines=fid.readlines()
    5064             fid.close()
    5065         except:
    5066             msg = 'Cannot open %s'%ordering_filename
    5067             raise Exception, msg
    5068 
    5069         reference_header = 'index, longitude, latitude\n'
    5070         reference_header_split = reference_header.split(',')
    5071         for i in range(3):
    5072             if not file_header[i] == reference_header_split[i]:
    5073                 msg = 'File must contain header: '+reference_header+'\n'
    5074                 raise Exception, msg
    5075 
    5076         if len(ordering_lines)<2:
    5077             msg = 'File must contain at least two points'
    5078             raise Exception, msg
    5079 
    5080    
    5081        
    5082         permutation = [int(line.split(',')[0]) for line in ordering_lines]
    5083    
    5084         latitudes = take(latitudes_urs, permutation)
    5085         longitudes = take(longitudes_urs, permutation)
     5096        latitudes = latitudes_urs
     5097        longitudes = longitudes_urs
    50865098       
    50875099        # Self check - can be removed to improve speed
     
    51165128        longitudes=longitudes_urs
    51175129        permutation = range(latitudes.shape[0])               
    5118        
    51195130
    51205131    msg='File is empty and or clipped region not in file region'
     
    51785189    if verbose: print 'Converting quantities'
    51795190    for j in range(len(times)):
    5180         for i, index in enumerate(permutation):
    5181             w = zscale*mux['HA'][index,j] + mean_stage
    5182             h=w-elevation[index]
     5191        for i in range(number_of_points):
     5192            w = zscale*mux['HA'][i,j] + mean_stage
     5193            h=w-elevation[i]
    51835194            stage[j,i] = w
    51845195
    5185             xmomentum[j,i] = mux['UA'][index,j]*h
    5186             ymomentum[j,i] = mux['VA'][index,j]*h
     5196            xmomentum[j,i] = mux['UA'][i,j]*h
     5197            ymomentum[j,i] = mux['VA'][i,j]*h
    51875198
    51885199    outfile.close()
Note: See TracChangeset for help on using the changeset viewer.