Changeset 3826


Ignore:
Timestamp:
Oct 19, 2006, 5:34:56 PM (17 years ago)
Author:
nick
Message:

updates to data_manager.py urs2sww

File:
1 edited

Legend:

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

    r3824 r3826  
    42304230    (points_num,)= unpack('i',f.read(4))
    42314231
     4232    #print 'points_num', points_num
     4233    #import sys; sys.exit()
    42324234    # nt, int - Number of time steps
    42334235    (time_step_count,)= unpack('i',f.read(4))
     
    42484250   
    42494251    lonlatdep = p_array.array('f')
     4252    points_num = 2914
    42504253    lonlatdep.read(f, columns * points_num)
    42514254    lonlatdep = array(lonlatdep, typecode=Float)   
    4252     lonlatdep = reshape(lonlatdep, ( points_num, columns))
     4255    lonlatdep = reshape(lonlatdep, (points_num, columns))
    42534256   
    42544257    #print "lonlatdep", lonlatdep
     
    42584261    lon_sorted.sort()
    42594262   
    4260     if not lon == lon_sorted:
     4263#    if not lon == lon_sorted:
     4264    if not allclose(lon, lon_sorted):
    42614265        msg = "Longitudes in mux file are not in ascending order"
    42624266        raise IOError, msg
    42634267    lat_sorted = lat[:]
    42644268    lat_sorted.sort()   
    4265     if not lat == lat_sorted:
     4269#    if not lat == lat_sorted:
     4270    if not allclose(lat, lat_sorted):
    42664271        msg = "Latitudes in mux file are not in ascending order"
    42674272   
    42684273    nc_file = Write_nc(quantity,
    42694274                       file_out,
    4270                         time_step_count,
    4271                         time_step,
    4272                         lon,
    4273                         lat)
     4275                       time_step_count,
     4276                       time_step,
     4277                       lon,
     4278                       lat)
    42744279
    42754280    for i in range(time_step_count):
     
    43484353    LAT = 1
    43494354    QUANTITY = 2
    4350     points_num = len(long_lat_dep)
    4351     lat = [long_lat_dep[0][LAT]]
    4352     this_rows_long = long_lat_dep[0][LONG]
    4353     i = 1 # Index of long_lat_dep
    43544355   
    4355     #Working out the lat's
    4356     while long_lat_dep[i][LONG] == this_rows_long and i < points_num:
    4357             lat.append(long_lat_dep[i][LAT])
    4358             i += 1
    4359            
    4360     lats_per_long = i       
    4361     long = [long_lat_dep[0][LONG]]
    4362     #Working out the longs
    4363     while i < points_num:
     4356#    print 'long_lat_dep', type(long_lat_dep), long_lat_dep.shape
     4357#    print 'long_lat_dep', long_lat_dep
     4358   
     4359    num_points = long_lat_dep.shape[0]
     4360    this_rows_long = long_lat_dep[0,LONG]
     4361   
     4362    # Count the length of unique latitudes
     4363    i = 0
     4364    while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
     4365        i += 1
     4366         
     4367    lat = long_lat_dep[:i, LAT]       
     4368    long = long_lat_dep[::i, LONG]     
     4369    lenlong = len(long)
     4370    lenlat = len(lat)
     4371#    print 'len lat', lat, len(lat)
     4372#    print 'len long', long, len(long)
     4373         
     4374    msg = 'Input data is not gridded'     
     4375    assert num_points % lenlat == 0, msg
     4376    assert num_points % lenlong == 0, msg
     4377         
     4378    # Test that data is gridded       
     4379    for i in range(lenlong):
    43644380        msg = 'Data is not gridded.  It must be for this operation'
    4365         assert long_lat_dep[i][LAT] == lat[i%lats_per_long], msg
    4366         if i%lats_per_long == 0:
    4367             long.append(long_lat_dep[i][LONG])
    4368         i += 1
     4381        first = i*lenlat
     4382        last = first + lenlat
     4383               
     4384        assert allclose(long_lat_dep[first:last,LAT], lat), msg
     4385        assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
    43694386   
    4370     msg = 'Our of range latitudes/longitudes'
     4387   
     4388#    print 'range long', min(long), max(long)
     4389#    print 'range lat', min(lat), max(lat)
     4390#    print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0])
     4391#    print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1])
     4392   
     4393   
     4394   
     4395    msg = 'Out of range latitudes/longitudes'
    43714396    for l in lat:assert -90 < l < 90 , msg
    43724397    for l in long:assert -180 < l < 180 , msg
    43734398
    4374     #changing quantity from lat being the fastest varying dimension to
     4399    # Changing quantity from lat being the fastest varying dimension to
    43754400    # long being the fastest varying dimension
    43764401    # FIXME - make this faster/do this a better way
    43774402    # use numeric transpose, after reshaping the quantity vector
    4378     quantity = zeros(len(long_lat_dep), Float)
    4379     lenlong = len(long)
    4380     lenlat = len(lat)
     4403#    quantity = zeros(len(long_lat_dep), Float)
     4404    quantity = zeros(num_points, Float)
     4405   
     4406#    print 'num',num_points
    43814407    for lat_i, _ in enumerate(lat):
    43824408        for long_i, _ in enumerate(long):
    43834409            q_index = lat_i*lenlong+long_i
    43844410            lld_index = long_i*lenlat+lat_i
    4385             quantity[q_index] = long_lat_dep[lld_index][QUANTITY]
     4411#            print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index
     4412            temp = long_lat_dep[lld_index, QUANTITY]
     4413            quantity[q_index] = temp
    43864414           
    43874415    return long, lat, quantity
Note: See TracChangeset for help on using the changeset viewer.