Changeset 5580 for anuga_work/production/perth/get_gauges.py

Timestamp:

Jul 25, 2008, 6:35:51 PM (16 years ago)

Author:

sexton

Message:

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

File:

• anuga_work/production/perth/get_gauges.py (modified) (4 diffs)

anuga_work/production/perth/get_gauges.py

@@ -1 +1 @@
 import os
 from os import sep
-#import create_boundary
 import project
 try:
@@ -18 +17 @@
     x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
     y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
-    print 'heloo', len(x), len(y)
     points=transpose(asarray([x.tolist(),y.tolist()]))
     time=fid.variables['time'][:]
@@ -32 +30 @@
         #plot(polygon[:,0],polygon[:,1])
         #show()
-        
+
+    basename='sts_gauge'
     quantity_names=['stage','xmomentum','ymomentum']
-    
-    fid_sts = open(project.boundaries_dir+sep+'sts_gauge.csv', 'w')
-    s = 'time, stage, xmomentum, ymomentum \n'
-    fid_sts.write(s)
     quantities = {}
     for i, name in enumerate(quantity_names):
         quantities[name] = fid.variables[name][:]
-        print len(quantities[name])
         if inside_points is not None:
             quantities[name] = take(quantities[name],indices,1)
-            for j in range(len(x)): 
+            for j in range(len(x)-1):
+                locx=int(x[j])
+                locy=int(y[j])
+                fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
+                s = 'time, stage, xmomentum, ymomentum \n'
+                fid_sts.write(s)
                 stage = quantities[name][:,j]
                 xmomentum = quantities[name][:,j]
                 ymomentum = quantities[name][:,j]
-                #print 'hello', len(stage)
-                #s = '%.2f, %.2f, %.2f, %.2f\n' %(time[j], stage, xmomentum, ymomentum)
-                #fid_sts.write(s)
+                for k in range(len(time)-1):
+                    s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
+                    fid_sts.write(s)
         else:
             msg = 'No gauges found inside polygon'
             raise msg
+            fid_sts.close()
     fid.close()
 
@@ -59 +59 @@
 
 polygon=project.poly_all
-points,quantities,elevation,time=get_sts_gauge_data(project.scenario_name,polygon,verbose=False)
-from Numeric import size
-print 'hello', len(quantities['stage'])
+points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
+
 
 assert len(points)==len(elevation)==len(quantities['stage'][0,:])
-
-##def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos):
-##    figname=figname = '%s_%s%s.png' %('tgs/urs_gauge',
-##                                      quantity_name,
-##                                      str(gauge_id))
-##    hold(False)
-##    plot(time/60.0,quantity,label='%s (elevation = %sm)'%('urs',str(elevation)))
-##    title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1]))
-##    xlabel('Time (minutes)')
-##    ylabel(quantity_name)
-##    legend()
-##    print 'saving figure here %s' %figname
-##    savefig(figname)
-##    
-##points=points[::5]
-##elevation=elevation[::5]
-##gauge_filename=project.boundaries_dir+sep+'urs_gauges'
-##d=','
-##if os.path.exists(gauge_filename+ '.txt'):
-##    print 'file with list of gauges already created'
-##else:
-##    fid = open(gauge_filename+ '.txt', 'w')
-##    fid.write('easting,northing,name,elevation\n')
-##    for i in range(len(points)):
-##        fid.write(str(points[i][0])+d+str(points[i][1])+d+str(i)+d+str(elevation[i])+'\n')
-##    fid.close()
-