Ignore:
Timestamp:
Jul 18, 2008, 4:59:46 PM (14 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.

File:
1 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()
Note: See TracChangeset for help on using the changeset viewer.