Changeset 5412


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

allow urs2sts to accept multiple sources

Location:
anuga_core/source/anuga/shallow_water
Files:
3 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   
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5372 r5412  
    61066106        return base_name, files
    61076107
    6108     def test_read_mux2_py(self):
     6108    def test_read_mux2_pyI(self):
    61096109        """Constant stage,momentum at each gauge
    61106110        """
     
    61316131        weights=ones(1,Float)
    61326132        #ensure that files are indeed mux2 files
    6133         times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
    6134         ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
     6133        times, latitudes, longitudes, elevation, stage, starttime =read_mux2_py([files[0]],weights)
     6134        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
    61356135        msg='ha and ua have different gauge meta data'
    6136         assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
    6137         va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
     6136        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
     6137        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity, starttime_va=read_mux2_py([files[2]],weights)
    61386138        msg='ha and va have different gauge meta data'
    6139         assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
     6139        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
    61406140
    61416141        self.delete_mux(files)
     
    61586158        assert allclose(yvelocity,va)
    61596159
    6160     def test_read_mux2_py2(self):
     6160    def test_read_mux2_pyII(self):
    61616161        """Spatially varing stage
    61626162        """
     
    61886188        weights=ones(1,Float)
    61896189        #ensure that files are indeed mux2 files
    6190         times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
    6191         ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
     6190        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
     6191        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
    61926192        msg='ha and ua have different gauge meta data'
    6193         assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
    6194         va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
     6193        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
     6194        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
    61956195        msg='ha and va have different gauge meta data'
    6196         assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
     6196        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
    61976197
    61986198
     
    62166216        assert allclose(yvelocity,va)
    62176217
    6218     def test_read_mux2_py3(self):
     6218    def test_read_mux2_pyIII(self):
    62196219        """Varying start and finsh times
    62206220        """
     
    62496249        weights=ones(1,Float)
    62506250        #ensure that files are indeed mux2 files
    6251         times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
    6252         ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
     6251        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
     6252        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
    62536253        msg='ha and ua have different gauge meta data'
    6254         assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
    6255         va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
     6254        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
     6255        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
    62566256        msg='ha and va have different gauge meta data'
    6257         assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
     6257        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
    62586258
    62596259        self.delete_mux(files)
     
    62886288        assert allclose(yvelocity,va)
    62896289
    6290     def test_urs2sts(self):
     6290    def test_urs2stsI(self):
     6291        """
     6292        Test single source
     6293        """
    62916294        tide=0
    62926295        time_step_count = 3
     
    63086311        ua=5*ones((n,time_step_count),Float)
    63096312        va=-10*ones((n,time_step_count),Float)
    6310         #-ve added to take into account mux file format where south is positive.
     6313
    63116314        base_name, files = self.write_mux2(lat_long_points,
    63126315                                      time_step_count, time_step,
     
    63956398        fid.close()
    63966399        self.delete_mux(files)
     6400        os.remove(sts_file)
     6401
     6402    def test_urs2stsII(self):
     6403        """
     6404        Test multiple sources
     6405        """
     6406        tide=0
     6407        time_step_count = 3
     6408        time_step = 2
     6409        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     6410        n=len(lat_long_points)
     6411        first_tstep=ones(n,Int)
     6412        first_tstep[0]+=1
     6413        first_tstep[2]+=1
     6414        last_tstep=(time_step_count)*ones(n,Int)
     6415        last_tstep[0]-=1
     6416
     6417        gauge_depth=20*ones(n,Float)
     6418        ha=2*ones((n,time_step_count),Float)
     6419        ha[0]=arange(0,time_step_count)
     6420        ha[1]=arange(time_step_count,2*time_step_count)
     6421        ha[2]=arange(2*time_step_count,3*time_step_count)
     6422        ha[3]=arange(3*time_step_count,4*time_step_count)
     6423        ua=5*ones((n,time_step_count),Float)
     6424        va=-10*ones((n,time_step_count),Float)
     6425   
     6426        base_nameI, filesI = self.write_mux2(lat_long_points,
     6427                                      time_step_count, time_step,
     6428                                      first_tstep, last_tstep,
     6429                                      depth=gauge_depth,
     6430                                      ha=ha,
     6431                                      ua=ua,
     6432                                      va=va)
     6433
     6434        base_nameII, filesII = self.write_mux2(lat_long_points,
     6435                                      time_step_count, time_step,
     6436                                      first_tstep, last_tstep,
     6437                                      depth=gauge_depth,
     6438                                      ha=ha,
     6439                                      ua=ua,
     6440                                      va=va)
     6441
     6442        urs2sts([base_nameI,base_nameII],weights=[1.0,1.0],mean_stage=tide,
     6443                verbose=False)
     6444
     6445        # now I want to check the sts file ...
     6446        sts_file = base_nameI + '.sts'
     6447
     6448        #Let's interigate the sww file
     6449        # Note, the sww info is not gridded.  It is point data.
     6450        fid = NetCDFFile(sts_file)
     6451
     6452        # Make x and y absolute
     6453        x = fid.variables['x'][:]
     6454        y = fid.variables['y'][:]
     6455
     6456        geo_reference = Geo_reference(NetCDFObject=fid)
     6457        points = geo_reference.get_absolute(map(None, x, y))
     6458        points = ensure_numeric(points)
     6459
     6460        x = points[:,0]
     6461        y = points[:,1]
     6462
     6463        #Check that first coordinate is correctly represented       
     6464        #Work out the UTM coordinates for first point
     6465        zone, e, n = redfearn(lat_long_points[0][0], lat_long_points[0][1])
     6466        assert allclose([x[0],y[0]], [e,n])
     6467
     6468        #Check the time vector
     6469        times = fid.variables['time'][:]
     6470
     6471        times_actual = []
     6472        for i in range(time_step_count):
     6473            times_actual.append(time_step * i)
     6474
     6475        assert allclose(ensure_numeric(times),
     6476                        ensure_numeric(times_actual))
     6477
     6478        #Check first value
     6479        stage = fid.variables['stage'][:]
     6480        xmomentum = fid.variables['xmomentum'][:]
     6481        ymomentum = fid.variables['ymomentum'][:]
     6482        elevation = fid.variables['elevation'][:]
     6483
     6484        # Set original data used to write mux file to be zero when gauges are
     6485        #not recdoring
     6486       
     6487        ha[0][0]=0.0
     6488        ha[0][time_step_count-1]=0.0
     6489        ha[2][0]=0.0
     6490        ua[0][0]=0.0
     6491        ua[0][time_step_count-1]=0.0
     6492        ua[2][0]=0.0
     6493        va[0][0]=0.0
     6494        va[0][time_step_count-1]=0.0
     6495        va[2][0]=0.0;
     6496
     6497        # The stage stored in the .sts file should be the sum of the stage
     6498        # in the two mux2 files because both have weights = 1. In this case
     6499        #the mux2 files are the same so stage == 2.0 * ha
     6500        assert allclose(2.0*transpose(ha),stage)  #Meters
     6501
     6502        #Check the momentums - ua
     6503        #momentum = velocity*(stage-elevation)
     6504        # elevation = - depth
     6505        #momentum = velocity_ua *(stage+depth)
     6506
     6507        depth=zeros((len(lat_long_points),time_step_count),Float)
     6508        for i in range(len(lat_long_points)):
     6509            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
     6510            #2.0*ha necessary because using two files with weights=1 are used
     6511
     6512        # The xmomentum stored in the .sts file should be the sum of the ua
     6513        # in the two mux2 files multiplied by the depth.
     6514        assert allclose(2.0*transpose(ua*depth),xmomentum)
     6515
     6516        #Check the momentums - va
     6517        #momentum = velocity*(stage-elevation)
     6518        # elevation = - depth
     6519        #momentum = velocity_va *(stage+depth)
     6520
     6521        # The ymomentum stored in the .sts file should be the sum of the va
     6522        # in the two mux2 files multiplied by the depth.
     6523        assert allclose(2.0*transpose(va*depth),ymomentum)
     6524
     6525        # check the elevation values.
     6526        # -ve since urs measures depth, sww meshers height,
     6527        assert allclose(-elevation, gauge_depth)  #Meters
     6528
     6529        fid.close()
     6530        self.delete_mux(filesI)
     6531        self.delete_mux(filesII)
    63976532        os.remove(sts_file)
    63986533
     
    84148549    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
    84158550#    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
    8416     suite = unittest.makeSuite(Test_Data_Manager,'test')
     8551    suite = unittest.makeSuite(Test_Data_Manager,'test_urs2stsII')
    84178552    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    84188553    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5375 r5412  
    2929   
    3030    Python call:
    31     read_mux2(filenames,weights,numSrc)
     31    read_mux2(numSrc,filenames,weights,file_params,write)
    3232
    3333    NOTE:
     
    8181  }
    8282
    83   if(PyList_Size(filenames) != pyweights->nd){
    84     PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
    85     return NULL;
     83  if(PyList_Size(filenames) != pyweights->dimensions[0]){
     84      PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
     85      return NULL;
    8686  }
    8787
     
    115115  }
    116116
    117   cdata=_read_mux2((int)pyweights->nd,muxFileNameArray,weights,(double*)file_params->data,(int)write);
     117  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)write);
    118118
    119119
     
    200200   /* note starting time */
    201201   time(&start_time);
    202    
    203    
     202
    204203   /* allocate space for the names and the weights and pointers to the data*/   
    205204   wt = (float *) malloc(numSrc*sizeof(float));
     
    274273      mytgs = (struct tgsrwg *)malloc( nsta0*sizeof(struct tgsrwg) );
    275274   } else {
    276      /* FIXME (Ole): What should happen in case the are no source files?*/
     275     /* FIXME (JJ): What should happen in case the are no source files?*/
    277276     /* If we exit here, tests will break */
    278277     // fprintf(stderr, "No source file specified\n");
     
    439438   //free(mytgs0);
    440439   
    441    if(numSrc>1)
    442    {
     440   //if(numSrc>1)
     441   //{
    443442     //free(data);
    444443     //free(mytgs);
    445    }
     444   //}
    446445   //can't free arrays because I only fill array by making pointer not copy
    447446   //for(isrc=0; isrc<numSrc;isrc++)
Note: See TracChangeset for help on using the changeset viewer.