Ignore:
Timestamp:
Jun 23, 2008, 4:26:45 PM (15 years ago)
Author:
jakeman
Message:

allow urs2sts to accept multiple sources

File:
1 edited

Legend:

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

    r5409 r5412  
    48384838NORTH_VELOCITY_MUX2_LABEL =  '_velocity-n-mux2'   
    48394839
    4840 def read_mux2_py(filenames,weights=None):
     4840def read_mux2_py(filenames,weights):
    48414841
    48424842    from Numeric import ones,Float,compress,zeros,arange
    48434843    from urs_ext import read_mux2
    48444844
    4845     if weights is None:
    4846         weights=ones(len(filenames),Float)/len(filenames) #default 1/numSrc
     4845    numSrc=len(filenames)
    48474846
    48484847    file_params=-1*ones(3,Float)#[nsta,dt,nt]
    48494848    write=0 #if true write txt files to current directory as well
    4850     data=read_mux2(1,filenames,weights,file_params,write)
     4849    data=read_mux2(numSrc,filenames,weights,file_params,write)
    48514850
    48524851    msg='File parameter values were not read in correctly from c file'
     
    48734872    elevation=zeros(data.shape[0],Float)
    48744873    quantity=zeros((data.shape[0],data.shape[1]-OFFSET),Float)
     4874    starttime=1e16
    48754875    for i in range(0,data.shape[0]):
    48764876        latitudes[i]=data[i][data.shape[1]-OFFSET]
     
    48784878        elevation[i]=-data[i][data.shape[1]-OFFSET+2]
    48794879        quantity[i]=data[i][:-OFFSET]
    4880 
    4881     return times, latitudes, longitudes, elevation, quantity
     4880        starttime=min(dt*data[i][data.shape[1]-OFFSET+3],starttime)
     4881       
     4882    return times, latitudes, longitudes, elevation, quantity, starttime
    48824883
    48834884def mux2sww_time(mux_times, mint, maxt):
     
    49004901
    49014902
    4902 def urs2sts(basename_in, basename_out = None, verbose = False, origin = None,
     4903def urs2sts(basename_in, basename_out = None, weights=None,
     4904            verbose = False, origin = None,
    49034905            mean_stage=0.0,zscale=1.0,
    49044906            minlat = None, maxlat = None,
     
    49064908    """Convert URS mux2 format for wave propagation to sts format
    49074909
    4908     Specify only basename_in and read files of the form
    4909     out_waveheight-z-mux2
    4910 
    49114910    Also convert latitude and longitude to UTM. All coordinates are
    49124911    assumed to be given in the GDA94 datum
     
    49184917    from Scientific.IO.NetCDF import NetCDFFile
    49194918    from Numeric import Float, Int, Int32, searchsorted, zeros, array
    4920     from Numeric import allclose, around
     4919    from Numeric import allclose, around,ones,Float
     4920    from types import ListType,StringType
     4921    from operator import __and__
     4922   
     4923    if not isinstance(basename_in, ListType):
     4924        if verbose: print 'Reading single source'
     4925        basename_in=[basename_in]
     4926
     4927    # Check that basename is a list of strings
     4928    if not reduce(__and__, map(lambda z:isinstance(z,StringType),basename_in)):
     4929        msg= 'basename_in must be a string or list of strings'
     4930        raise Exception, msg
     4931
     4932    # Find the number of sources to be used
     4933    numSrc=len(basename_in)
     4934
     4935    # A weight must be specified for each source
     4936    if weights is None:
     4937        weights=ones(numSrc,Float)/numSrc
     4938    else:
     4939        weights = ensure_numeric(weights)
     4940        msg = 'When combining multiple sources must specify a weight for '+\
     4941              'mux2 source file'
     4942        assert len(weights)== numSrc, msg
    49214943
    49224944    precision = Float
     
    49384960
    49394961    if basename_out is None:
    4940         stsname = basename_in + '.sts'
     4962        stsname = basename_in[0] + '.sts'
    49414963    else:
    49424964        stsname = basename_out + '.sts'
    49434965
    4944     files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
    4945                 basename_in + EAST_VELOCITY_MUX2_LABEL,
    4946                 basename_in + NORTH_VELOCITY_MUX2_LABEL]
     4966    files_in=[[],[],[]]
     4967    for files in basename_in:
     4968        files_in[0].append(files + WAVEHEIGHT_MUX2_LABEL),
     4969        files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL)
     4970        files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL)
     4971       
     4972    #files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
     4973    #            basename_in + EAST_VELOCITY_MUX2_LABEL,
     4974    #            basename_in + NORTH_VELOCITY_MUX2_LABEL]
     4975   
    49474976    quantities = ['HA','UA','VA']
    49484977
    4949     for file_in in files_in:
    4950         if (os.access(file_in, os.F_OK) == 0):
    4951             msg = 'File %s does not exist or is not accessible' %file_in
    4952             raise IOError, msg
     4978    # For each source file check that there exists three files ending with:
     4979    # WAVEHEIGHT_MUX2_LABEL,
     4980    # EAST_VELOCITY_MUX2_LABEL, and
     4981    # NORTH_VELOCITY_MUX2_LABEL
     4982    for i in range(len(quantities)):
     4983        for file_in in files_in[i]:
     4984            if (os.access(file_in, os.F_OK) == 0):
     4985                msg = 'File %s does not exist or is not accessible' %file_in
     4986                raise IOError, msg
    49534987
    49544988    #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
     
    49564990    if (verbose): print 'reading mux2 file'
    49574991    mux={}
    4958     for quantity, file in map(None, quantities, files_in):
    4959         times_urs, latitudes_urs, longitudes_urs, elevation, mux[quantity] = read_mux2_py([file])
     4992    for i,quantity in enumerate(quantities):
     4993        times_urs, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights)
    49604994        if quantity!=quantities[0]:
    49614995            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
     
    49644998            assert allclose(longitudes_urs,longitudes_urs_old),msg
    49654999            assert allclose(elevation,elevation_old),msg
     5000            assert allclose(starttime,starttime_old)
    49665001        times_urs_old=times_urs
    49675002        latitudes_urs_old=latitudes_urs
    49685003        longitudes_urs_old=longitudes_urs
    49695004        elevation_old=elevation
     5005        starttime_old=starttime
    49705006       
    49715007    if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
    4972         if verbose: print 'Cliiping urs data'
     5008        if verbose: print 'Cliping urs data'
    49735009        latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)
    49745010        longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)
     
    49945030   
    49955031    # Create new file
    4996     starttime = times[0]
     5032    #starttime = times[0]
    49975033    sts = Write_sts()
    4998     sts.store_header(outfile, times,number_of_points, description=description,
     5034    sts.store_header(outfile, times+starttime,number_of_points, description=description,
    49995035                     verbose=verbose,sts_precision=Float)
    50005036   
Note: See TracChangeset for help on using the changeset viewer.