Ignore:
Timestamp:
Aug 5, 2008, 1:17:37 PM (14 years ago)
Author:
kristy
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/perth/get_gauges.py

    r5592 r5607  
    1010from Scientific.IO.NetCDF import NetCDFFile
    1111
    12 def get_sts_gauge_data(filename,polygon,verbose=False):
     12def get_sts_gauge_data(filename,verbose=False):
    1313    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    1414        compress,zeros,fabs,take
     
    1717    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    1818    points=transpose(asarray([x.tolist(),y.tolist()]))
    19     time=fid.variables['time'][:]
     19    time=fid.variables['time'][:]+fid.starttime
    2020    elevation=fid.variables['elevation'][:]
    21     indices=inside_polygon(points, polygon, closed=True, verbose=False)
    22     inside_points=take(points,indices,0)
    23     elevation=take(elevation,indices,0)
    24     if verbose:
    25         inside_points_array=ensure_numeric(inside_points)#for plotting only
    26         #plot(points[:,0],points[:,1],'o')
    27         #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
    28         polygon=ensure_numeric(polygon)
    29         #plot(polygon[:,0],polygon[:,1])
    30         #show()
    31 
    32    
     21       
    3322    basename='sts_gauge'
    3423    quantity_names=['stage','xmomentum','ymomentum']
     
    3625    for i, name in enumerate(quantity_names):
    3726        quantities[name] = fid.variables[name][:]
    38 ##        if inside_points is not None:
    39 ##            quantities[name] = take(quantities[name],indices,1)
    40 ##        else:
    41 ##            msg = 'No gauges found inside polygon'
    42 ##            raise msg
    43 ##           
    44     for j in range(len(x)-1):
     27
     28    maxname = 'max_sts_stage.csv'
     29    fid_max = open(project.boundaries_dir+sep+maxname,'w')
     30    s = 'x, y, max_stage \n'
     31    fid_max.write(s)   
     32    for j in range(len(x)):
    4533        locx=int(x[j])
    4634        locy=int(y[j])
     35        stage = quantities['stage'][:,j]
     36        xmomentum = quantities['xmomentum'][:,j]
     37        ymomentum = quantities['ymomentum'][:,j]
     38
     39        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
     40        fid_max.write(s)
     41       
    4742        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
    4843        s = 'time, stage, xmomentum, ymomentum \n'
    4944        fid_sts.write(s)
    50         stage = quantities['stage'][:,j]
    51         xmomentum = quantities['xmomentum'][:,j]
    52         ymomentum = quantities['ymomentum'][:,j]
     45       
    5346        for k in range(len(time)-1):
    54             s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
     47            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
    5548            fid_sts.write(s)
    5649
     
    5952    fid.close()
    6053       
    61     return inside_points,quantities,elevation,time
     54    return quantities,elevation,time
    6255
    63 polygon=project.poly_all
    64 points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
    65 print quantities['stage'][0,:]
    66 print len(points), len(elevation), len(quantities['stage'][0,:])
    67 assert len(points)==len(elevation)==len(quantities['stage'][0,:])
     56quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
     57
     58print len(elevation), len(quantities['stage'][0,:])
Note: See TracChangeset for help on using the changeset viewer.