Changeset 5580


Ignore:
Timestamp:
Jul 25, 2008, 6:35:51 PM (11 years ago)
Author:
sexton
Message:

update script to write time series to file from sts file (for Perth)

File:
1 edited

Legend:

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

    r5576 r5580  
    11import os
    22from os import sep
    3 #import create_boundary
    43import project
    54try:
     
    1817    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    1918    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    20     print 'heloo', len(x), len(y)
    2119    points=transpose(asarray([x.tolist(),y.tolist()]))
    2220    time=fid.variables['time'][:]
     
    3230        #plot(polygon[:,0],polygon[:,1])
    3331        #show()
    34        
     32
     33    basename='sts_gauge'
    3534    quantity_names=['stage','xmomentum','ymomentum']
    36    
    37     fid_sts = open(project.boundaries_dir+sep+'sts_gauge.csv', 'w')
    38     s = 'time, stage, xmomentum, ymomentum \n'
    39     fid_sts.write(s)
    4035    quantities = {}
    4136    for i, name in enumerate(quantity_names):
    4237        quantities[name] = fid.variables[name][:]
    43         print len(quantities[name])
    4438        if inside_points is not None:
    4539            quantities[name] = take(quantities[name],indices,1)
    46             for j in range(len(x)):
     40            for j in range(len(x)-1):
     41                locx=int(x[j])
     42                locy=int(y[j])
     43                fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
     44                s = 'time, stage, xmomentum, ymomentum \n'
     45                fid_sts.write(s)
    4746                stage = quantities[name][:,j]
    4847                xmomentum = quantities[name][:,j]
    4948                ymomentum = quantities[name][:,j]
    50                 #print 'hello', len(stage)
    51                 #s = '%.2f, %.2f, %.2f, %.2f\n' %(time[j], stage, xmomentum, ymomentum)
    52                 #fid_sts.write(s)
     49                for k in range(len(time)-1):
     50                    s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
     51                    fid_sts.write(s)
    5352        else:
    5453            msg = 'No gauges found inside polygon'
    5554            raise msg
     55            fid_sts.close()
    5656    fid.close()
    5757
     
    5959
    6060polygon=project.poly_all
    61 points,quantities,elevation,time=get_sts_gauge_data(project.scenario_name,polygon,verbose=False)
    62 from Numeric import size
    63 print 'hello', len(quantities['stage'])
     61points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
     62
    6463
    6564assert len(points)==len(elevation)==len(quantities['stage'][0,:])
    66 
    67 ##def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos):
    68 ##    figname=figname = '%s_%s%s.png' %('tgs/urs_gauge',
    69 ##                                      quantity_name,
    70 ##                                      str(gauge_id))
    71 ##    hold(False)
    72 ##    plot(time/60.0,quantity,label='%s (elevation = %sm)'%('urs',str(elevation)))
    73 ##    title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1]))
    74 ##    xlabel('Time (minutes)')
    75 ##    ylabel(quantity_name)
    76 ##    legend()
    77 ##    print 'saving figure here %s' %figname
    78 ##    savefig(figname)
    79 ##   
    80 ##points=points[::5]
    81 ##elevation=elevation[::5]
    82 ##gauge_filename=project.boundaries_dir+sep+'urs_gauges'
    83 ##d=','
    84 ##if os.path.exists(gauge_filename+ '.txt'):
    85 ##    print 'file with list of gauges already created'
    86 ##else:
    87 ##    fid = open(gauge_filename+ '.txt', 'w')
    88 ##    fid.write('easting,northing,name,elevation\n')
    89 ##    for i in range(len(points)):
    90 ##        fid.write(str(points[i][0])+d+str(points[i][1])+d+str(i)+d+str(elevation[i])+'\n')
    91 ##    fid.close()
    92 
Note: See TracChangeset for help on using the changeset viewer.