Changeset 5537


Ignore:
Timestamp:
Jul 18, 2008, 4:59:46 PM (11 years ago)
Author:
ole
Message:

Work on permutations in urs_ext.c
It does not work yet, but the array is being passed in and the algorithm
has been described.

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

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

    r5534 r5537  
    48474847    """Access the mux files specified in the filenames list. Combine the
    48484848       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
     4849       If permutation is None or empty extract timeseries data for all gauges within the
    48504850       files.
    48514851         
     
    48714871        verbose=0
    48724872       
     4873       
     4874    #msg = 'I got the permutation vector:' + str(permutation)
     4875    #assert permutation is not None, msg
     4876    if permutation is None:
     4877        permutation = ensure_numeric([], Float)   
     4878       
    48734879    # Call underlying C implementation urs2sts_ext.c   
    4874     data=read_mux2(numSrc,filenames,weights,file_params,verbose)
     4880    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    48754881
    48764882    msg='File parameter values were not read in correctly from c file'
    48774883    assert len(compress(file_params>0,file_params))!=0,msg
    4878     msg='The number of stations specifed in the c array and in the file are inconsitent'
    4879     assert file_params[0]==data.shape[0],msg
     4884   
     4885    msg='The number of stations specifed in the c array and in the file are inconsistent'
     4886    assert file_params[0] >= len(permutation)
     4887   
     4888    msg='The number of stations returned is inconsistent with the requested number'   
     4889    assert len(permutation) == 0 or len(permutation) == data.shape[0], msg
    48804890   
    48814891    nsta=int(file_params[0])
    48824892    msg='Must have at least one station'
    48834893    assert nsta>0,msg
     4894   
    48844895    dt=file_params[1]
    48854896    msg='Must have a postive timestep'
    48864897    assert dt>0,msg
     4898   
    48874899    nt=int(file_params[2])
    48884900    msg='Must have at least one gauge value'
     
    49624974              other lines: index,longitude,latitude
    49634975             
    4964 
    4965               If ordering is None all points are taken in the order they
     4976              If ordering is None or ordering file is empty then
     4977              all points are taken in the order they
    49664978              appear in the mux2 file.
    49674979   
     
    50215033                raise IOError, msg
    50225034
    5023     # Read MUX2 files
    5024     if (verbose): print 'reading mux2 file'
    5025     mux={}
    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     
    5038         if quantity!=quantities[0]:
    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
    5046         times_old=times
    5047         latitudes_urs_old=latitudes_urs
    5048         longitudes_urs_old=longitudes_urs
    5049         elevation_old=elevation
    5050         starttime_old=starttime
    5051 
    5052        
     5035    # Establish permutation array
    50535036    if ordering_filename is not None:
     5037   
    50545038        # Read ordering file
    5055        
    50565039        try:
    50575040            fid=open(ordering_filename, 'r')
     
    50745057            raise Exception, msg
    50755058
    5076    
    5077        
    5078         permutation = [int(line.split(',')[0]) for line in ordering_lines]
    5079    
    5080         latitudes = take(latitudes_urs, permutation)
    5081         longitudes = take(longitudes_urs, permutation)
     5059        permutation = ensure_numeric([int(line.split(',')[0]) for line in ordering_lines])
     5060    else:
     5061        # Use empty array to signify 'all' points
     5062        # We could have used 'None' but it got too hard in the C-code ;-)
     5063        permutation = ensure_numeric([], Float)
     5064
     5065               
     5066    # Read MUX2 files
     5067    if (verbose): print 'reading mux2 file'
     5068    mux={}
     5069    for i, quantity in enumerate(quantities):
     5070   
     5071        # For each quantity read the associated list of source mux2 file with
     5072        # extention associated with that quantity
     5073        times,\
     5074        latitudes,\
     5075        longitudes,\
     5076        elevation,\
     5077        mux[quantity],\
     5078        starttime = read_mux2_py(files_in[i], weights, permutation, verbose=verbose)
     5079   
     5080        # Check that all quantities have consistent time and space information     
     5081        if quantity!=quantities[0]:
     5082            msg='%s, %s and %s have inconsistent gauge data'\
     5083                %(files_in[0], files_in[1], files_in[2])
     5084            assert allclose(times, times_old), msg
     5085            assert allclose(latitudes, latitudes_old), msg
     5086            assert allclose(longitudes, longitudes_old), msg
     5087            assert allclose(elevation, elevation_old), msg
     5088            assert allclose(starttime, starttime_old), msg
     5089        times_old=times
     5090        latitudes_old=latitudes
     5091        longitudes_old=longitudes
     5092        elevation_old=elevation
     5093        starttime_old=starttime
     5094
     5095
     5096       
     5097       
    50825098       
    50835099        # Self check - can be removed to improve speed
    5084         ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]               
    5085         ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]       
    5086        
    5087         msg = 'Longitudes specified in ordering file do not match those found in mux files'
    5088         msg += 'I got %s instead of %s (only beginning shown)'\
    5089             %(str(longitudes[:10]) + '...',
    5090               str(ref_longitudes[:10]) + '...')       
    5091         assert allclose(longitudes, ref_longitudes), msg       
    5092        
    5093         msg = 'Latitudes specified in ordering file do not match those found in mux files. '
    5094         msg += 'I got %s instead of %s (only beginning shown)'\
    5095             %(str(latitudes[:10]) + '...',
    5096               str(ref_latitudes[:10]) + '...')               
    5097         assert allclose(latitudes, ref_latitudes), msg
    5098        
    5099     else:
    5100         latitudes=latitudes_urs
    5101         longitudes=longitudes_urs
    5102         permutation = range(latitudes.shape[0])               
    5103        
    5104        
     5100
     5101        #ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]               
     5102        #ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]       
     5103        #
     5104        #msg = 'Longitudes specified in ordering file do not match those found in mux files'
     5105        #msg += 'I got %s instead of %s (only beginning shown)'\
     5106        #    %(str(longitudes[:10]) + '...',
     5107        #      str(ref_longitudes[:10]) + '...')       
     5108        #assert allclose(longitudes, ref_longitudes), msg       
     5109        #
     5110        #msg = 'Latitudes specified in ordering file do not match those found in mux files. '
     5111        #msg += 'I got %s instead of %s (only beginning shown)'\
     5112        #    %(str(latitudes[:10]) + '...',
     5113        #      str(ref_latitudes[:10]) + '...')               
     5114        #assert allclose(latitudes, ref_latitudes), msg
     5115       
     5116               
    51055117       
    51065118    # Store timeseries in NetCDF STS file   
     
    51085120    assert len(latitudes>0),msg
    51095121
    5110     number_of_points = latitudes.shape[0]
    5111     number_of_times = times.shape[0]
    5112     number_of_latitudes = latitudes.shape[0]
    5113     number_of_longitudes = longitudes.shape[0]
     5122    number_of_points = latitudes.shape[0]      # Number of stations retrieved
     5123    number_of_times = times.shape[0]           # Number of timesteps
     5124    number_of_latitudes = latitudes.shape[0]   # Number latitudes
     5125    number_of_longitudes = longitudes.shape[0] # Number longitudes
    51145126
    51155127    # NetCDF file definition
     
    51205132   
    51215133    # Create new file
    5122     #starttime = times[0]
    51235134    sts = Write_sts()
    51245135    sts.store_header(outfile,
     
    51315142    # Store
    51325143    from anuga.coordinate_transforms.redfearn import redfearn
    5133     x = zeros(number_of_points, Float)  #Easting
    5134     y = zeros(number_of_points, Float)  #Northing
     5144    x = zeros(number_of_points, Float)  # Easting
     5145    y = zeros(number_of_points, Float)  # Northing
    51355146
    51365147    # Check zone boundaries
     
    51655176    if verbose: print 'Converting quantities'
    51665177    for j in range(len(times)):
    5167         for i, index in enumerate(permutation):
    5168             w = zscale*mux['HA'][index,j] + mean_stage
    5169             h=w-elevation[index]
     5178        for i in range(number_of_points):
     5179            w = zscale*mux['HA'][i,j] + mean_stage
     5180            h=w-elevation[i]
    51705181            stage[j,i] = w
    51715182
    5172             xmomentum[j,i] = mux['UA'][index,j]*h
    5173             ymomentum[j,i] = mux['VA'][index,j]*h
     5183            xmomentum[j,i] = mux['UA'][i,j]*h
     5184            ymomentum[j,i] = mux['VA'][i,j]*h
    51745185
    51755186    outfile.close()
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5531 r5537  
    61076107        return base_name, files
    61086108
    6109     def test_read_mux2_pyI(self):
    6110         """Constant stage,momentum at each gauge
     6109    def test_urs2sts_read_mux2_pyI(self):
     6110        """test_urs2sts_read_mux2_pyI(self):
     6111        Constant stage,momentum at each gauge
    61116112        """
    61126113        tide = 1
     
    61596160        assert allclose(yvelocity,va)
    61606161
    6161     def test_read_mux2_pyII(self):
     6162    def test_urs2sts_read_mux2_pyII(self):
    61626163        """Spatially varing stage
    61636164        """
     
    62176218        assert allclose(yvelocity,va)
    62186219
    6219     def test_read_mux2_pyIII(self):
     6220    def test_urs2sts_read_mux2_pyIII(self):
    62206221        """Varying start and finsh times
    62216222        """
     
    65876588        permutation = [3,0,2]
    65886589
    6589         # Do it wrongly at first and check that exception is being raised
     6590        _, ordering_filename = tempfile.mkstemp('')
     6591        order_fid = open(ordering_filename, 'w') 
     6592        order_fid.write('index, longitude, latitude\n')
     6593        for index in permutation:
     6594            order_fid.write('%d, %f, %f\n' %(index,
     6595                                             lat_long_points[index][1],
     6596                                             lat_long_points[index][0]))
     6597        order_fid.close()
     6598       
     6599           
     6600
     6601                                               
     6602        # Call urs2sts with multiple mux files
     6603        urs2sts([base_nameI, base_nameII],
     6604                basename_out=base_nameI,
     6605                ordering_filename=ordering_filename,
     6606                weights=[1.0,1.0],
     6607                mean_stage=tide,
     6608                verbose=False)
     6609
     6610        # now I want to check the sts file ...
     6611        sts_file = base_nameI + '.sts'
     6612
     6613        #Let's interrogate the sts file
     6614        # Note, the sts info is not gridded.  It is point data.
     6615        fid = NetCDFFile(sts_file)
     6616
     6617        # Make x and y absolute
     6618        x = fid.variables['x'][:]
     6619        y = fid.variables['y'][:]
     6620
     6621        geo_reference = Geo_reference(NetCDFObject=fid)
     6622        points = geo_reference.get_absolute(map(None, x, y))
     6623        points = ensure_numeric(points)
     6624
     6625        x = points[:,0]
     6626        y = points[:,1]
     6627
     6628        for i, index in enumerate(permutation):
     6629            # Check that STS points are stored in the correct order
     6630           
     6631            # Work out the UTM coordinates sts point i
     6632            zone, e, n = redfearn(lat_long_points[index][0],
     6633                                  lat_long_points[index][1])             
     6634           
     6635            assert allclose([x[i],y[i]], [e,n])
     6636           
     6637                       
     6638        # Check the time vector
     6639        times = fid.variables['time'][:]
     6640
     6641        times_actual = []
     6642        for i in range(time_step_count):
     6643            times_actual.append(time_step * i)
     6644
     6645        assert allclose(ensure_numeric(times),
     6646                        ensure_numeric(times_actual))
     6647                       
     6648
     6649        # Check sts values
     6650        stage = fid.variables['stage'][:]
     6651        xmomentum = fid.variables['xmomentum'][:]
     6652        ymomentum = fid.variables['ymomentum'][:]
     6653        elevation = fid.variables['elevation'][:]
     6654
     6655        # Set original data used to write mux file to be zero when gauges are
     6656        # not recdoring
     6657       
     6658        ha[0][0]=0.0
     6659        ha[0][time_step_count-1]=0.0
     6660        ha[2][0]=0.0
     6661        ua[0][0]=0.0
     6662        ua[0][time_step_count-1]=0.0
     6663        ua[2][0]=0.0
     6664        va[0][0]=0.0
     6665        va[0][time_step_count-1]=0.0
     6666        va[2][0]=0.0;
     6667
     6668       
     6669        # The stage stored in the .sts file should be the sum of the stage
     6670        # in the two mux2 files because both have weights = 1. In this case
     6671        # the mux2 files are the same so stage == 2.0 * ha
     6672        #print 2.0*transpose(ha) - stage
     6673       
     6674        ha_permutation = take(ha, permutation)
     6675        ua_permutation = take(ua, permutation)         
     6676        va_permutation = take(va, permutation)                 
     6677        gauge_depth_permutation = take(gauge_depth, permutation)                         
     6678
     6679       
     6680        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
     6681
     6682        #Check the momentums - ua
     6683        #momentum = velocity*(stage-elevation)
     6684        # elevation = - depth
     6685        #momentum = velocity_ua *(stage+depth)
     6686
     6687        depth=zeros((len(lat_long_points),time_step_count),Float)
     6688        for i in range(len(lat_long_points)):
     6689            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
     6690            #2.0*ha necessary because using two files with weights=1 are used
     6691           
     6692        depth_permutation = take(depth, permutation)                     
     6693       
     6694
     6695        # The xmomentum stored in the .sts file should be the sum of the ua
     6696        # in the two mux2 files multiplied by the depth.
     6697        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
     6698
     6699        #Check the momentums - va
     6700        #momentum = velocity*(stage-elevation)
     6701        # elevation = - depth
     6702        #momentum = velocity_va *(stage+depth)
     6703
     6704        # The ymomentum stored in the .sts file should be the sum of the va
     6705        # in the two mux2 files multiplied by the depth.
     6706        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
     6707
     6708        # check the elevation values.
     6709        # -ve since urs measures depth, sww meshers height,
     6710        assert allclose(-gauge_depth_permutation, elevation)  #Meters
     6711
     6712        fid.close()
     6713        self.delete_mux(filesI)
     6714        self.delete_mux(filesII)
     6715        os.remove(sts_file)
     6716       
     6717
     6718       
     6719       
     6720       
     6721    def test_urs2sts_ordering_exception(self):
     6722        """Test that inconsistent lats and lons in ordering file are caught.
     6723        """
     6724       
     6725        tide=0
     6726        time_step_count = 3
     6727        time_step = 2
     6728        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     6729        n=len(lat_long_points)
     6730        first_tstep=ones(n,Int)
     6731        first_tstep[0]+=1
     6732        first_tstep[2]+=1
     6733        last_tstep=(time_step_count)*ones(n,Int)
     6734        last_tstep[0]-=1
     6735
     6736        gauge_depth=20*ones(n,Float)
     6737        ha=2*ones((n,time_step_count),Float)
     6738        ha[0]=arange(0,time_step_count)
     6739        ha[1]=arange(time_step_count,2*time_step_count)
     6740        ha[2]=arange(2*time_step_count,3*time_step_count)
     6741        ha[3]=arange(3*time_step_count,4*time_step_count)
     6742        ua=5*ones((n,time_step_count),Float)
     6743        va=-10*ones((n,time_step_count),Float)
     6744
     6745        # Create two identical mux files to be combined by urs2sts
     6746        base_nameI, filesI = self.write_mux2(lat_long_points,
     6747                                             time_step_count, time_step,
     6748                                             first_tstep, last_tstep,
     6749                                             depth=gauge_depth,
     6750                                             ha=ha,
     6751                                             ua=ua,
     6752                                             va=va)
     6753
     6754        base_nameII, filesII = self.write_mux2(lat_long_points,
     6755                                               time_step_count, time_step,
     6756                                               first_tstep, last_tstep,
     6757                                               depth=gauge_depth,
     6758                                               ha=ha,
     6759                                               ua=ua,
     6760                                               va=va)
     6761
     6762                                               
     6763        # Create ordering file
     6764        permutation = [3,0,2]
     6765
     6766        # Do it wrongly and check that exception is being raised
    65906767        _, ordering_filename = tempfile.mkstemp('')
    65916768        order_fid = open(ordering_filename, 'w') 
     
    66136790
    66146791       
    6615            
    6616         # Then do it right
    6617         _, ordering_filename = tempfile.mkstemp('')
    6618         order_fid = open(ordering_filename, 'w') 
    6619         order_fid.write('index, longitude, latitude\n')
    6620         for index in permutation:
    6621             order_fid.write('%d, %f, %f\n' %(index,
    6622                                              lat_long_points[index][1],
    6623                                              lat_long_points[index][0]))
    6624         order_fid.close()
    6625        
    6626            
    6627 
    6628                                                
    6629         # Call urs2sts with multiple mux files
    6630         urs2sts([base_nameI, base_nameII],
    6631                 basename_out=base_nameI,
    6632                 ordering_filename=ordering_filename,
    6633                 weights=[1.0,1.0],
    6634                 mean_stage=tide,
    6635                 verbose=False)
    6636 
    6637         # now I want to check the sts file ...
    6638         sts_file = base_nameI + '.sts'
    6639 
    6640         #Let's interrogate the sts file
    6641         # Note, the sts info is not gridded.  It is point data.
    6642         fid = NetCDFFile(sts_file)
    6643 
    6644         # Make x and y absolute
    6645         x = fid.variables['x'][:]
    6646         y = fid.variables['y'][:]
    6647 
    6648         geo_reference = Geo_reference(NetCDFObject=fid)
    6649         points = geo_reference.get_absolute(map(None, x, y))
    6650         points = ensure_numeric(points)
    6651 
    6652         x = points[:,0]
    6653         y = points[:,1]
    6654 
    6655         for i, index in enumerate(permutation):
    6656             # Check that STS points are stored in the correct order
    6657            
    6658             # Work out the UTM coordinates sts point i
    6659             zone, e, n = redfearn(lat_long_points[index][0],
    6660                                   lat_long_points[index][1])             
    6661            
    6662             assert allclose([x[i],y[i]], [e,n])
    6663            
    6664                        
    6665         # Check the time vector
    6666         times = fid.variables['time'][:]
    6667 
    6668         times_actual = []
    6669         for i in range(time_step_count):
    6670             times_actual.append(time_step * i)
    6671 
    6672         assert allclose(ensure_numeric(times),
    6673                         ensure_numeric(times_actual))
    6674                        
    6675 
    6676         # Check sts values
    6677         stage = fid.variables['stage'][:]
    6678         xmomentum = fid.variables['xmomentum'][:]
    6679         ymomentum = fid.variables['ymomentum'][:]
    6680         elevation = fid.variables['elevation'][:]
    6681 
    6682         # Set original data used to write mux file to be zero when gauges are
    6683         # not recdoring
    6684        
    6685         ha[0][0]=0.0
    6686         ha[0][time_step_count-1]=0.0
    6687         ha[2][0]=0.0
    6688         ua[0][0]=0.0
    6689         ua[0][time_step_count-1]=0.0
    6690         ua[2][0]=0.0
    6691         va[0][0]=0.0
    6692         va[0][time_step_count-1]=0.0
    6693         va[2][0]=0.0;
    6694 
    6695        
    6696         # The stage stored in the .sts file should be the sum of the stage
    6697         # in the two mux2 files because both have weights = 1. In this case
    6698         # the mux2 files are the same so stage == 2.0 * ha
    6699         #print 2.0*transpose(ha) - stage
    6700        
    6701         ha_permutation = take(ha, permutation)
    6702         ua_permutation = take(ua, permutation)         
    6703         va_permutation = take(va, permutation)                 
    6704         gauge_depth_permutation = take(gauge_depth, permutation)                         
    6705 
    6706        
    6707         assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
    6708 
    6709         #Check the momentums - ua
    6710         #momentum = velocity*(stage-elevation)
    6711         # elevation = - depth
    6712         #momentum = velocity_ua *(stage+depth)
    6713 
    6714         depth=zeros((len(lat_long_points),time_step_count),Float)
    6715         for i in range(len(lat_long_points)):
    6716             depth[i]=gauge_depth[i]+tide+2.0*ha[i]
    6717             #2.0*ha necessary because using two files with weights=1 are used
    6718            
    6719         depth_permutation = take(depth, permutation)                     
    6720        
    6721 
    6722         # The xmomentum stored in the .sts file should be the sum of the ua
    6723         # in the two mux2 files multiplied by the depth.
    6724         assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
    6725 
    6726         #Check the momentums - va
    6727         #momentum = velocity*(stage-elevation)
    6728         # elevation = - depth
    6729         #momentum = velocity_va *(stage+depth)
    6730 
    6731         # The ymomentum stored in the .sts file should be the sum of the va
    6732         # in the two mux2 files multiplied by the depth.
    6733         assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
    6734 
    6735         # check the elevation values.
    6736         # -ve since urs measures depth, sww meshers height,
    6737         assert allclose(-gauge_depth_permutation, elevation)  #Meters
    6738 
    6739         fid.close()
    6740         self.delete_mux(filesI)
    6741         self.delete_mux(filesII)
    6742         os.remove(sts_file)
    6743        
    6744 
    6745     def test_file_boundary_stsI(self):
     6792       
     6793       
     6794    def test_urs2sts_file_boundary_stsI(self):
    67466795        from anuga.shallow_water import Domain
    67476796        from anuga.shallow_water import Reflective_boundary
     
    68186867        os.remove(meshname)
    68196868
    6820     def test_file_boundary_stsII(self):
     6869    def test_urs2sts_file_boundary_stsII(self):
    68216870        """ mux2 file has points not included in boundary
    68226871         mux2 gauges are not stored with the same order as they are
     
    69016950        os.remove(meshname)
    69026951
    6903     def test_file_boundary_stsIII(self):
     6952    def test_urs2sts_file_boundary_stsIII_ordering(self):
    69046953        """Read correct points from ordering file
    69056954        """
     
    88918940if __name__ == "__main__":
    88928941
    8893     suite = unittest.makeSuite(Test_Data_Manager,'test')
    8894     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts0')
    8895     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
     8942    #suite = unittest.makeSuite(Test_Data_Manager,'test')
     8943    suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_read_mux2_pyI')
     8944    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts')
    88968945    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    88978946    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
     
    89068955    else:
    89078956        pass
    8908     runner = unittest.TextTestRunner() # verbosity=2)
     8957    runner = unittest.TextTestRunner(verbosity=2)
    89098958    runner.run(suite)
    89108959
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5531 r5537  
    2828long getNumData(const int *fros, const int *lros, const int nsta);
    2929char isdata(float);
    30 float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose);
     30float** _read_mux2(int numSrc,
     31                   char **muxFileNameArray,
     32                   float *weights,
     33                   double *params,
     34                   int number_of_selected_stations,
     35                   int *permutation,
     36                   int verbose);
    3137
    3238PyObject *read_mux2(PyObject *self, PyObject *args){
     
    3440
    3541    Python call:
    36     read_mux2(numSrc,filenames,weights,file_params,verbose)
     42    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    3743
    3844    NOTE:
     
    4349    PyArrayObject *pyweights;
    4450    PyArrayObject *file_params;
     51    PyArrayObject *permutation;   
    4552    PyArrayObject *pydata;
    4653    PyObject *fname;
     
    5360    int verbose;
    5461    int nsta0;
     62    int number_of_selected_stations;   
    5563    int nt;
    5664    double dt;
     
    6674
    6775    // Convert Python arguments to C
    68     if (!PyArg_ParseTuple(args, "iOOOi",
    69         &numSrc, &filenames, &pyweights, &file_params, &verbose))
     76    if (!PyArg_ParseTuple(args, "iOOOOi",
     77                          &numSrc, &filenames, &pyweights, &file_params, &permutation, &verbose))
    7078    {
    7179            PyErr_SetString(PyExc_RuntimeError,
     
    106114    }
    107115
    108     //printf("Checkpoint 2 for urs2sts_ext.c\n");   
     116    printf("Checkpoint 2 for urs2sts_ext.c\n");   
    109117    //printf("numSrc = %d\n", numSrc);
    110118    for (i = 0; i < numSrc; i++)
     
    128136        if (muxFileNameArray[i] == NULL)
    129137        {
    130             printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
     138            PyErr_SetString(PyExc_ValueError,
     139                            "ERROR: Memory for muxFileNameArray could not be allocated.\n");
    131140            return NULL;
    132141        }
     
    137146    {
    138147        PyErr_SetString(PyExc_ValueError,
    139             "file_params must be one-dimensional and of type double");
     148                        "file_params must be one-dimensional and of type double");
    140149        return NULL;
    141150    }
    142151
    143     //printf("Checkpoint 4 for urs2sts_ext.c\n");       
     152    printf("Checkpoint 4 for urs2sts_ext.c\n");       
    144153    // Create array for weights which are passed to read_mux2
    145154    weights = (float*) malloc(numSrc*sizeof(float));
     
    149158    }
    150159   
    151     //printf("Checkpoint 5 for urs2sts_ext.c\n");           
     160   
     161    number_of_selected_stations = (int) permutation->dimensions[0];
     162    printf("Checkpoint 5 for urs2sts_ext.c\n");       
    152163    // Read in mux2 data from file
    153164    cdata = _read_mux2(numSrc,
    154165                       muxFileNameArray,
    155166                       weights,
    156                        (double*)file_params->data,
     167                       (double*)file_params->data,
     168                       number_of_selected_stations, // Desired number of stations
     169                       (int*) permutation->data, // Ordering of selected stations
    157170                       verbose);
    158171
     172    if (!cdata)
     173    {
     174        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
     175        return NULL;
     176    }       
     177                       
     178                       
    159179    //printf("Checkpoint 6 for urs2sts_ext.c\n");               
    160180    // Allocate space for return vector
     
    163183    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
    164184
    165     //printf("Checkpoint 7 for urs2sts_ext.c\n");                   
     185    printf("Checkpoint 7 for urs2sts_ext.c\n");                   
     186   
    166187    // Find min and max start times of all gauges
    167188    start_tstep = nt + 1;
    168189    finish_tstep = -1;
    169     for (i = 0; i < nsta0; i++)
    170     {
     190    for (i = 0; i < number_of_selected_stations; i++)
     191    {
     192        printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
     193        printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]); 
     194       
    171195        if ((int)cdata[i][nt + 3] < start_tstep)
    172196        {
     
    179203    }
    180204
    181     //printf("Checkpoint 8 for urs2sts_ext.c\n");                       
     205    printf("Checkpoint 8 for urs2sts_ext.c\n");                       
    182206    if ((start_tstep > nt) | (finish_tstep < 0))
    183207    {
    184         printf("ERROR: Gauge data has incorrect start and finsh times\n");
     208        printf("ERROR: Gauge data has incorrect start and finish times:\n");
     209        printf("   start_tstep = %d, max_number_of_steps = %d\n", start_tstep, nt);
     210        printf("   finish_tstep = %d, min_number_of_steps = %d\n", finish_tstep, 0);   
     211        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
    185212        return NULL;
    186213    }
     
    188215    if (start_tstep >= finish_tstep)
    189216    {
    190         printf("ERROR: Gauge data has non-postive_length\n");
     217        PyErr_SetString(PyExc_ValueError, "ERROR: Gauge data has non-postive_length");
    191218        return NULL;
    192219    }
    193220
    194221    num_ts = finish_tstep - start_tstep + 1;
    195     dimensions[0] = nsta0;
     222    dimensions[0] = number_of_selected_stations;
    196223    dimensions[1] = num_ts + POFFSET;
    197224    pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    198225    if(pydata == NULL)
    199226    {
    200         printf("ERROR: Memory for pydata array could not be allocated.\n");
    201         return NULL;
    202     }
    203 
    204     //printf("Checkpoint 9 for urs2sts_ext.c\n");                           
     227        PyErr_SetString(PyExc_ValueError, "ERROR: Memory for pydata array could not be allocated.");
     228        return NULL;
     229    }
     230
     231    printf("Checkpoint 9 for urs2sts_ext.c\n");                           
    205232    // Each gauge begins and ends recording at different times. When a gauge is
    206233    // not recording but at least one other gauge is. Pad the non-recording gauge
    207234    // array with zeros.
    208     for (i = 0; i < nsta0; i++)
     235    for (i = 0; i < number_of_selected_stations; i++)
    209236    {
    210237        time = 0;
     
    232259    }
    233260
    234     //printf("Checkpoint 10 for urs2sts_ext.c\n");                               
     261    printf("Checkpoint 10 for urs2sts_ext.c\n");                               
    235262    free(weights);
    236263   
     
    239266    free(muxFileNameArray);
    240267   
    241     for (i = 0; i < nsta0; ++i)
     268    for (i = 0; i < number_of_selected_stations; ++i)
    242269    {
    243270        free(cdata[i]);
     
    245272    free(cdata);
    246273
    247     //printf("Checkpoint 11 for urs2sts_ext.c\n");                                   
     274    printf("Checkpoint 11 for urs2sts_ext.c\n");                                   
    248275    return  PyArray_Return(pydata);
    249276}
     
    253280                   float *weights,
    254281                   double *params,
     282                   int number_of_selected_stations,
     283                   int *permutation,
    255284                   int verbose)
    256285{
     
    299328    }
    300329
    301     /* loop over all sources, read headers and check compatibility */
     330    printf("P1\n");
     331    /* Loop over all sources, read headers and check compatibility */
    302332    for (isrc = 0; isrc < numSrc; isrc++)
    303333    {
    304334        muxFileName = muxFileNameArray[isrc];
    305335
    306                 /* open the mux file */
     336        /* open the mux file */
    307337        if((fp = fopen(muxFileName, "r")) == NULL)
    308338        {
    309339            fprintf(stderr, "cannot open file %s\n", muxFileName);
    310             exit(-1)
     340            return NULL
    311341        }
    312342       
     
    329359            if(nsta != nsta0)
    330360            {
    331                 fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
     361                fprintf(stderr,"%s has different number of stations to %s\n",
     362                        muxFileName,
     363                        muxFileNameArray[0]);
    332364                fclose(fp);
    333                 exit(-1);   
     365                return NULL;   
    334366            }
    335367
     
    340372                if (mytgs[ista].dt != mytgs0[ista].dt)
    341373                {
    342                     fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
     374                    fprintf(stderr,"%s has different sampling rate to %s\n",
     375                            muxFileName,
     376                            muxFileNameArray[0]);
    343377                    fclose(fp);
    344                     exit(-1);                
     378                    return NULL;                   
    345379                }   
    346380                if (mytgs[ista].nt != mytgs0[ista].nt)
    347381                {
    348                     fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
     382                    fprintf(stderr,"%s has different series length to %s\n",
     383                            muxFileName,
     384                            muxFileNameArray[0]);
    349385                    fclose(fp);
    350                     exit(-1);                
     386                    return NULL;                   
    351387                }
    352388            }
    353389        }
    354390
    355         /* read the start and stop times for this source */
     391        /* Read the start and stop times for this source */
    356392        fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    357393        fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    358394
    359         /* compute the size of the data block for this source */
     395        /* Compute the size of the data block for this source */
    360396        numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
    361397
    362         /* Burbidge: Added a sanity check here */
     398        /* Sanity check */
    363399        if (numData < 0)
    364400        {
    365401            fprintf(stderr,"Size of data block appears to be negative!\n");
    366             exit(-1);
     402            return NULL;           
    367403        }
    368404
     
    370406    }
    371407
     408    printf("P2\n");   
    372409    params[0] = (double)nsta0;
    373410    params[1] = (double)mytgs0[0].dt;
    374411    params[2] = (double)mytgs0[0].nt;
    375412
    376     /* make array(s) to hold the demuxed data */
    377     //FIXME: This can be reduced to only contain stations given in permutation file
    378     sts_data = (float**)malloc(nsta0*sizeof(float *));
     413    // Apply rule that an empty permutation file means 'take all stations'
     414    // We can change this later by passing in None instead of the empty permutation.
     415   
     416    printf("number_of_selected_stations B4 = %d\n", number_of_selected_stations);   
     417    if (number_of_selected_stations == 0)
     418    {
     419        number_of_selected_stations = nsta0; 
     420     
     421        printf("Creating identity permutation vector of length = %d\n", nsta0);
     422        // Create the Identity permutation vector
     423        permutation = (int *) malloc(number_of_selected_stations*sizeof(int));
     424        for (i = 0; i < number_of_selected_stations; i++)
     425        {
     426          permutation[i] = i; 
     427        }
     428     
     429    }
     430   
     431    printf("P3\n");       
     432    printf("Number of selected stations = %d\n", number_of_selected_stations);
     433    /* Make array(s) to hold demuxed data for stations given in the permutation file */
     434    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
    379435    if (sts_data == NULL)
    380436    {
    381437        printf("ERROR: Memory for sts_data could not be allocated.\n");
    382         exit(-1);
    383     }
    384 
    385     len_sts_data = mytgs0[0].nt + POFFSET;
    386     for (ista = 0; ista < nsta0; ista++)
    387     {
     438        return NULL;
     439    }
     440
     441    printf("P4\n");           
     442    // For each selected station, allocate space for its data
     443    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries (I think)
     444    for (i = 0; i < number_of_selected_stations; i++)
     445    {
     446        printf("P4.1, i=%d\n", i);               
     447   
    388448        // Initialise sts_data to zero
    389         sts_data[ista] = (float*)calloc(len_sts_data, sizeof(float));
    390         if (sts_data[ista] == NULL)
     449        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
     450        if (sts_data[i] == NULL)
    391451        {
    392452            printf("ERROR: Memory for sts_data could not be allocated.\n");
    393             exit(-1);
    394         }
    395 
    396         sts_data[ista][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
    397         sts_data[ista][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
    398         sts_data[ista][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
    399         sts_data[ista][mytgs0[0].nt + 3] = (float)fros[ista];
    400         sts_data[ista][mytgs0[0].nt + 4] = (float)lros[ista];
    401     }
    402 
     453            return NULL;
     454        }
     455
     456       
     457        printf("P4.2, i=%d\n", i);                     
     458        ista = permutation[i]; // Get global index into mux data
     459       
     460
     461        printf("P4.3, i=%d\n", i);                     
     462        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
     463        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
     464        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
     465        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
     466        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
     467    }
     468
     469    printf("P5\n");               
    403470    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
    404471
     
    415482            //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
    416483            fprintf(stderr, "cannot open file %s\n", muxFileName);
    417             exit(-1); 
     484            return NULL;                           
    418485        }
    419486
    420487        if (verbose){
    421             printf("Reading mux file %s\n",muxFileName);
     488            printf("Reading mux file %s\n", muxFileName);
    422489        }
    423490
     
    430497        fclose(fp);
    431498
    432         /* loop over stations */
    433         /* This is where we should only take stations with indices in the permutation array
    434            (and in that order) */
    435        
    436         //for i in range(len(permution)):
    437         //     ista = permutation[i]
     499        // loop over stations present in the permutation array
    438500        //     use ista with mux data
    439501        //     use i with the processed data to be returned         
    440         i=0;
    441502       
    442         for(ista = 0; ista < nsta0; ista++)
     503        for(i = 0; i < number_of_selected_stations; i++)
    443504        {               
     505       
     506            ista = permutation[i]; // Get global index into mux data   
     507           
    444508            /* fill the data0 array from the mux file, and weight it */
    445             fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros + isrc*nsta0, lros + isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
     509            fillDataArray(ista,
     510                          nsta0,
     511                          mytgs0[ista].nt,
     512                          mytgs0[ista].ig,
     513                          fros + isrc*nsta0,
     514                          lros + isrc*nsta0,
     515                          temp_sts_data,
     516                          &istart,
     517                          &istop,
     518                          muxData);
    446519
    447520            /* weight appropriately and add */
    448521            for(k = 0; k < mytgs0[ista].nt; k++)
    449522            {
    450                 if((isdata(sts_data[ista][k])) && isdata(temp_sts_data[k]))
     523                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
    451524                {
    452                     sts_data[ista][k] += temp_sts_data[k] * weights[isrc];
     525                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
    453526                    //printf("%d,%d,%f\n",ista,k,sts_data[ista][k]);
    454527                }
    455528                else
    456529                {
    457                     sts_data[ista][k] = NODATA;
     530                    sts_data[i][k] = NODATA;
    458531                }
    459532            }
     
    461534    }
    462535
     536    printf("P6\n");                   
    463537    free(muxData);
    464538    free(temp_sts_data);
Note: See TracChangeset for help on using the changeset viewer.